Add code for SingleGeneAnalysis

r.scripts/sc_SingleGeneAnalysis.R

+library(Seurat)
+library(pastecs)
+setwd("/work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/RUN/Pr/sc-TissueMapper_Pr")
+load("./analysis/sc10x.Pr.ALL.cluster.NOStress.IDepi+st+ne.Merge.Rda")
+
+for (i in c("ESR1","ESR2","GPER1")){
+sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL")
+sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=names(sc10x.Group@data[rownames(sc10x.Group@data)==i,][sc10x.Group@data[rownames(sc10x.Group@data)==i,]==0]),ident.use="DropOut")
+sc10x.Group.sub <- SubsetData(object=sc10x.Group,ident.use="ALL")
+postscript(paste0("./WR/Hist_",i,".eps"))
+hist(sc10x.Group.sub@data[rownames(sc10x.Group.sub@data)==i,],breaks=50,col="grey",xlab=paste0(i," UMI Counts"),main=paste0(i," Expression Minus Dropouts"),prob=TRUE,plot=TRUE)
+d <- density(sc10x.Group.sub@data[rownames(sc10x.Group.sub@data)==i,])
+lines(d,lwd=5)
+ts_y<-ts(d$y)
+tp<-turnpoints(ts_y)
+pit <- d$x[tp$pits]
+abline(v=pit,col="green",lwd=2.5)
+dev.off()
+
+hi <- names(sc10x.Group@data[rownames(sc10x.Group@data)==i,][sc10x.Group@data[rownames(sc10x.Group@data)==i,]>pit])
+sc10x.Group@ident <- plyr::mapvalues(sc10x.Group@ident,from="ALL",to="Low")
+sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=hi,ident.use="High")
+sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("DropOut","Low","High"))
+
+postscript(paste0("./WR/Plot_",i,".eps"))
+GenePlot(object=sc10x.Group,gene1="tSNE_1",gene2=i,col.use=c("black","red","green"))
+title(main=i,line=2.5)
+abline(h=pit,col="green",lwd=2.5)
+legend("bottomright",legend=c("DropOut","Low","High"),bg="white",pch=16,col=c("black","red","green"))
+dev.off()
+
+postscript(paste0("./WR/Plot_",i,".vs.COL1A1.eps"))
+GenePlot(object=sc10x.Group,gene1="COL1A1",gene2=i,col.use=c("black","red","green"))
+title(main=i,line=2.5)
+abline(h=pit,col="green",lwd=2.5)
+legend("bottomright",legend=c("DropOut","Low","High"),bg="white",pch=16,col=c("black","red","green"))
+dev.off()
+
+write.table(table(sc10x.Group@ident,sc10x.Group@meta.data$'SubClust.ID+NE'),file=paste0("./WR/Table_",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
+
+sc10x.Group <- StashIdent(object=sc10x.Group,save.name=paste0("gene_",i))
+
+if (i=="ESR1"){
+  ESR1.pit <- pit
+}else if (i=="ESR2"){
+  ESR2.pit <- pit
+}
+}
+
+sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL")
+postscript(paste0("./WR/Plot_ESR1.vs.ESR2.eps"))
+GenePlot(object=sc10x.Group,gene1="ESR1",gene2="ESR2",col.use=c("black","red","green"))
+title(main="ESR1 vs ESR2",line=2.5)
+abline(v=ESR1.pit,col="green",lwd=2.5)
+abline(h=ESR2.pit,col="green",lwd=2.5)
+dev.off()
+