Skip to content
Snippets Groups Projects
Commit 222475c0 authored by Gervaise Henry's avatar Gervaise Henry :cowboy:
Browse files

Add multi-CCA, fix RunTSNE and FindClusters to use cca.aligned, add Seurat Rename function

parent 318a9ff7
No related merge requests found
......@@ -255,7 +255,7 @@ scPC <- function(sc10x,lx=0.15,hx=3.5,ly=0.75,cc=FALSE,pc=50,hpc=0.75,file="pre.
return(results)
}
scCluster <- function(sc10x,pc.use=10,res.use=0.1,folder="pre.stress"){
scCluster <- function(sc10x,pc.use=10,res.use=0.1,folder="pre.stress",red="pca"){
#Cluster
#Inputs:
......@@ -268,7 +268,7 @@ scCluster <- function(sc10x,pc.use=10,res.use=0.1,folder="pre.stress"){
#Seurat object
#Calculate tSNE
sc10x <- RunTSNE(object=sc10x,dims.use=1:pc.use,do.fast=TRUE)
sc10x <- RunTSNE(object=sc10x,reduction.use=red,dims.use=1:pc.use,do.fast=TRUE)
#Create tSNE grouped by samples
postscript(paste0("./analysis/tSNE/",folder,"/tSNE_Sample.eps"))
......@@ -280,7 +280,7 @@ scCluster <- function(sc10x,pc.use=10,res.use=0.1,folder="pre.stress"){
#Cluster and save resolution
gc()
sc10x <- FindClusters(object=sc10x,reduction.type="pca",dims.use=1:pc.use,resolution=res.use,print.output=0,save.SNN=TRUE,force.recalc=TRUE)
sc10x <- FindClusters(object=sc10x,reduction.type=red,dims.use=1:pc.use,resolution=res.use,print.output=0,save.SNN=TRUE,force.recalc=TRUE)
sc10x <- BuildClusterTree(sc10x,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
sc10x <- StashIdent(object=sc10x,save.name=paste0("res",res.use))
......@@ -1006,10 +1006,101 @@ scCCA <- function(sc10x.1,sc10x.2,nm.1="D17",nm.2="D27",cc=FALSE){
#MetageneBicorPlot(sc10x,dims.eval=1:50,grouping.var="patient")
#dev.off()
gc()
sc10x <- AlignSubspace(sc10x,reduction.type="cca",grouping.var="patient",dims.align=1:25)
gc()
postscript("./analysis/cca/CCA.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
DimPlot(sc10x,reduction.use="cca.aligned",group.by="patient",do.return=FALSE,pt.size=0.05)
dev.off()
postscript("./analysis/cca/Violin.CCA1.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot <- VlnPlot(sc10x,features.plot="ACC1",group.by="patient",size.title.use=20,point.size.use=0.05)
plot(plot)
dev.off()
postscript("./analysis/cca/Violin.CCA2.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot <- VlnPlot(sc10x,features.plot="ACC2",group.by="patient",size.title.use=20,point.size.use=0.05)
plot(plot)
dev.off()
results <- list(
sc10x=sc10x,
genes.hvg.cca=genes.hvg.Comb
)
return(results)
}
sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=FALSE){
gc()
if (cc==TRUE){
sc10x.1 <- ScaleData(object=sc10x.1,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=50)
} else {
sc10x.1 <- ScaleData(object=sc10x.1,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=50)
}
gc()
gc()
if (cc==TRUE){
sc10x.2 <- ScaleData(object=sc10x.2,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=50)
} else {
sc10x.2 <- ScaleData(object=sc10x.2,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=50)
}
gc()
gc()
if (cc==TRUE){
sc10x.3 <- ScaleData(object=sc10x.3,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=50)
} else {
sc10x.3 <- ScaleData(object=sc10x.3,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=50)
}
gc()
sc10x.1 <- FindVariableGenes(sc10x.1,do.plot=FALSE)
sc10x.2 <- FindVariableGenes(sc10x.2,do.plot=FALSE)
sc10x.3 <- FindVariableGenes(sc10x.3,do.plot=FALSE)
genes.hvg.1 <- head(rownames(sc10x.1@hvg.info),1000)
genes.hvg.2 <- head(rownames(sc10x.2@hvg.info),1000)
genes.hvg.3 <- head(rownames(sc10x.3@hvg.info),1000)
genes.hvg.Comb <- unique(c(genes.hvg.1,genes.hvg.2,genes.hvg.3))
genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.1@scale.data))
genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.2@scale.data))
genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.3@scale.data))
sc10x.1@meta.data$patient <- nm.1
sc10x.2@meta.data$patient <- nm.2
sc10x.3@meta.data$patient <- nm.3
sc10x.1 <- RenameCells(sc10x.1,nm.1)
sc10x.2 <- RenameCells(sc10x.2,nm.2)
sc10x.3 <- RenameCells(sc10x.3,nm.3)
sc10x.list <- list(sc10x.1,sc10x.2,sc10x.3)
sc10x <- RunMultiCCA(sc10x.list,genes.use=genes.hvg.Comb,num.cc=50)
rm(sc10x.1)
rm(sc10x.2)
rm(sc10x.3)
rm(sc10x.list)
postscript("./analysis/cca/CCA.eps",paper="special",width=10,height=5,horizontal=FALSE)
DimPlot(sc10x,reduction.use="cca",group.by="patient",do.return=FALSE,pt.size=0.05)
dev.off()
postscript("./analysis/cca/Violin.CCA1.Raw.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot <- VlnPlot(sc10x,features.plot="CC1",group.by="patient",size.title.use=20,point.size.use=0.05)
plot(plot)
dev.off()
postscript("./analysis/cca/Violin.CCA2.Raw.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot <- VlnPlot(sc10x,features.plot="CC2",group.by="patient",size.title.use=20,point.size.use=0.05)
plot(plot)
dev.off()
#postscript("./analysis/cca/Bicor.eps",paper="special",width=10,height=5,horizontal=FALSE)
#MetageneBicorPlot(sc10x,dims.eval=1:50,grouping.var="patient")
#dev.off()
gc()
sc10x <- AlignSubspace(sc10x,reduction.type="cca",grouping.var="patient",dims.align=1:30)
gc()
postscript("./analysis/cca/CCA.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
DimPlot(sc10x,reduction.use="cca.aligned",group.by="patient",do.return=FALSE,pt.size=0.05)
dev.off()
postscript("./analysis/cca/Violin.CCA1.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot <- VlnPlot(sc10x,features.plot="ACC1",group.by="patient",size.title.use=20,point.size.use=0.05)
plot(plot)
......@@ -1135,3 +1226,74 @@ scPseudotime <- function(sc10x,i,ds){
BEAM_res <- BEAM_res[,c("gene_short_name","pval","qval")]
#plot_genes_branched_heatmap(scMonocle[rownames(subset(BEAM_res,qval<1e-4)),],branch_point=1,num_clusters=5,cores=50,use_gene_short_name=TRUE,show_rownames=TRUE)
}
RenameCells <- function(object, add.cell.id = NULL, new.names = NULL,
for.merge = FALSE) {
if (missing(add.cell.id) && missing(new.names)) {
stop("One of 'add.cell.id' and 'new.names' must be set")
}
if (!missing(add.cell.id) && !missing(new.names)) {
stop("Only one of 'add.cell.id' and 'new.names' must be set")
}
if (!missing(add.cell.id)) {
new.cell.names <- paste(add.cell.id, object@cell.names, sep = "_")
new.rawdata.names <- paste(add.cell.id, colnames(object@raw.data), sep = "_")
} else {
if(ncol(object@raw.data) != ncol(object@data)) {
stop("raw.data contains a different number of cells than data")
}
if(any(colnames(object@raw.data) != colnames(object@data))){
stop("cells in raw.data are different than the cells in data")
}
if (length(new.names) == length(object@cell.names)) {
new.cell.names <- new.names
new.rawdata.names <- new.names
} else {
stop("the length of 'new.names' (", length(new.names), ") must be the ",
"same as the length of 'object@cell.names' (",
length(object@cell.names), ")")
}
}
colnames(object@raw.data) <- new.rawdata.names
rownames(object@meta.data) <- new.cell.names
object@cell.names <- new.cell.names
if (for.merge) {
return(object)
}
colnames(object@data) <- new.cell.names
if (!is.null(object@scale.data)) {
colnames(object@scale.data) <- new.cell.names
}
names(object@ident) <- new.cell.names
if (length(object@dr) > 0) {
for (dr in names(object@dr)) {
rownames(object@dr[[dr]]@cell.embeddings) <- new.cell.names
}
}
if (nrow(object@snn) == length(new.cell.names)) {
colnames(object@snn) <- new.cell.names
rownames(object@snn) <- new.cell.names
}
if (!is.null(object@kmeans)) {
if (!is.null(object@kmeans@gene.kmeans.obj)) {
colnames(object@kmeans@gene.kmeans.obj$centers) <- new.cell.names
}
if (!is.null(object@kmeans@cell.kmeans.obj)) {
names(object@kmeans@cell.kmeans.obj$cluster) <- new.cell.names
}
}
# rownames(object@spatial@mix.probs) <- new.cell.names
return(object)
}
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment