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

Remove out of date subcluster QuSAGE

parent 2096510b
No related merge requests found
Hide whitespace changes
Inline Side-by-side
r.scripts/sc_QuSAGE.LineageSubClust.R deleted 100644 → 0
+ 0
− 244
View file @ 2096510b
gc()
.libPaths("/home2/ghenry/R/x86_64-pc-linux-gnu-library/3.3")
library(optparse)
library(Seurat)
library(readr)
library(qusage)
#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("--r",action="store",default=0.1,type='numeric',help="Resolution"),
make_option("--st",action="store",default=TRUE,type='logical',help="Analyze Data With Stress Pops Removed?"),
make_option("--s",action="store",default=5000,type='integer',help="Number Of Cells To Downsample To")
)
opt=parse_args(OptionParser(option_list=option_list))
rm(option_list)
Project.Name <- opt$p
if (opt$st==TRUE){
sub.folder <- "/NoStress"
sub.file <- ".NOStress"
} else {
sub.folder <- ""
sub.file <- ""
}
#load data
setwd("../")
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".EpiSubcluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda"))
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".StSubcluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda"))
#create folder structure
if (!dir.exists(paste0("./analysis/",opt$g))){
dir.create(paste0("./analysis/",opt$g))
}
if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){
dir.create(paste0("./analysis/",opt$g,"/global"))
}
if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder))){
dir.create(paste0("./analysis/",opt$g,"/global",sub.folder))
}
if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){
dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))
}
#load resolution
sc10x.Group.Epi <- SetAllIdent(object=sc10x.Group.Epi,id=paste0("res",opt$r))
sc10x.Group.St <- SetAllIdent(object=sc10x.Group.St,id=paste0("res",opt$r))
for (x in c("Epi","St")){
#downsample if number given
if (x=="Epi"){
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group.Epi@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group.Epi@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group.Epi@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group.Epi@ident)))
labels <- labels[rnd]
rm(rnd)
} else {
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group.St@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group.St@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group.St@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group.St@ident)))
labels <- labels[rnd]
rm(rnd)
}}
#QuSAGE with LGEA epi markers
Number.Clusters <- length(unique(sc10x.Group.Epi@ident))
C <- list()
cC <- list()
for (i in 1:Number.Clusters){
t <- labels
t[!(t %in% paste0("Cluster_",i))] <- "REST"
C[i] <- list(i=t)
rm(t)
cC[i] <- paste0("Cluster_",i,"-REST")
}
rm(labels)
gene.set2 <- read_csv("./genesets/Basal cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hBC"
gene.set <- c(gene.set2)
gene.set2 <- read_csv("./genesets/Normal AT2 cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hAT2"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- read_csv("./genesets/Club_Goblet cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hCGC"
gene.set <- c(gene.set,gene.set2)
rm(gene.set2)
#gene.set.c2 <- read.gmt("./genesets/c2.all.v6.1.symbols.gmt")
#gene.set.c5 <- read.gmt("./genesets/c5.all.v6.1.symbols.gmt")
gene.set.BioGPS <- read.gmt("./genesets/gene_set_library_up_crisp-BioGPS.gmt")
gene.set.TISSUES <- read.gmt("./genesets/gene_set_library_crisp-TISSUES.gmt")
for (i in 1:Number.Clusters){
assign(paste0("qs.results.Epi.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
#assign(paste0("qs.results.c2.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.c2))
#assign(paste0("qs.results.c5.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.c5))
assign(paste0("qs.results.Epi.BioGPS.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.BioGPS))
assign(paste0("qs.results.Epi.TISSUES.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.TISSUES))
}
rm(eset)
rm(C)
rm(cC)
#generate correlation figures
max.x.rg <- 0
min.x.rg <- 0
max.y.rg <- 0
for (i in 1:Number.Clusters){
qs <- get(paste0("qs.results.",i))
if (max(qs$path.mean)>max.x.rg){
max.x.rg <- max(qs$path.mean)
if (min(qs$path.mean)<min.x.rg){
min.x.rg <- min(qs$path.mean)
}}
if (max(qs$path.PDF)>max.y.rg){
max.y.rg <- max(qs$path.PDF)
}}
for (i in 1:4){
qs <- get(paste0("qs.results.Epi.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_Epi_",i,".LGEA.eps"))
plotDensityCurves(qs,path.index=1:3,col=1:Num,main=paste0("Cluster ",i),xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation")
legend("topright",legend=c("hBC","hAT2","hCGC"),lty=1,col=1:3,lwd=1,cex=1,ncol=1)
dev.off()
qs <- get(paste0("qs.results.Epi.BioGPS.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_Epi_",i,".CI.BioGPS.eps"))
plotCIs(qs,path.index=1:numPathways(qs))
dev.off()
tab <- qsTable(qs,number=numPathways(qs),sort.by="logFC")
tab <- tab[tab$FDR<0.05,]
tab <- tab[order(-tab[,2]),]
write.table(tab,file=paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_Epi_",i,".Table.BioGPS.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
qs <- get(paste0("qs.results.Epi.TISSUES.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_Epi_",i,".CI.TISSUES.eps"))
plotCIs(qs,path.index=1:numPathways(qs))
dev.off()
tab <- qsTable(qs,number=numPathways(qs),sort.by="logFC")
tab <- tab[tab$FDR<0.05,]
tab <- tab[order(-tab[,2]),]
write.table(tab,file=paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_Epi_",i,".Table.TISSUES.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
}
#QuSAGE with LGEA epi markers
Number.Clusters <- length(unique(sc10x.Group.St@ident))
C <- list()
cC <- list()
for (i in 1:Number.Clusters){
t <- labels
t[!(t %in% paste0("Cluster_",i))] <- "REST"
C[i] <- list(i=t)
rm(t)
cC[i] <- paste0("Cluster_",i,"-REST")
}
rm(labels)
gene.set2 <- read_csv("./genesets/Basal cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hBC"
gene.set <- c(gene.set2)
gene.set2 <- read_csv("./genesets/Normal AT2 cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hAT2"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- read_csv("./genesets/Club_Goblet cells-signature-genes.csv")
gene.set2 <- gene.set2[,2]
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "hCGC"
gene.set <- c(gene.set,gene.set2)
rm(gene.set2)
#gene.set.c2 <- read.gmt("./genesets/c2.all.v6.1.symbols.gmt")
#gene.set.c5 <- read.gmt("./genesets/c5.all.v6.1.symbols.gmt")
gene.set.BioGPS <- read.gmt("./genesets/gene_set_library_up_crisp-BioGPS.gmt")
gene.set.TISSUES <- read.gmt("./genesets/gene_set_library_crisp-TISSUES.gmt")
for (i in 1:Number.Clusters){
assign(paste0("qs.results.St.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
#assign(paste0("qs.results.c2.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.c2))
#assign(paste0("qs.results.c5.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.c5))
assign(paste0("qs.results.St.BioGPS.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.BioGPS))
assign(paste0("qs.results.St.TISSUES.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.TISSUES))
}
rm(eset)
rm(C)
rm(cC)
#generate correlation figures
max.x.rg <- 0
min.x.rg <- 0
max.y.rg <- 0
for (i in 1:Number.Clusters){
qs <- get(paste0("qs.results.",i))
if (max(qs$path.mean)>max.x.rg){
max.x.rg <- max(qs$path.mean)
if (min(qs$path.mean)<min.x.rg){
min.x.rg <- min(qs$path.mean)
}}
if (max(qs$path.PDF)>max.y.rg){
max.y.rg <- max(qs$path.PDF)
}}
for (i in 1:4){
qs <- get(paste0("qs.results.St.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".LGEA.eps"))
plotDensityCurves(qs,path.index=1:3,col=1:3,main=paste0("Cluster ",i),xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation")
legend("topright",legend=c("hBC","hAT2","hCGC"),lty=1,col=1:3,lwd=1,cex=1,ncol=1)
dev.off()
qs <- get(paste0("qs.results.St.BioGPS.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".CI.BioGPS.eps"))
plotCIs(qs,path.index=1:numPathways(qs))
dev.off()
tab <- qsTable(qs,number=numPathways(qs),sort.by="logFC")
tab <- tab[tab$FDR<0.05,]
tab <- tab[order(-tab[,2]),]
write.table(tab,file=paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".Table.BioGPS.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
qs <- get(paste0("qs.results.St.TISSUES.",i))
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".CI.TISSUES.eps"))
plotCIs(qs,path.index=1:numPathways(qs))
dev.off()
tab <- qsTable(qs,number=numPathways(qs),sort.by="logFC")
tab <- tab[tab$FDR<0.05,]
tab <- tab[order(-tab[,2]),]
write.table(tab,file=paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".Table.TISSUES.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
}
save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LinSubClust.Rda"))
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