diff --git a/r.scripts/sc-TissueMapper.R b/r.scripts/sc-TissueMapper.R index 196bb868929a39fa62e848ef433c1fd63dcb4cd5..71d98579bb582df4326eda8683d27572aa1f90ac 100644 --- a/r.scripts/sc-TissueMapper.R +++ b/r.scripts/sc-TissueMapper.R @@ -103,17 +103,14 @@ scCellCycle <- function(sc10x,sub=FALSE){ 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) - tryCatch({ - 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(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() - },error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")}) - + 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(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(paste0(folder,"PCA_cc.Raw.eps")) plot <- PCAPlot(object=sc10x,do.return=TRUE) @@ -349,6 +346,8 @@ scStress <- function(sc10x,stg="go",res.use=1,pc.use=10,cut=0.95){ sc10x.Stress <- t(sc10x.Stress) sc10x.Stress <- sc10x.Stress[,apply(sc10x.Stress,2,var)!=0] sc10x.Stress.pca <- prcomp(sc10x.Stress,center=TRUE,scale.=TRUE) + sc10x.Stress.pca.pc1var <- round((sc10x.Stress.pca$sdev[1]^2)/sum(sc10x.Stress.pca$sdev^2)*100,0) + sc10x.Stress.pca.pc2var <- round((sc10x.Stress.pca$sdev[2]^2)/sum(sc10x.Stress.pca$sdev^2)*100,0) sc10x.Stress.pca <- sc10x.Stress.pca$x[,1:2] colnames(x=sc10x.Stress.pca) <- paste0("Stress",1:2) if (skewness(sc10x.Stress.pca[,1])<0){ @@ -366,6 +365,8 @@ scStress <- function(sc10x,stg="go",res.use=1,pc.use=10,cut=0.95){ plot <- plot+geom_density2d(color="black",bins=25,alpha=0.5) 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))) + plot <- plot+xlab(paste0("Stress PC1 (",sc10x.Stress.pca.pc1var,"%)")) + plot <- plot+ylab(paste0("Stress PC2 (",sc10x.Stress.pca.pc2var,"%)")) plot(plot) dev.off() @@ -400,11 +401,13 @@ scStress <- function(sc10x,stg="go",res.use=1,pc.use=10,cut=0.95){ 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))) plot <- plot+geom_vline(xintercept=cut.x,color="red",lwd=2.5) + plot <- plot+xlab(paste0("Stress PC1 (",sc10x.Stress.pca.pc1var,"%)")) + plot <- plot+ylab(paste0("Stress PC2 (",sc10x.Stress.pca.pc2var,"%)")) plot(plot) dev.off() #Subsample all cells (+Stress) to better visualize their clustering - if (ncol(sc10x@data)>2500){ + if (ncol(sc10x@data)<2500){ rnd <- sample(1:ncol(sc10x@data),2500) } else { rnd <- 1:ncol(sc10x@data) @@ -848,6 +851,8 @@ scNE <- function(sc10x,neg="EurUro",cut=0.95){ sc10x.NE <- t(sc10x.NE) sc10x.NE <- sc10x.NE[,apply(sc10x.NE,2,var)!=0] sc10x.NE.pca <- prcomp(sc10x.NE,center=TRUE,scale.=TRUE) + sc10x.NE.pca.pc1var <- round((sc10x.NE.pca$sdev[1]^2)/sum(sc10x.NE.pca$sdev^2)*100,0) + sc10x.NE.pca.pc2var <- round((sc10x.NE.pca$sdev[2]^2)/sum(sc10x.NE.pca$sdev^2)*100,0) sc10x.NE.pca <- sc10x.NE.pca$x[,1:2] colnames(x=sc10x.NE.pca) <- paste0("NE",1:2) if (skewness(sc10x.NE.pca[,1])<0){ @@ -864,6 +869,8 @@ scNE <- function(sc10x,neg="EurUro",cut=0.95){ plot <- DimPlot(object=sc10x,reduction.use="NE",pt.size=2.5,do.return=TRUE,vector.friendly=FALSE) 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))) + plot <- plot+xlab(paste0("NE PC1 (",sc10x.NE.pca.pc1var,"%)")) + plot <- plot+ylab(paste0("NE PC2 (",sc10x.NE.pca.pc2var,"%)")) plot(plot) dev.off() @@ -897,11 +904,13 @@ scNE <- function(sc10x,neg="EurUro",cut=0.95){ 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))) plot <- plot+geom_vline(xintercept=cut.x,color="red",lwd=2.5) + plot <- plot+xlab(paste0("NE PC1 (",sc10x.NE.pca.pc1var,"%)")) + plot <- plot+ylab(paste0("NE PC2 (",sc10x.NE.pca.pc2var,"%)")) plot(plot) dev.off() #Subsample all cells (+NE) to better visualize their clustering - if (ncol(sc10x@data)>2500){ + if (ncol(sc10x@data)<2500){ rnd <- sample(1:ncol(sc10x@data),2500) } else { rnd <- 1:ncol(sc10x@data)