diff --git a/bash.scripts/sc_TissueMapper-Pd.sh b/bash.scripts/sc_TissueMapper-Pd.sh index 4c4ef806410f3e4bcf5c2f3ea1f5e3138cfcb31c..ba634e4d1d16a18814001c81718de640c0db6f56 100644 --- a/bash.scripts/sc_TissueMapper-Pd.sh +++ b/bash.scripts/sc_TissueMapper-Pd.sh @@ -13,3 +13,4 @@ module load R/3.4.1-gccmkl_20181025 sh ./sc_LinkData.sh Pd Rscript ../r.scripts/sc-TissueMapper_RUN.Pd.R +Rscript ../r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R diff --git a/r.scripts/StressCompare.R b/r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R similarity index 93% rename from r.scripts/StressCompare.R rename to r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R index 9353312c6b14f97fbd2a21418f8ce970333fc8f2..fc23ed8ce0e316175feeb4840d48c6b99a85b61e 100755 --- a/r.scripts/StressCompare.R +++ b/r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R @@ -9,6 +9,12 @@ library(RColorBrewer) library(VennDiagram) library(Cairo) +options(bitmapType="cairo") + +if (!dir.exists("../analysis/stress")){ + dir.create("../analysis/stress") +} + load("../analysis/sc10x.Stress.Rda") sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL") gene.exp <- as.data.frame(as.matrix(sc10x.Stress@raw.data[,colnames(sc10x.Stress@data)])) @@ -31,7 +37,7 @@ genes.stress.vdb <- data.frame(vdBrink.Stress=genes.stress.vdb) genes.homolog <- read_delim("../genesets/Ensemble.mus-hum.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) genes.stress.vdb <- list("vdBrink.Stress"=merge(genes.stress.vdb,genes.homolog[,c(2,4)],by.x="vdBrink.Stress",by.y="Gene name")[,2]) -venn.diagram(list(go=unlist(genes.stress.go),dws=unlist(genes.stress.dws),vdb=unlist(genes.stress.vdb)),"../stress/Venn.tiff") +venn.diagram(list(go=unlist(genes.stress.go),dws=unlist(genes.stress.dws),vdb=unlist(genes.stress.vdb)),"../analysis/stress/Venn.tiff") for (i in c("go","dws","vdb")){ @@ -53,7 +59,7 @@ for (i in c("go","dws","vdb")){ histo.ntstress <- hist(gene.stress.pct[cell.index.ntstress],breaks=seq(0,ceiling(max(gene.stress.pct)*100)/100,0.005),main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes",plot=FALSE) #abline(v=.0575,col="red") - Cairo(paste0("../stress/Hist.Stress.",i,".png"),width=1000,height=500,bg="white") + Cairo(paste0("../analysis/stress/Hist.Stress.",i,".png"),width=1000,height=500,bg="white") plot(histo.ntstress$mids,histo.ntstress$density,col=rgb(0,1,0,1/5),xlim=c(0,ceiling(max(gene.stress.pct)*100)/100),ylim=c(0,ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10),type="h",lend="square",lwd=10,xlab="Fraction of transcriptome associated to stress genes",ylab="Density") par(new=TRUE) plot(histo.stress$mids,histo.stress$density,col=rgb(1,0,0,1/5),xlim=c(0,ceiling(max(gene.stress.pct)*100)/100),ylim=c(0,ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10),type="h",lend="square",lwd=10,axes=FALSE,xlab="",ylab="") @@ -116,13 +122,13 @@ comp <- cbind(vdb,go[match(rownames(vdb),rownames(go))],dws[match(rownames(vdb), colnames(comp) <- c("vdb","go","dws") comp <- data.frame(comp) -Cairo(paste0("./Stress.vdb.vs.go.png"),width=1000,height=500,bg="white") +Cairo(paste0("../analysis/stress/Stress.vdb.vs.go.png"),width=1000,height=500,bg="white") ggplot(comp,aes(x=vdb,y=go))+geom_point(size=0.01)+geom_vline(xintercept=cut.vdb.0.95)+geom_hline(yintercept=cut.go.0.95)+annotate("text",x=max(comp$vdb)/4*3,y=max(comp$go)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb>cut.vdb.0.95,][comp[comp$vdb>cut.vdb.0.95,"go"]>cut.go.0.95,])/nrow(comp[comp$vdb>cut.vdb.0.95,])*100,1),"%"))+annotate("text",x=max(comp$vdb)/100,y=max(comp$go)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb<=cut.vdb.0.95,][comp[comp$vdb<=cut.vdb.0.95,"go"]<=cut.go.0.95,])/nrow(comp[comp$vdb<=cut.vdb.0.95,])*100,1),"%")) dev.off() -Cairo(paste0("./Stress.vdb.vs.dws.png"),width=1000,height=500,bg="white") +Cairo(paste0("../analysis/stress/Stress.vdb.vs.dws.png"),width=1000,height=500,bg="white") ggplot(comp,aes(x=vdb,y=dws))+geom_point(size=0.01)+geom_vline(xintercept=cut.vdb.0.95)+geom_hline(yintercept=cut.dws.0.95)+annotate("text",x=max(comp$vdb)/4*3,y=max(comp$dws)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb>cut.vdb.0.95,][comp[comp$vdb>cut.vdb.0.95,"dws"]>cut.dws.0.95,])/nrow(comp[comp$vdb>cut.vdb.0.95,])*100,1),"%"))+annotate("text",x=max(comp$vdb)/100,y=max(comp$dws)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb<=cut.vdb.0.95,][comp[comp$vdb<=cut.vdb.0.95,"dws"]<=cut.dws.0.95,])/nrow(comp[comp$vdb<=cut.vdb.0.95,])*100,1),"%")) dev.off() -Cairo(paste0("./Stress.go.vs.dws.png"),width=1000,height=500,bg="white") +Cairo(paste0("../analysis/stress/Stress.go.vs.dws.png"),width=1000,height=500,bg="white") ggplot(comp,aes(x=go,y=dws))+geom_point(size=0.01)+geom_vline(xintercept=cut.go.0.95)+geom_hline(yintercept=cut.dws.0.95)+annotate("text",x=max(comp$go)/4*3,y=max(comp$dws)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$go>cut.go.0.95,][comp[comp$go>cut.go.0.95,"dws"]>cut.dws.0.95,])/nrow(comp[comp$go>cut.go.0.95,])*100,1),"%"))+annotate("text",x=max(comp$go)/100,y=max(comp$dws)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$go<=cut.go.0.95,][comp[comp$go<=cut.go.0.95,"dws"]<=cut.dws.0.95,])/nrow(comp[comp$go<=cut.go.0.95,])*100,1),"%")) dev.off()