add sc_LineageSubClust.R test script

r.scripts/sc_LineageSubClust.R

+gc()
+.libPaths("/home2/ghenry/R/x86_64-pc-linux-gnu-library/3.3")
+library(optparse)
+library(Seurat)
+
+#retrive command line options
+option_list=list(
+  make_option("--p",action="store",default="Pr",type='character',help="Project Name"),
+  make_option("--g",action="store",default="ALL",type='character',help="Group To Analyze"),
+  make_option("--st",action="store",default=TRUE,type='logical',help="Analyze Data With Stress Pops Removed?"),
+  make_option("--pc",action="store",default=10,type='integer',help="Number Of Principal Components")
+)
+opt=parse_args(OptionParser(option_list=option_list))
+rm(option_list)
+Project.Name <- opt$p
+PC.Max <- opt$pc
+if (opt$st==TRUE){
+  sub.folder <- "/NoStress"
+  sub.file <- ".NOStress"
+} else {
+  sub.folder <- ""
+  sub.file <- ""
+  }
+
+#load data
+setwd("../")
+load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData"))
+load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda"))
+sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"))
+rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"))
+sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Lineage")
+sc10x.Group.Epi <- SubsetData(object=sc10x.Group,ident.use="Epi")
+sc10x.Group.St <- SubsetData(object=sc10x.Group,ident.use="St")
+sc10x.Group.All <- sc10x.Group
+
+sc10x.Group.Epi <- FindVariableGenes(object=sc10x.Group.Epi,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
+sc10x.Group.St <- FindVariableGenes(object=sc10x.Group.St,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
+
+sc10x.Group.Epi <- RunPCA(object=sc10x.Group.Epi,pc.genes=sc10x.Group.Epi@var.genes,do.print=FALSE,pcs.compute=50)
+sc10x.Group.St <- RunPCA(object=sc10x.Group.St,pc.genes=sc10x.Group.St@var.genes,do.print=FALSE,pcs.compute=50)
+
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.Epi.eps"))
+plot <- PCElbowPlot(object=sc10x.Group.Epi,num.pc=50)
+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))
+plot(plot)
+dev.off()
+
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.St.eps"))
+plot <- PCElbowPlot(object=sc10x.Group.St,num.pc=50)
+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))
+plot(plot)
+dev.off()
+
+sc10x.Group.Epi <- RunTSNE(object=sc10x.Group.Epi,dims.use=1:PC.Max,do.fast=TRUE)
+postscript(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Sample.Epi.eps"))
+plot <- TSNEPlot(object=sc10x.Group.Epi,group.by="samples",pt.size=2.5,do.return=TRUE)
+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)
+dev.off()
+
+sc10x.Group.St <- RunTSNE(object=sc10x.Group.St,dims.use=1:PC.Max,do.fast=TRUE)
+postscript(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Sample.St.eps"))
+plot <- TSNEPlot(object=sc10x.Group.St,group.by="samples",pt.size=2.5,do.return=TRUE)
+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)
+dev.off()
+
+for (j in (c("Epi","St"))){
+  if (j=="Epi"){
+    sc10x.Group <- sc10x.Group.Epi
+  } else {
+    sc10x.Group <- sc10x.Group.St
+  }
+  for (i  in c(5,2,1,0.5,0.2,0.1)){
+    gc()  
+    sc10x.Group <- FindClusters(object=sc10x.Group,reduction.type="pca",dims.use=1:PC.Max,resolution=i,print.output=0,save.SNN=TRUE,force.recalc=TRUE)
+    sc10x.Group <- BuildClusterTree(sc10x.Group,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
+    sc10x.Group <- StashIdent(object=sc10x.Group,save.name=paste0("res",i))
+    postscript(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Global.",j,".res",i,".eps"))
+    plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE)
+    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)
+    dev.off()
+}}
\ No newline at end of file