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

Update D27 run code with current settings

parent 67aaa246
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
r.scripts/sc-TissueMapper_RUN.D27.R
+ 14
− 84
View file @ 2f30f839
...@@ -77,29 +77,29 @@ if (!dir.exists("./analysis/deg")){ ...@@ -77,29 +77,29 @@ if (!dir.exists("./analysis/deg")){
#Retrieve command-line options #Retrieve command-line options
option_list=list( option_list=list(
make_option("--p",action="store",default="D27PrF",type='character',help="Project Name"), make_option("--p",action="store",default="D17PrF",type='character',help="Project Name"),
make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"), make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"),
make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"), make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"),
make_option("--hg",action="store",default=2500,type='integer',help="Threshold for cells with maximum genes"), make_option("--hg",action="store",default=2500,type='integer',help="Threshold for cells with maximum genes"),
make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"), make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"),
make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"), make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"),
make_option("--lx",action="store",default=0.15,type='numeric',help="x low threshold for hvg selection"), make_option("--lx",action="store",default=0.0125,type='numeric',help="x low threshold for hvg selection"),
make_option("--hx",action="store",default=35,type='numeric',help="x high threshold for hvg selection"), make_option("--hx",action="store",default=3,type='numeric',help="x high threshold for hvg selection"),
make_option("--ly",action="store",default=0.5,type='numeric',help="y low threshold for hvg selection"), make_option("--ly",action="store",default=0.5,type='numeric',help="y low threshold for hvg selection"),
make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"), make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"),
make_option("--pc",action="store",default=25,type='integer',help="Number of PCs to cacluate"), make_option("--pc",action="store",default=25,type='integer',help="Number of PCs to cacluate"),
make_option("--hpc",action="store",default=0.95,type='numeric',help="Max variance cutoff for PCs to use, pre-stress"), make_option("--hpc",action="store",default=0.7,type='numeric',help="Max variance cutoff for PCs to use, pre-stress"),
make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"), make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"),
make_option("--st",action="store",default="TRUE",type='logical',help="Remove stressed cells?"), make_option("--st",action="store",default="TRUE",type='logical',help="Remove stressed cells?"),
make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"), make_option("--stg",action="store",default="go",type='character',help="Geneset to use for stress ID"),
make_option("--hpc.poststress",action="store",default=0.95,type='numeric',help="Max variance cutoff for PCs to use, post-stress"), make_option("--hpc.poststress",action="store",default=0.7,type='numeric',help="Max variance cutoff for PCs to use, post-stress"),
make_option("--res.poststress",action="store",default=0.5,type='numeric',help="Resolution to cluster, post-stress"), make_option("--res.poststress",action="store",default=0.5,type='numeric',help="Resolution to cluster, post-stress"),
make_option("--ds",action="store",default=0,type='integer',help="Number of cells to downsample"), make_option("--ds",action="store",default=0,type='integer',help="Number of cells to downsample"),
make_option("--hpc.epi",action="store",default=0.75,type='numeric',help="Max variance cutoff for PCs to use, Epi"), make_option("--hpc.epi",action="store",default=0.7,type='numeric',help="Max variance cutoff for PCs to use, Epi"),
make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"), make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"),
make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"), make_option("--hpc.st",action="store",default=0.7,type='numeric',help="Max variance cutoff for PCs to use, St"),
make_option("--res.st",action="store",default=0.2,type='numeric',help="Resolution to cluster, St") make_option("--res.st",action="store",default=0.2,type='numeric',help="Resolution to cluster, St")
) )
opt=parse_args(OptionParser(option_list=option_list)) opt=parse_args(OptionParser(option_list=option_list))
rm(option_list) rm(option_list)
if (opt$lm==0){opt$lm=-Inf} if (opt$lm==0){opt$lm=-Inf}
...@@ -182,6 +182,7 @@ genes.hvg.epi <- results[[2]] ...@@ -182,6 +182,7 @@ genes.hvg.epi <- results[[2]]
pc.use.epi <- results[[3]] pc.use.epi <- results[[3]]
rm(results) rm(results)
sc10x.Epi <- scCluster(sc10x.Epi,pc.use=pc.use.epi,res.use=0.1,folder="epi")
sc10x.Epi <- scCluster(sc10x.Epi,pc.use=pc.use.epi,res.use=opt$res.epi,folder="epi") sc10x.Epi <- scCluster(sc10x.Epi,pc.use=pc.use.epi,res.use=opt$res.epi,folder="epi")
gene.set1 <- read_delim("./genesets/genes.deg.BE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set1 <- read_delim("./genesets/genes.deg.BE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
...@@ -194,12 +195,12 @@ gene.set1 <- gene.set1[1] ...@@ -194,12 +195,12 @@ gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1) gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "LE" names(gene.set1) <- "LE"
gene.set <- c(gene.set,gene.set1) gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE_SCGB.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set1 <- read_delim("./genesets/genes.deg.OE1.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1] gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1) gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "OE_SCGB" names(gene.set1) <- "OE_SCGB"
gene.set <- c(gene.set,gene.set1) gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE_KRT13.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set1 <- read_delim("./genesets/genes.deg.OE2.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1] gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1) gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "OE_KRT13" names(gene.set1) <- "OE_KRT13"
...@@ -218,6 +219,7 @@ genes.hvg.st <- results[[2]] ...@@ -218,6 +219,7 @@ genes.hvg.st <- results[[2]]
pc.use.st <- results[[3]] pc.use.st <- results[[3]]
rm(results) rm(results)
sc10x.St <- scCluster(sc10x.St,pc.use=pc.use.st,res.use=0.1,folder="st")
sc10x.St <- scCluster(sc10x.St,pc.use=pc.use.st,res.use=opt$res.st,folder="st") sc10x.St <- scCluster(sc10x.St,pc.use=pc.use.st,res.use=opt$res.st,folder="st")
gene.set1 <- read_delim("./genesets/genes.deg.Fib.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set1 <- read_delim("./genesets/genes.deg.Fib.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
...@@ -271,78 +273,6 @@ scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc") ...@@ -271,78 +273,6 @@ scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc")
scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc_NE") scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc_NE")
scTables(sc10x,i.1="Merge_Epi.dws.sc_St.dws.sc_NE",i.2="Merge_Epi.dws.sc_St.dws.sc") scTables(sc10x,i.1="Merge_Epi.dws.sc_St.dws.sc_NE",i.2="Merge_Epi.dws.sc_St.dws.sc")
genes.deg.Stress <- scDEG(sc10x.Stress,i="Stress",g.1="Stress",g.2="ALL",pct=0.5,t=5)
genes.deg.Epi <- scDEG(sc10x,i="Lin",g.1="Epi",g.2="St",t=2)
genes.deg.St <- scDEG(sc10x,i="Lin",g.1="St",g.2="Epi",t=2)
genes.deg.BE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sc",g.1="BE",g.2=c("LE","OE_SCGB","OE_KRT13"),pct=0.25,t=2)
genes.deg.LE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sc",g.1="LE",g.2=c("BE","LE","OE_SCGB"),pct=0.25,t=2)
genes.deg.OE_SCGB <- scDEG(sc10x.Epi.NE,i="Epi.dws.sc",g.1="OE_SCGB",g.2=c("BE","LE","OE_KRT13"),pct=0.25,t=2)
genes.deg.OE_KRT13 <- scDEG(sc10x.Epi.NE,i="Epi.dws.sc",g.1="OE_KRT13",g.2=c("BE","LE","OE_SCGB"),pct=0.25,t=2)
genes.deg.NE <- scDEG(sc10x.Epi.NE,i="NE",g.1="NE",g.2="ALL",pct=0.25,t=2)
genes.deg.Fib <- scDEG(sc10x.St,i="St.dws.sc",g.1="Fib",g.2=c("SM","Endo","Leu"),pct=0.25,t=2)
genes.deg.SM <- scDEG(sc10x.St,i="St.dws.sc",g.1="SM",g.2=c("Fib","Endo","Leu"),pct=0.25,t=2)
genes.deg.Endo <- scDEG(sc10x.St,i="St.dws.sc",g.1="Endo",g.2=c("Fib","SM","Leu"),pct=0.25,t=2)
genes.deg.Leu <- scDEG(sc10x.St,i="St.dws.sc",g.1="Leu",g.2=c("Fib","SM","Endo"),pct=0.25,t=2)
genes.deg.BE.unique <- setdiff(rownames(genes.deg.BE),Reduce(union,list(rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.LE.unique <- setdiff(rownames(genes.deg.LE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.OE_SCGB.unique <- setdiff(rownames(genes.deg.OE_SCGB),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.OE_KRT13.unique <- setdiff(rownames(genes.deg.OE_KRT13),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.NE.unique <- setdiff(rownames(genes.deg.NE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.Fib.unique <- setdiff(rownames(genes.deg.Fib),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.SM.unique <- setdiff(rownames(genes.deg.SM),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.Endo.unique <- setdiff(rownames(genes.deg.Endo),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Leu))))
genes.deg.Leu.unique <- setdiff(rownames(genes.deg.Leu),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE_SCGB),rownames(genes.deg.OE_KRT13),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo))))
genes.deg.5 <- c(genes.deg.BE.unique[1:5],genes.deg.LE.unique[1:5],genes.deg.OE_SCGB.unique[1:5],genes.deg.OE_KRT13.unique[1:5],genes.deg.NE.unique[1:5],genes.deg.Fib.unique[1:5],genes.deg.SM.unique[1:5],genes.deg.Endo.unique[1:5],genes.deg.Leu.unique[1:5])
genes.deg.5 <- rev(genes.deg.5)
genes.deg.10 <- c(genes.deg.BE.unique[1:10],genes.deg.LE.unique[1:10],genes.deg.OE_SCGB.unique[1:10],genes.deg.OE_KRT13.unique[1:10],genes.deg.NE.unique[1:10],genes.deg.Fib.unique[1:10],genes.deg.SM.unique[1:10],genes.deg.Endo.unique[1:10],genes.deg.Leu.unique[1:10])
genes.deg.10 <- rev(genes.deg.10)
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go_NE")
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go_NE")
sc10x@ident <- factor(sc10x@ident,levels=c("BE","LE","OE_SCGB","OE_KRT13","NE","Fib","SM","Endo","Leu"))
postscript("./analysis/deg/Dot.DEG.eps",paper="special",width=10,height=5,horizontal=FALSE)
DotPlot(sc10x,genes.deg.5,x.lab.rot=TRUE,plot.legend=TRUE,dot.scale=4)
dev.off()
postscript("./analysis/deg/Heatmap.DEG.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot(DoHeatmap(sc10x,genes.use=genes.deg.10,slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
dev.off()
for (i in ls(pattern="^genes.deg*unique")){
postscript(paste0("./analysis/deg/Violin.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
plot(VlnPlot(sc10x,features.plot=get(i)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
dev.off()
postscript(paste0("./analysis/deg/Ridge.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
plot(RidgePlot(sc10x,features.plot=get(i)[1:10],nCol=5,size.x.use=10,size.title.use=15))
dev.off()
postscript(paste0("./analysis/deg/Heatmap.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
plot(DoHeatmap(sc10x,genes.use=get(i),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
dev.off()
}
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="Stress")
postscript("./analysis/deg/Violin.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot(VlnPlot(sc10x.Stress,features.plot=rownames(genes.deg.Stress)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
dev.off()
postscript("./analysis/deg/Ridge.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot(RidgePlot(sc10x.Stress,features.plot=rownames(genes.deg.Stress)[1:10],nCol=5,size.x.use=10,size.title.use=15))
dev.off()
postscript("./analysis/deg/Heatmap.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
plot(DoHeatmap(sc10x.Stress,genes.use=rownames(genes.deg.Stress),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
dev.off()
for (i in ls(pattern="^genes.deg")){
write.table(get(i),file=paste0("./analysis/deg/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
}
save(list=ls(pattern="^genes.deg"),file="./analysis/DEG.Rda")
rm(list=ls(pattern="^genes.deg"))
save(list=ls(pattern="^sc10x"),file="./analysis/sc10x.Rda") save(list=ls(pattern="^sc10x"),file="./analysis/sc10x.Rda")
rm(list=ls(pattern="^sc10x")) rm(list=ls(pattern="^sc10x"))
save.image(file="./analysis/Data.RData") save.image(file="./analysis/Data.RData")
\ No newline at end of file
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