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

Finalize? stress comparisson (van den Brink) code

parent 52315a1a
2 merge requests!6Develop,!5Stress compare
Hide whitespace changes
Inline Side-by-side
r.scripts/StressCompare.R 100644 → 100755
+ 129
− 124
View file @ a684c35f
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(VennDiagram)
library(Cairo)
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)]))
genes.stress.go <- read_csv("../genesets/DEG_C2.CGP.M10970.txt")
genes.stress.go <- genes.stress[2:nrow(genes.stress.go),]
genes.stress.go <- as.list(genes.stress.go)
genes.stress.dws <- read_delim("../genesets/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress.dws <- genes.stress.dws[1]
colnames(genes.stress.dws) <- "scDWS.Stress"
genes.stress.dws <- as.list(genes.stress.dws)
genes.stress.vdb <- read_delim("./vanderBrink.Stress.mus.tsv","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
colnames(genes.stress.vdb) <- "vdBrink.Stress"
genes.homolog <- read_delim("./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)),"./Venn.tiff")
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
gene.stress.pct <- apply(gene.exp,2,function(x) sum(x[rownames(gene.exp) %in% unlist(genes.stress)])/sum(x))
#histo <- hist(gene.stress.pct,breaks=100,prob=TRUE,plot=TRUE,main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes")
#abline(v=.0575,col="red")
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL+Stress")
cell.index.stress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="Stress"])
cell.index.stress <- which(colnames(gene.exp) %in% cell.index.stress)
cell.index.ntstress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="ALL"])
cell.index.ntstress <- which(colnames(gene.exp) %in% cell.index.ntstress)
histo.stress <- hist(gene.stress.pct[cell.index.stress],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")
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("./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="")
abline(v=0.0575,col="red")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks>0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks>0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks<=0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks<=0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks>0.035][!is.na(histo.stress$counts[histo.stress$breaks>0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks<=0.035][!is.na(histo.stress$counts[histo.stress$breaks<=0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,"Not Stressed",col="green")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,"Stressed",col="red")
dev.off()
}
rm(i)
rm(genes.stress)
rm(gene.stress.pct)
rm(cell.index.stress)
rm(cell.index.ntstress)
rm(histo.stress)
rm(histo.ntstress)
for (cut in c(0.9,0.95)){
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL")
sc10x.Stress.dta <- as.data.frame(as.matrix(sc10x.Stress@data[rownames(sc10x.Stress@scale.data) %in% unlist(genes.stress),]))
#Run PCA of subsetted stress genes
sc10x.Stress.dta <- t(sc10x.Stress.dta)
sc10x.Stress.dta <- sc10x.Stress.dta[,apply(sc10x.Stress.dta,2,var)!=0]
sc10x.Stress.dta.pca <- prcomp(sc10x.Stress.dta,center=TRUE,scale.=TRUE)
sc10x.Stress.dta.pca <- sc10x.Stress.dta.pca$x[,1:2]
colnames(x=sc10x.Stress.dta.pca) <- paste0("Stress",1:2)
if (skewness(sc10x.Stress.dta.pca[,1])<0){
sc10x.Stress.dta.pca[,1] <- (-sc10x.Stress.dta.pca[,1])
}
if (skewness(sc10x.Stress.dta.pca[,2])<0){
sc10x.Stress.dta.pca[,2] <- (-sc10x.Stress.dta.pca[,2])
}
sc10x.dta <- sc10x.Stress
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="cell.embeddings",new.data=sc10x.Stress.dta.pca)
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="key",new.data="Stress")
#CDS
cdf <- ecdf(GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
cut.x <- quantile(cdf,probs=cut)
assign(i,GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
assign(paste0("cut.",i,".",cut),cut.x)
}}
rm(cut)
rm(i)
rm(genes.stress)
rm(sc10x.Stress.dta)
rm(sc10x.Stress.dta.pca)
rm(sc10x.dta)
rm(cdf)
rm(cut.x)
comp <- cbind(vdb,go[match(rownames(vdb),rownames(go))],dws[match(rownames(vdb),rownames(dws))])
colnames(comp) <- c("vdb","go","dws")
comp <- data.frame(comp)
Cairo(paste0("./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.9)+geom_hline(yintercept=cut.go.0.9)+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.9,][comp[comp$vdb>cut.vdb.0.9,"go"]>cut.go.0.9,])/nrow(comp[comp$vdb>cut.vdb.0.9,])*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.9,][comp[comp$vdb<=cut.vdb.0.9,"go"]<=cut.go.0.9,])/nrow(comp[comp$vdb<=cut.vdb.0.9,])*100,1),"%"))
dev.off()
Cairo(paste0("./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.9)+geom_hline(yintercept=cut.dws.0.9)+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.9,][comp[comp$vdb>cut.vdb.0.9,"dws"]>cut.dws.0.9,])/nrow(comp[comp$vdb>cut.vdb.0.9,])*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.9,][comp[comp$vdb<=cut.vdb.0.9,"dws"]<=cut.dws.0.9,])/nrow(comp[comp$vdb<=cut.vdb.0.9,])*100,1),"%"))
dev.off()
Cairo(paste0("./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.9)+geom_hline(yintercept=cut.dws.0.9)+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.9,][comp[comp$go>cut.go.0.9,"dws"]>cut.dws.0.9,])/nrow(comp[comp$go>cut.go.0.9,])*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.9,][comp[comp$go<=cut.go.0.9,"dws"]<=cut.dws.0.9,])/nrow(comp[comp$go<=cut.go.0.9,])*100,1),"%"))
dev.off()
rm(comp)
\ No newline at end of file
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(VennDiagram)
library(Cairo)
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)]))
genes.stress.go <- read_csv("../genesets/DEG_C2.CGP.M10970.txt")
genes.stress.go <- genes.stress.go[2:nrow(genes.stress.go),]
genes.stress.go <- as.list(genes.stress.go)
genes.stress.dws <- read_delim("../genesets/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress.dws <- genes.stress.dws[1]
colnames(genes.stress.dws) <- "scDWS.Stress"
genes.stress.dws <- as.list(genes.stress.dws)
#genes.stress.vdb <- read_delim("./vanderBrink.Stress.mus.tsv","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
#colnames(genes.stress.vdb) <- "vdBrink.Stress"
genes.stress.vdb <- c("Actg1__chr11","Ankrd1__chr19","Arid5a__chr1","Atf3__chr1","Atf4__chr15","Bag3__chr7","Bhlhe40__chr6","Brd2__chr17","Btg1__chr10","Btg2__chr1","Ccnl1__chr3","Ccrn4l__chr3","Cebpb__chr2","Cebpd__chr16","Cebpg__chr7","Csrnp1__chr9","Cxcl1__chr5","Cyr61__chr3","Dcn__chr10","Ddx3x__chrX","Ddx5__chr11","Des__chr1","Dnaja1__chr4","Dnajb1__chr8","Dnajb4__chr3","Dusp1__chr17","Dusp8__chr7","Egr1__chr18","Egr2__chr10","Eif1__chr11","Eif5__chr12","Erf__chr7","Errfi1__chr4","Fam132b__chr1","Fos__chr12","Fosb__chr7","Fosl2__chr5","Gadd45a__chr6","Gcc1__chr6","Gem__chr4","H3f3b__chr11","Hipk3__chr2","Hsp90aa1__chr12","Hsp90ab1__chr17","Hspa1a__chr17","Hspa1b__chr17","Hspa5__chr2","Hspa8__chr9","Hspb1__chr5","Hsph1__chr5","Id3__chr4","Idi1__chr13","Ier2__chr8","Ier3__chr17","Ifrd1__chr12","Il6__chr5","Irf1__chr11","Irf8__chr8","Itpkc__chr7","Jun__chr4","Junb__chr8","Jund__chr8","Klf2__chr8","Klf4__chr4","Klf6__chr13","Klf9__chr19","Litaf__chr16","Lmna__chr3","Maff__chr15","Mafk__chr5","Mcl1__chr3","Midn__chr10","Mir22hg__chr11","Mt1__chr8","Mt2__chr8","Myadm__chr7","Myc__chr15","Myd88__chr9","Nckap5l__chr15","Ncoa7__chr10","Nfkbia__chr12","Nfkbiz__chr16","Nop58__chr1","Nppc__chr1","Nr4a1__chr15","Odc1__chr12","Osgin1__chr8","Oxnad1__chr14","Pcf11__chr7","Pde4b__chr4","Per1__chr11","Phlda1__chr10","Pnp__chr14","Pnrc1__chr4","Ppp1cc__chr5","Ppp1r15a__chr7","Pxdc1__chr13","Rap1b__chr10","Rassf1__chr9","Rhob__chr12","Rhoh__chr5","Ripk1__chr13","Sat1__chrX","Sbno2__chr10","Sdc4__chr2","Serpine1__chr5","Skil__chr3","Slc10a6__chr5","Slc38a2__chr15","Slc41a1__chr1","Socs3__chr11","Sqstm1__chr11","Srf__chr17","Srsf5__chr12","Srsf7__chr17","Stat3__chr11","Tagln2__chr1","Tiparp__chr3","Tnfaip3__chr10","Tnfaip6__chr2","Tpm3__chr3","Tppp3__chr8","Tra2a__chr6","Tra2b__chr16","Trib1__chr15","Tubb4b__chr2","Tubb6__chr18","Ubc__chr5","Usp2__chr9","Wac__chr18","Zc3h12a__chr4","Zfand5__chr19","Zfp36__chr7","Zfp36l1__chr12","Zfp36l2__chr17","Zyx__chr6","Gadd45g__chr13","Hspe1__chr1","Ier5__chr1","Kcne4__chr1")
genes.stress.vdb <- gsub("\\__.*","",genes.stress.vdb)
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")
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
gene.stress.pct <- apply(gene.exp,2,function(x) sum(x[rownames(gene.exp) %in% unlist(genes.stress)])/sum(x))
#histo <- hist(gene.stress.pct,breaks=100,prob=TRUE,plot=TRUE,main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes")
#abline(v=.0575,col="red")
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL+Stress")
cell.index.stress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="Stress"])
cell.index.stress <- which(colnames(gene.exp) %in% cell.index.stress)
cell.index.ntstress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="ALL"])
cell.index.ntstress <- which(colnames(gene.exp) %in% cell.index.ntstress)
histo.stress <- hist(gene.stress.pct[cell.index.stress],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")
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")
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="")
abline(v=0.0575,col="red")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks>0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks>0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks<=0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks<=0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks>0.035][!is.na(histo.stress$counts[histo.stress$breaks>0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks<=0.035][!is.na(histo.stress$counts[histo.stress$breaks<=0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,"Not Stressed",col="green")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,"Stressed",col="red")
dev.off()
}
rm(i)
rm(genes.stress)
rm(gene.stress.pct)
rm(cell.index.stress)
rm(cell.index.ntstress)
rm(histo.stress)
rm(histo.ntstress)
for (cut in c(0.9,0.95)){
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL")
sc10x.Stress.dta <- as.data.frame(as.matrix(sc10x.Stress@data[rownames(sc10x.Stress@scale.data) %in% unlist(genes.stress),]))
#Run PCA of subsetted stress genes
sc10x.Stress.dta <- t(sc10x.Stress.dta)
sc10x.Stress.dta <- sc10x.Stress.dta[,apply(sc10x.Stress.dta,2,var)!=0]
sc10x.Stress.dta.pca <- prcomp(sc10x.Stress.dta,center=TRUE,scale.=TRUE)
sc10x.Stress.dta.pca <- sc10x.Stress.dta.pca$x[,1:2]
colnames(x=sc10x.Stress.dta.pca) <- paste0("Stress",1:2)
if (skewness(sc10x.Stress.dta.pca[,1])<0){
sc10x.Stress.dta.pca[,1] <- (-sc10x.Stress.dta.pca[,1])
}
if (skewness(sc10x.Stress.dta.pca[,2])<0){
sc10x.Stress.dta.pca[,2] <- (-sc10x.Stress.dta.pca[,2])
}
sc10x.dta <- sc10x.Stress
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="cell.embeddings",new.data=sc10x.Stress.dta.pca)
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="key",new.data="Stress")
#CDS
cdf <- ecdf(GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
cut.x <- quantile(cdf,probs=cut)
assign(i,GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
assign(paste0("cut.",i,".",cut),cut.x)
}}
rm(cut)
rm(i)
rm(genes.stress)
rm(sc10x.Stress.dta)
rm(sc10x.Stress.dta.pca)
rm(sc10x.dta)
rm(cdf)
rm(cut.x)
comp <- cbind(vdb,go[match(rownames(vdb),rownames(go))],dws[match(rownames(vdb),rownames(dws))])
colnames(comp) <- c("vdb","go","dws")
comp <- data.frame(comp)
Cairo(paste0("./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")
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")
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()
rm(comp)
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