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

Create Pd analysis workflow and fix some of the outputs

parent b7016e66
Branches
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
r.scripts/sc-TissueMapper.R
+ 55
− 32
View file @ 80d84d6b
......@@ -75,7 +75,7 @@ scSubset <- function(sc10x,i="ALL",g="ALL"){
}
scCellCycle <- function(sc10x){
scCellCycle <- function(sc10x,sub=FALSE){
#Runs Seurat based PCA analysis for cell cycle ID
#Inputs:
......@@ -86,6 +86,15 @@ scCellCycle <- function(sc10x){
#results[2] = s genes
#results[3] = g2m genes
#Make sub-folders if necessary
if (sub==FALSE){
folder <- "./analysis/qc/cc/"
} else {
folder <- paste0("./analysis/qc/cc/",sub,"/")
if (!dir.exists(folder)){
dir.create(folder)
}}
#score cell cycle and generate figures
genes.cc <- readLines(con="./genesets/regev_lab_cell_cycle_genes.txt")
genes.s <- genes.cc[1:43]
......@@ -94,16 +103,16 @@ scCellCycle <- function(sc10x){
sc10x <- ScaleData(object=sc10x,display.progress=FALSE,do.par=TRUE,num.cores=45)
sc10x <- CellCycleScoring(object=sc10x,s.genes=genes.s,g2m.genes=genes.g2m,set.ident=TRUE)
postscript("./analysis/qc/cc/Ridge_cc.Raw.eps")
postscript(paste0(folder,"Ridge_cc.Raw.eps"))
plot <- RidgePlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),y.log=TRUE,nCol=2,do.return=TRUE)
plot(plot)
dev.off()
postscript("./analysis/qc/cc/Violin_cc.Raw.eps")
postscript(paste0(folder,"Violin_cc.Raw.eps"))
plot <- VlnPlot(object=sc10x,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),nCol=2,point.size.use=1,size.title.use=20,x.lab.rot=TRUE)
plot(plot)
dev.off()
sc10x <- RunPCA(object=sc10x,pc.genes=c(genes.s,genes.g2m),do.print=FALSE,pcs.store=2)
postscript("./analysis/qc/cc/PCA_cc.Raw.eps")
postscript(paste0(folder,"PCA_cc.Raw.eps"))
plot <- PCAPlot(object=sc10x,do.return=TRUE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
......@@ -113,7 +122,7 @@ scCellCycle <- function(sc10x){
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
gc()
sc10x <- RunPCA(object=sc10x,pc.genes=c(genes.s,genes.g2m),do.print=FALSE,pcs.store=2)
postscript("./analysis/qc/cc/PCA_cc.Norm.eps")
postscript(paste0(folder,"PCA_cc.Norm.eps"))
plot <- PCAPlot(object=sc10x,do.return=TRUE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
......@@ -129,7 +138,7 @@ scCellCycle <- function(sc10x){
}
scQC <- function(sc10x,lg=500,hg=2500,lm=-Inf,hm=0.1){
scQC <- function(sc10x,lg=500,hg=2500,lm=-Inf,hm=0.1,sub=FALSE){
#QC and filter Seurat object
#Inputs:
......@@ -146,18 +155,27 @@ scQC <- function(sc10x,lg=500,hg=2500,lm=-Inf,hm=0.1){
#result[4] = filtered cell count
#result[5] = filtered gene count
#Make sub-folders if necessary
if (sub==FALSE){
folder <- "./analysis/qc/"
} else {
folder <- paste0("./analysis/qc/",sub,"/")
if (!dir.exists(folder)){
dir.create(folder)
}}
#Calculate %mito genes
mito.genes <- grep(pattern="^MT-",x=rownames(x=sc10x@data),value=TRUE)
percent.mito <- colSums(sc10x@raw.data[mito.genes,])/colSums(sc10x@raw.data)
percent.mito <- colSums(as.matrix(sc10x@raw.data[mito.genes,]))/colSums(as.matrix(sc10x@raw.data))
sc10x <- AddMetaData(object=sc10x,metadata=percent.mito,col.name="percent.mito")
#Plot raw stats
sc10x <- SetAllIdent(object=sc10x,id="ALL")
postscript("./analysis/qc/Violin_qc.raw.eps")
postscript(paste0(folder,"Violin_qc.raw.eps"))
plot <- VlnPlot(object=sc10x,features.plot=c("nGene","nUMI","percent.mito"),nCol=3,do.return=TRUE)
plot(plot)
dev.off()
postscript("./analysis/qc/Plot_qc.raw.eps")
postscript(paste0(folder,"Plot_qc.raw.eps"))
par(mfrow=c(1,2))
GenePlot(object=sc10x,gene1="nUMI",gene2="percent.mito")
GenePlot(object=sc10x,gene1="nUMI",gene2="nGene")
......@@ -170,11 +188,11 @@ scQC <- function(sc10x,lg=500,hg=2500,lm=-Inf,hm=0.1){
sc10x <- NormalizeData(object=sc10x,normalization.method="LogNormalize",scale.factor=10000)
#Plot filtered stats
postscript("./analysis/qc/Violin_qc.filtered.eps")
postscript(paste0(folder,"Violin_qc.filtered.eps"))
plot <- VlnPlot(object=sc10x,features.plot=c("nGene","nUMI","percent.mito"),nCol=3,do.return=TRUE)
plot(plot)
dev.off()
postscript("./analysis/qc/Plot_qc.filtered.eps")
postscript(paste0(folder,"Plot_qc.filtered.eps"))
par(mfrow=c(1,2))
GenePlot(object=sc10x,gene1="nUMI",gene2="percent.mito")
GenePlot(object=sc10x,gene1="nUMI",gene2="nGene")
......@@ -549,12 +567,13 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){
for (j in 1:number.clusters){
qs <- get(paste0("results.",j))
if (j==1){
plotDensityCurves(qs,path.index=i,col=j,main=names(gs)[i],xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=2,cex.main=2.5)
plotDensityCurves(qs,path.index=i,col=viridis(number.clusters)[j],main=names(gs)[i],xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
} else {
plotDensityCurves(qs,path.index=i,add=TRUE,col=j,lwd=2)
plotDensityCurves(qs,path.index=i,add=TRUE,col=viridis(number.clusters)[j],lwd=5)
}}
leg <- paste0("Cluster ",1:number.clusters)
legend("topright",legend=leg,lty=1,col=1:number.clusters,lwd=2,cex=2,ncol=2,bty="n",pt.cex=1)
legend("topright",legend=leg,lty=1,col=viridis(number.clusters),lwd=5,cex=2,ncol=1,bty="n",pt.cex=2)
box(lwd=5)
dev.off()
}
......@@ -564,11 +583,12 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){
postscript(paste0("./analysis/cor/QuSAGE_",nm,"_Clust",i,".eps"))
for (j in 1:length(gs)){
if (j==1){
plotDensityCurves(qs,path.index=j,col=j,main=paste0("Cluster ",i),xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=2,cex.main=2.5)
plotDensityCurves(qs,path.index=j,col=j,main=paste0("Cluster ",i),xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=j,lwd=2)
plotDensityCurves(qs,path.index=j,add=TRUE,col=j,lwd=5)
}}
legend("topright",legend=names(gs),lty=1,col=1:length(gs),lwd=2,cex=2,ncol=2,bty="n",pt.cex=1)
legend("topright",legend=names(gs),lty=1,col=1:length(gs),lwd=5,cex=2,ncol=1,bty="n",pt.cex=2)
box(lwd=5)
dev.off()
}
......@@ -676,12 +696,13 @@ scQuSAGEsm <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){
for (j in clusters){
qs <- get(paste0("results.",j))
if (j==clusters[1]){
plotDensityCurves(qs,path.index=i,col=which(clusters == j),main=names(gs)[i],xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=2,cex.main=2.5)
plotDensityCurves(qs,path.index=i,col=which(clusters == j),main=names(gs)[i],xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
} else {
plotDensityCurves(qs,path.index=i,add=TRUE,col=which(clusters == j),lwd=2)
plotDensityCurves(qs,path.index=i,add=TRUE,col=which(clusters == j),lwd=5)
}}
leg <- clusters
legend("topright",legend=leg,lty=1,col=1:length(clusters),lwd=2,cex=2,ncol=2,bty="n",pt.cex=1)
legend("topright",legend=leg,lty=1,col=1:length(clusters),lwd=5,cex=2,ncol=1,bty="n",pt.cex=2)
dev.off()
}
......@@ -691,11 +712,12 @@ scQuSAGEsm <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){
postscript(paste0("./analysis/cor/QuSAGE_",nm,".",folder,"_",i,".eps"))
for (j in 1:length(gs)){
if (j==1){
plotDensityCurves(qs,path.index=j,col=j,main=i,xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=2,cex.main=2.5)
plotDensityCurves(qs,path.index=j,col=j,main=i,xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=j,lwd=2)
plotDensityCurves(qs,path.index=j,add=TRUE,col=j,lwd=5)
}}
legend("topright",legend=names(gs),lty=1,col=1:length(gs),lwd=2,cex=2,ncol=2,bty="n",pt.cex=1)
legend("topright",legend=names(gs),lty=1,col=1:length(gs),lwd=5,cex=2,ncol=1,bty="n",pt.cex=2)
dev.off()
}
......@@ -1049,7 +1071,7 @@ scCCA <- function(sc10x.1,sc10x.2,nm.1="D17",nm.2="D27",cc=FALSE){
}
sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=FALSE){
sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=FALSE,cca=25,acca=15){
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=45)
......@@ -1077,9 +1099,9 @@ sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=
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.1 <- head(rownames(sc10x.1@hvg.info),2000)
genes.hvg.2 <- head(rownames(sc10x.2@hvg.info),2000)
genes.hvg.3 <- head(rownames(sc10x.3@hvg.info),2000)
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))
......@@ -1092,7 +1114,7 @@ sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=
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)
sc10x <- RunMultiCCA(sc10x.list,genes.use=genes.hvg.Comb,num.cc=cca)
rm(sc10x.1)
rm(sc10x.2)
rm(sc10x.3)
......@@ -1109,12 +1131,12 @@ sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=
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()
postscript("./analysis/cca/Bicor.eps",paper="special",width=10,height=5,horizontal=FALSE)
MetageneBicorPlot(sc10x,dims.eval=1:cca,grouping.var="patient")
dev.off()
gc()
sc10x <- AlignSubspace(sc10x,reduction.type="cca",grouping.var="patient",dims.align=1:25)
sc10x <- AlignSubspace(sc10x,reduction.type="cca",grouping.var="patient",dims.align=1:acca)
gc()
postscript("./analysis/cca/CCA.Aligned.eps",paper="special",width=10,height=5,horizontal=FALSE)
......@@ -1316,3 +1338,4 @@ RenameCells <- function(object, add.cell.id = NULL, new.names = NULL,
return(object)
}
r.scripts/sc-TissueMapper_RUN.Pd.R 0 → 100644
+ 441
− 0
View file @ 80d84d6b
This diff is collapsed.
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