Add % var to PC plots

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)