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

Remove outdated (non-SubClust) QuSAGE .R files

parent 459fe2d0
Branches downsampled_fastqs
No related merge requests found
Hide whitespace changes
Inline Side-by-side
r.scripts/sc_QuSAGE.Epi.R deleted 100755 → 0
+ 0
− 306
View file @ 459fe2d0
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=20000,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,".Cumulative.RData"))
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDlineage.Rda"))
sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDlineage"))
rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDlineage"))
#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 <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
#Get Epi clusters
Epi <- table(sc10x.Group@ident,sc10x.Group@meta.data$Lineage)
Epi <- Epi[,"Epi"]
Epi <- Epi[Epi>0]
Epi <- as.integer(names(Epi))
#downsample if number given
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group@ident)))
labels <- labels[rnd]
rm(rnd)
#QuSAGE
Number.Clusters <- length(unique(sc10x.Group@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_delim("./genesets/DEG_BE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "BE"
gene.set <- c(gene.set2)
gene.set2 <- read_delim("./genesets/DEG_LE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "LE"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- read_delim("./genesets/DEG_OE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "OE"
gene.set <- c(gene.set,gene.set2)
rm(gene.set2)
for (i in Epi){
assign(paste0("qs.results.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
}
rm(eset)
rm(C)
rm(cC)
#generate ID table
rm(qs.Cor.table)
for (i in Epi){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
qs <- qs[(qs$FDR < 0.05) & (qs$log.fold.change > 0),1:2]
if (nrow(qs)==0){
qs <- NULL
qs$pathway.name <- "?"
qs$log.fold.change <- 0
qs <- as.data.frame(qs)
}
qs$Cluster <- i
if (i==Epi[1]){
qs.Cor.table <- qs
} else {
qs.Cor.table <- rbind(qs.Cor.table,qs)
}}
rownames(qs.Cor.table) <- NULL
rm(qs)
Epi.qs.Cor.table <- qs.Cor.table[qs.Cor.table[,3]==1,][which.max(qs.Cor.table[qs.Cor.table[,3]==1,][,2]),]
for (i in 2:Number.Clusters){
Epi.qs.Cor.table <- rbind(Epi.qs.Cor.table,qs.Cor.table[qs.Cor.table[,3]==i,][which.max(qs.Cor.table[qs.Cor.table[,3]==i,][,2]),])
}
#generate correlation figures
max.x.rg <- 0
min.x.rg <- 0
max.y.rg <- 0
for (i in Epi){
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 (j in 1:3){
if (j==1){
pop <- "BE"
} else if (j==2) {
pop <- "LE"
} else if (j==3) {
pop <- "OE"
}
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".eps"))
for (i in Epi){
qs <- get(paste0("qs.results.",i))
if (i==Epi[1]){
plotDensityCurves(qs,path.index=j,col=i,main=pop,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")
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=i)
}}
leg <- paste0("Cluster ",Epi)
legend("topright",legend=leg,lty=1,col=Epi,lwd=1,cex=1,ncol=2)
dev.off()
}
#save EpiPop labels
sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Lineage")
St.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="St"])
UK.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="?"])
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=St.cells,ident.use="St")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=UK.cells,ident.use="?")
merge.cluster <- NULL
for (i in Epi){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
merge.cluster <- c(merge.cluster,as.character(factor(qs[which.max(qs$log.fold.change),1])))
}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=Epi,to=merge.cluster)
if ("?" %in% names(summary(sc10x.Group@ident))){
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("BE","LE","OE","St","?"))
} else {
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("BE","LE","OE","St"))
}
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="EpiPop")
#generate EpiPop tSNE
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_EpiPop.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()
#label EpiSubPops and generate tSNE
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
Ident <- table(sc10x.Group@ident,sc10x.Group@meta.data$EpiPop)
tryCatch({
UK <- Ident[,"?"]
UK <- UK[UK>0]
UK <- as.integer(names(UK))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
St <- Ident[,"St"]
St <- St[St>0]
St <- as.integer(names(St))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
BE <- Ident[,"BE"]
BE <- BE[BE>0]
BE <- as.integer(names(BE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
LE <- Ident[,"LE"]
LE <- LE[LE>0]
LE <- as.integer(names(LE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
OE <- Ident[,"OE"]
OE <- OE[OE>0]
OE <- as.integer(names(OE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
merge.cluster <- NULL
for (i in 1:Number.Clusters){
if (exists("UK")==TRUE){
for (j in 1:length(UK)){
if (UK[j]==i){
if (length(UK)>1){
merge.cluster <- c(merge.cluster,paste0("?",j))
} else {
merge.cluster <- c(merge.cluster,"?")
}}}}
if (exists("BE")==TRUE){
for (j in 1:length(BE)){
if (BE[j]==i){
if (length(BE)>1){
merge.cluster <- c(merge.cluster,paste0("BE",j))
} else {
merge.cluster <- c(merge.cluster,"BE")
}}}}
if (exists("LE")==TRUE){
for (j in 1:length(LE)){
if (LE[j]==i){
if (length(LE)>1){
merge.cluster <- c(merge.cluster,paste0("LE",j))
} else {
merge.cluster <- c(merge.cluster,"LE")
}}}}
if (exists("OE")==TRUE){
for (j in 1:length(OE)){
if (OE[j]==i){
if (length(OE)>1){
merge.cluster <- c(merge.cluster,paste0("OE",j))
} else {
merge.cluster <- c(merge.cluster,"OE")
}}}}
if (exists("St")==TRUE){
for (j in 1:length(St)){
if (St[j]==i){
if (length(St)>1){
merge.cluster <- c(merge.cluster,paste0("St",j))
} else {
merge.cluster <- c(merge.cluster,"St")
}}}}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster)
clust.order <- c(grep("BE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("LE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("OE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("St",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="EpiSubPop")
rm(clust.order)
rm(Ident)
rm(St)
rm(BE)
rm(LE)
rm(OE)
rm(merge.cluster)
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_EpiSubPop.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()
#make generic variables of run specific variables
rm(plot)
rm(UK)
rm(UK.cells)
rm(St.cells)
rm(Epi)
rm(qs)
rm(gene.set)
rm(i)
rm(j)
rm(leg)
rm(max.x.rg)
rm(min.x.rg)
rm(max.y.rg)
rm(pop)
rm(Number.Clusters)
Group.List <- c(Group.List,"EpiPop","EpiSubPop")
assign(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi"),sc10x.Group)
rm(sc10x.Group)
opt.epipop <- opt
rm(opt)
#save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files
save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.EpiPop.Rda"))
rm(list=ls(pattern="^qs."))
save(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.epipop$g,".cluster",sub.file,".IDepi.Rda"))
rm(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi"))
rm(sub.folder)
rm(sub.file)
save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData"))
r.scripts/sc_QuSAGE.LGEA.Epi.R deleted 100755 → 0
+ 0
− 310
View file @ 459fe2d0
gc()
.libPaths("/home2/ghenry/R/x86_64-pc-linux-gnu-library/3.3")
library(optparse)
library(Seurat)
library(readr)
library(readxl)
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=20000,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,".Cumulative.RData"))
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne.Rda"))
sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne"))
rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne"))
#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 <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
#Get Epi clusters
Epi <- table(sc10x.Group@ident,sc10x.Group@meta.data$Lineage)
Epi <- Epi[,"Epi"]
Epi <- Epi[Epi>0]
Epi <- as.integer(names(Epi))
#downsample if number given
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group@ident)))
labels <- labels[rnd]
rm(rnd)
#QuSAGE with LGEA epi markers
Number.Clusters <- length(unique(sc10x.Group@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)
for (i in Epi){
assign(paste0("qs.results.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
}
rm(eset)
rm(C)
rm(cC)
#generate ID table
rm(qs.Cor.table)
for (i in Epi){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
qs <- qs[(qs$FDR < 0.05) & (qs$log.fold.change > 0),1:2]
if (nrow(qs)==0){
qs <- NULL
qs$pathway.name <- "?"
qs$log.fold.change <- 0
}
qs$Cluster <- i
if (i==Epi[1]){
qs.Cor.table <- qs
} else {
qs.Cor.table <- rbind(qs.Cor.table,qs)
}
}
rownames(qs.Cor.table) <- NULL
rm(qs)
#generate correlation figures
max.x.rg <- 0
min.x.rg <- 0
max.y.rg <- 0
for (i in Epi){
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 (j in 1:3){
if (j==1){
pop <- "hBC"
} else if (j==2) {
pop <- "hAT2"
} else if (j==3) {
pop <- "hCGC"
}
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".eps"))
for (i in Epi){
qs <- get(paste0("qs.results.",i))
if (i==Epi[1]){
plotDensityCurves(qs,path.index=j,col=i,main=pop,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")
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=i)
}}
leg <- paste0("Cluster ",Epi)
legend("topright",legend=leg,lty=1,col=Epi,lwd=1,cex=1,ncol=2)
dev.off()
}
#save LGEA Epi labels
sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Lineage")
tryCatch({
St.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="St"])
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
UK.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="?"])
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=St.cells,ident.use="St")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=UK.cells,ident.use="?")
merge.cluster <- NULL
for (i in Epi){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
merge.cluster <- c(merge.cluster,as.character(factor(qs[which.max(qs[qs[qs$log.fold.change>0,]$p.Value<0.05,]$log.fold.change),1])))
}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=Epi,to=merge.cluster)
clust.order <- c(grep("hBC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("hAT2",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("hCGC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("St",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGEA.Epi")
rm(clust.order)
#generate EpiPop tSNE
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_LGEA.Epi.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()
#label LGEA EpiSubPops and generate tSNE
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
Ident <- table(sc10x.Group@ident,sc10x.Group@meta.data$LGEA.Epi)
tryCatch({
UK <- Ident[,"?"]
UK <- UK[UK>0]
UK <- as.integer(names(UK))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
St <- Ident[,"St"]
St <- St[St>0]
St <- as.integer(names(St))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
hBC <- Ident[,"hBC"]
hBC <- hBC[hBC>0]
hBC <- as.integer(names(hBC))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
hAT2 <- Ident[,"hAT2"]
hAT2 <- hAT2[hAT2>0]
hAT2 <- as.integer(names(hAT2))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
hCGC <- Ident[,"hCGC"]
hCGC <- hCGC[hCGC>0]
hCGC <- as.integer(names(hCGC))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
merge.cluster <- NULL
for (i in 1:Number.Clusters){
if (exists("UK")==TRUE){
for (j in 1:length(UK)){
if (UK[j]==i){
if (length(UK)>1){
merge.cluster <- c(merge.cluster,paste0("?"))
} else {
merge.cluster <- c(merge.cluster,"?")
}}}}
if (exists("St")==TRUE){
for (j in 1:length(St)){
if (St[j]==i){
if (length(St)>1){
merge.cluster <- c(merge.cluster,"St")
} else {
merge.cluster <- c(merge.cluster,"St")
}}}}
if (exists("hBC")==TRUE){
for (j in 1:length(hBC)){
if (hBC[j]==i){
if (length(hBC)>1){
merge.cluster <- c(merge.cluster,paste0("hBC",j))
} else {
merge.cluster <- c(merge.cluster,"hBC")
}}}}
if (exists("hAT2")==TRUE){
for (j in 1:length(hAT2)){
if (hAT2[j]==i){
if (length(hAT2)>1){
merge.cluster <- c(merge.cluster,paste0("hAT2",j))
} else {
merge.cluster <- c(merge.cluster,"hAT2")
}}}}
if (exists("hCGC")==TRUE){
for (j in 1:length(hCGC)){
if (hCGC[j]==i){
if (length(hCGC)>1){
merge.cluster <- c(merge.cluster,paste0("hCGC",j))
} else {
merge.cluster <- c(merge.cluster,"hCGC")
}}}}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster)
clust.order <- c(grep("hBC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("hAT2",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("hCGC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("St",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGEA.EpiSubPop")
rm(clust.order)
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_LGEA.EpiSubPop.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()
rm(Ident)
rm(UK)
rm(Epi)
rm(St)
rm(hBC)
rm(hAT2)
rm(hCGC)
rm(St.cells)
rm(UK.cells)
rm(merge.cluster)
#make generic variables of run specific variables
rm(plot)
rm(Epi)
rm(qs)
rm(gene.set)
rm(i)
rm(j)
rm(leg)
rm(max.x.rg)
rm(min.x.rg)
rm(max.y.rg)
rm(pop)
rm(Number.Clusters)
Group.List <- c(Group.List,"LGEA.Epi","LGEA.EpiSubPop")
assign(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"),sc10x.Group)
rm(sc10x.Group)
opt.lgeaepipop <- opt
rm(opt)
#save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files
save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.EpiPop.Rda"))
rm(list=ls(pattern="^qs."))
save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda"))
rm(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"))
rm(sub.folder)
rm(sub.file)
save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData"))
r.scripts/sc_QuSAGE.LGEA.St.R deleted 100755 → 0
+ 0
− 383
View file @ 459fe2d0
gc()
.libPaths("/home2/ghenry/R/x86_64-pc-linux-gnu-library/3.3")
library(optparse)
library(Seurat)
library(readr)
library(qusage)
library(readxl)
#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=20000,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,".Cumulative.RData"))
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda"))
sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"))
rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"))
#create folder structure
if (!dir.exists(paste0("./anaysis/",opt$g))){
dir.create(paste0("./analysis/",opt$g))
}
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 <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
#Get St clusters
St <- table(sc10x.Group@ident,sc10x.Group@meta.data$Lineage)
if ("NewColumn" %in% colnames(St)){
St <- St[,c("NewColumn","St")]
St <- rowSums(St)
} else {
St <- St[,"St"]
}
St <- St[St>0]
St <- as.integer(names(St))
#downsample if number given
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group@ident)))
labels <- labels[rnd]
rm(rnd)
#QuSAGE using LGEA human otholog St markers
Number.Clusters <- length(unique(sc10x.Group@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.set3 <- read_excel("./genesets/journal.pcbi.1004575.s026.XLSX",skip=2)
gene.set3 <- gene.set3[,c(3,9)]
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C1 : Proliferative Fibroblast",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "ProlFib"
gene.set <- c(gene.set2)
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C2: Myofibroblast / Smooth Muscle-like",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "MyoFibSM"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C3: Pericyte",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "Peri"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C5: Matrix Fibroblast",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "MatFib"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C7: Endothelial Cells",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "EndoC"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C8: Myeloid / Immune Cells",]
gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1])
names(gene.set2) <- "MylImm"
gene.set <- c(gene.set,gene.set2)
rm(gene.set2)
rm(gene.set3)
for (i in 1:Number.Clusters){
assign(paste0("qs.results.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
}
rm(eset)
rm(C)
rm(cC)
#generate ID table
UK.clust <- as.numeric(Lin.qs.Cor.table[Lin.qs.Cor.table[,1]=="?",][,3])
rm(qs.Cor.table)
for (i in c(St,UK.clust)){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
qs <- qs[(qs$FDR < 0.05) & (qs$log.fold.change > 0),1:2]
if (nrow(qs)==0){
qs <- NULL
qs$pathway.name <- "?"
qs$log.fold.change <- 0
}
qs$Cluster <- i
if (i==St[1]){
qs.Cor.table <- qs
} else {
qs.Cor.table <- rbind(qs.Cor.table,qs)
}
}
rownames(qs.Cor.table) <- NULL
rm(qs)
#generate correlation figures
max.x.rg <- 0
min.x.rg <- 0
max.y.rg <- 0
for (i in c(St,UK.clust)){
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 (j in 1:6){
if (j==1){
pop <- "ProlFib"
} else if (j==2) {
pop <- "MyoFibSM"
} else if (j==3) {
pop <- "Peri"
} else if (j==4){
pop <- "MatFib"
} else if (j==5) {
pop <- "EndoC"
} else if (j==6) {
pop <- "MylImm"
}
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".eps"))
for (i in St){
qs <- get(paste0("qs.results.",i))
if (i==St[1]){
plotDensityCurves(qs,path.index=j,col=i,main=pop,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")
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=i)
}}
leg <- paste0("Cluster ",St)
legend("topright",legend=leg,lty=1,col=St,lwd=1,cex=1,ncol=2)
dev.off()
}
#save LGEA St labels
merge.cluster <- NULL
for (i in c(St,UK.clust)){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
if (as.numeric(as.character(factor(qs[which.max(qs$log.fold.change),2])))>0){
merge.cluster <- c(merge.cluster,as.character(factor(qs[which.max(qs$log.fold.change),1])))
} else {
merge.cluster <- c(merge.cluster,"?St")
}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=c(St,UK.clust),to=merge.cluster)
ProlFib.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="ProlFib"])
MyoFibSM.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="MyoFibSM"])
Peri.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="Peri"])
MatFib.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="MatFib"])
EndoC.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="EndoC"])
MylImm.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="MylImm"])
UK.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="?St"])
sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Lineage")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=ProlFib.cells,ident.use="ProlFib")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=MyoFibSM.cells,ident.use="MyoFibSM")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=Peri.cells,ident.use="Peri")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=MatFib.cells,ident.use="MatFib")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=EndoC.cells,ident.use="EndoC")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=MylImm.cells,ident.use="MylImm")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=UK.cells,ident.use="?St")
clust.order <- c(grep("Epi",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("ProlFib",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MatFib",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MyoFibSM",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MylImm",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("EndoC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("Peri",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGEA.St")
rm(clust.order)
rm(ProlFib.cells)
rm(MyoFibSM.cells)
rm(Peri.cells)
rm(MatFib.cells)
rm(EndoC.cells)
rm(MylImm.cells)
rm(UK.cells)
#generate LGEA StPops tSNE
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_LGEA.St.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()
#label LGEA StSubPops and generate tSNE
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
Ident <- table(sc10x.Group@ident,sc10x.Group@meta.data$LGEA.St)
tryCatch({
UK <- Ident[,"?"]
UK <- UK[UK>0]
UK <- as.integer(names(UK))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
Epi <- Ident[,"Epi"]
Epi <- Epi[Epi>0]
Epi <- as.integer(names(Epi))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
ProlFib <- Ident[,"ProlFib"]
ProlFib <- ProlFib[ProlFib>0]
ProlFib <- as.integer(names(ProlFib))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
MyoFibSM <- Ident[,"MyoFibSM"]
MyoFibSM <- MyoFibSM[MyoFibSM>0]
MyoFibSM <- as.integer(names(MyoFibSM))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
Peri <- Ident[,"Peri"]
Peri <- Peri[Peri>0]
Peri <- as.integer(names(Peri))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
MatFib <- Ident[,"MatFib"]
MatFib <- MatFib[MatFib>0]
MatFib <- as.integer(names(MatFib))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
EndoC <- Ident[,"EndoC"]
EndoC <- EndoC[EndoC>0]
EndoC <- as.integer(names(EndoC))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
MylImm <- Ident[,"MylImm"]
MylImm <- MylImm[MylImm>0]
MylImm <- as.integer(names(MylImm))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
merge.cluster <- NULL
for (i in 1:Number.Clusters){
if (exists("UK")==TRUE){
for (j in 1:length(UK)){
if (UK[j]==i){
if (length(UK)>1){
merge.cluster <- c(merge.cluster,paste0("?",j))
} else {
merge.cluster <- c(merge.cluster,"?")
}}}}
if (exists("Epi")==TRUE){
for (j in 1:length(Epi)){
if (Epi[j]==i){
if (length(Epi)>1){
merge.cluster <- c(merge.cluster,"Epi")
} else {
merge.cluster <- c(merge.cluster,"Epi")
}}}}
if (exists("ProlFib")==TRUE){
for (j in 1:length(ProlFib)){
if (ProlFib[j]==i){
if (length(ProlFib)>1){
merge.cluster <- c(merge.cluster,paste0("ProlFib",j))
} else {
merge.cluster <- c(merge.cluster,"ProlFib")
}}}}
if (exists("MyoFibSM")==TRUE){
for (j in 1:length(MyoFibSM)){
if (MyoFibSM[j]==i){
if (length(MyoFibSM)>1){
merge.cluster <- c(merge.cluster,paste0("MyoFibSM",j))
} else {
merge.cluster <- c(merge.cluster,"MyoFibSM")
}}}}
if (exists("Peri")==TRUE){
for (j in 1:length(Peri)){
if (Peri[j]==i){
if (length(Peri)>1){
merge.cluster <- c(merge.cluster,paste0("Peri",j))
} else {
merge.cluster <- c(merge.cluster,"Peri")
}}}}
if (exists("MatFib")==TRUE){
for (j in 1:length(MatFib)){
if (MatFib[j]==i){
if (length(MatFib)>1){
merge.cluster <- c(merge.cluster,paste0("MatFib",j))
} else {
merge.cluster <- c(merge.cluster,"MatFib")
}}}}
if (exists("EndoC")==TRUE){
for (j in 1:length(EndoC)){
if (EndoC[j]==i){
if (length(EndoC)>1){
merge.cluster <- c(merge.cluster,paste0("EndoC",j))
} else {
merge.cluster <- c(merge.cluster,"EndoC")
}}}}
if (exists("MylImm")==TRUE){
for (j in 1:length(MylImm)){
if (MylImm[j]==i){
if (length(MylImm)>1){
merge.cluster <- c(merge.cluster,paste0("MylImm",j))
} else {
merge.cluster <- c(merge.cluster,"MylImm")
}}}}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster)
clust.order <- c(grep("Epi",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("ProlFib",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MatFib",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MyoFibSM",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("MylImm",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("EndoC",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("Peri",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGES.StSubPop")
rm(clust.order)
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_LGES.StSubPop.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()
rm(Ident)
rm(Epi)
rm(ProlFib)
rm(MyoFibSM)
rm(Peri)
rm(MatFib)
rm(EndoC)
rm(MylImm)
rm(UK)
rm(merge.cluster)
#make generic variables of run specific variables
rm(UK.clust)
rm(plot)
rm(St)
rm(qs)
rm(gene.set)
rm(i)
rm(j)
rm(leg)
rm(max.x.rg)
rm(min.x.rg)
rm(max.y.rg)
rm(pop)
rm(Number.Clusters)
Group.List <- c(Group.List,"LGEA.St","LGEA.StSubPop")
assign(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"),sc10x.Group)
rm(sc10x.Group)
opt.lgeastpop <- opt
rm(opt)
#save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files
save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.StPop.Rda"))
rm(list=ls(pattern="^qs."))
save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda"))
rm(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"))
rm(sub.folder)
rm(sub.file)
save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData"))
r.scripts/sc_QuSAGE.St.R deleted 100755 → 0
+ 0
− 337
View file @ 459fe2d0
.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=20000,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,".Cumulative.RData"))
load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi.Rda"))
sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi"))
rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi"))
#create folder structure
if (!dir.exists(paste0("./analysis/",opt$g))){
dir.create(paste0("./analysis/",opt$g))
}
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 <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
#Get St clusters
St <- table(sc10x.Group@ident,sc10x.Group@meta.data$Lineage)
St <- St[,"St"]
St <- St[St>0]
St <- as.integer(names(St))
#downsample if number given
if (!is.na(opt$s)){
rnd <- sample(1:ncol(sc10x.Group@data),opt$s)
} else {
rnd <- 1:ncol(sc10x.Group@data)
}
eset <- as.data.frame(as.matrix(sc10x.Group@data))
eset <- eset[,rnd]
labels <- paste0("Cluster_",as.vector(factor(sc10x.Group@ident)))
labels <- labels[rnd]
rm(rnd)
#QuSAGE using GO_ENDOTHELIAL_CELL_DIFFERENTIATION,GO_SMOOTH_MUSCLE_CELL_DIFFERENTIATION,GO_REGULATION_OF_FIBROBLAST_PROLIFERATION
Number.Clusters <- length(unique(sc10x.Group@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_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set2 <- as.list(gene.set2[-1,])
names(gene.set2) <- "Endo"
gene.set <- c(gene.set2)
gene.set2 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set2 <- as.list(gene.set2[-1,])
names(gene.set2) <- "SM"
gene.set <- c(gene.set,gene.set2)
gene.set2 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set2 <- as.list(gene.set2[-1,])
names(gene.set2) <- "Fib"
gene.set <- c(gene.set,gene.set2)
rm(gene.set2)
for (i in 1:Number.Clusters){
assign(paste0("qs.results.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set))
}
rm(eset)
rm(C)
rm(cC)
#generate ID table
UK.clust <- as.numeric(Lin.qs.Cor.table[Lin.qs.Cor.table[,1]=="?",][,3])
rm(qs.Cor.table)
for (i in c(St,UK.clust)){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
qs <- qs[(qs$FDR < 0.05) & (qs$log.fold.change > 0),1:2]
if (nrow(qs)==0){
qs <- NULL
qs$pathway.name <- "?"
qs$log.fold.change <- 0
}
qs$Cluster <- i
if (i==St[1]){
qs.Cor.table <- qs
} else {
qs.Cor.table <- rbind(qs.Cor.table,qs)
}}
rownames(qs.Cor.table) <- NULL
rm(qs)
St.qs.Cor.table <- qs.Cor.table[qs.Cor.table[,3]==1,][which.max(qs.Cor.table[qs.Cor.table[,3]==1,][,2]),]
for (i in 2:Number.Clusters){
St.qs.Cor.table <- rbind(St.qs.Cor.table,qs.Cor.table[qs.Cor.table[,3]==i,][which.max(qs.Cor.table[qs.Cor.table[,3]==i,][,2]),])
}
#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 (j in 1:3){
if (j==1){
pop <- "Endo"
} else if (j==2) {
pop <- "SM"
} else if (j==3) {
pop <- "Fib"
}
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".eps"))
for (i in St){
qs <- get(paste0("qs.results.",i))
if (i==St[1]){
plotDensityCurves(qs,path.index=j,col=i,main=pop,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")
} else {
plotDensityCurves(qs,path.index=j,add=TRUE,col=i)
}}
leg <- paste0("Cluster ",St)
legend("topright",legend=leg,lty=1,col=St,lwd=1,cex=1,ncol=2)
dev.off()
}
#save StPop labels
merge.cluster <- NULL
for (i in c(St,UK.clust)){
qs <- qsTable(get(paste0("qs.results.",i)),number=length(gene.set))
if (as.numeric(as.character(factor(qs[which.max(qs$log.fold.change),2])))>0){
merge.cluster <- c(merge.cluster,as.character(factor(qs[which.max(qs$log.fold.change),1])))
} else {
merge.cluster <- c(merge.cluster,"?St")
}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=c(St,UK.clust),to=merge.cluster)
Endo.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="Endo"])
SM.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="SM"])
Fib.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="Fib"])
UK.cells <- names(sc10x.Group@ident[sc10x.Group@ident=="?St"])
sc10x.Group <- SetAllIdent(object=sc10x.Group,id="EpiPop")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=Endo.cells,ident.use="Endo")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=SM.cells,ident.use="SM")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=Fib.cells,ident.use="Fib")
sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=UK.cells,ident.use="?St")
if ("?St" %in% names(summary(sc10x.Group@ident))){
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("BE","LE","OE","Fib","SM","Endo","?St"))
} else {
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("BE","LE","OE","Fib","SM","Endo"))
}
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="ALLPop")
rm(Endo.cells)
rm(SM.cells)
rm(Fib.cells)
rm(UK.cells)
#generate ALLPop tSNE
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_ALLPop.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()
#label StSubPops and generate tSNE
sc10x.Group <- SetAllIdent(object=sc10x.Group,id=paste0("res",opt$r))
Ident <- table(sc10x.Group@ident,sc10x.Group@meta.data$ALLPop)
tryCatch({
BE <- Ident[,"BE"]
BE <- BE[BE>0]
BE <- as.integer(names(BE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
LE <- Ident[,"LE"]
LE <- LE[LE>0]
LE <- as.integer(names(LE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
OE <- Ident[,"OE"]
OE <- OE[OE>0]
OE <- as.integer(names(OE))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
Endo <- Ident[,"Endo"]
Endo <- Endo[Endo>0]
Endo <- as.integer(names(Endo))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
SM <- Ident[,"SM"]
SM <- SM[SM>0]
SM <- as.integer(names(SM))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
Fib <- Ident[,"Fib"]
Fib <- Fib[Fib>0]
Fib <- as.integer(names(Fib))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
tryCatch({
UKSt <- Ident[,"?St"]
UKSt <- UKSt[UKSt>0]
UKSt <- as.integer(names(UKSt))
},error=function(e){cat("ERROR : ",conditionMessage(e),"/\n")})
merge.cluster <- NULL
for (i in 1:Number.Clusters){
if (exists("BE")==TRUE){
for (j in 1:length(BE)){
if (BE[j]==i){
if (length(BE)>1){
merge.cluster <- c(merge.cluster,paste0("BE",j))
} else {
merge.cluster <- c(merge.cluster,"BE")
}}}}
if (exists("LE")==TRUE){
for (j in 1:length(LE)){
if (LE[j]==i){
if (length(LE)>1){
merge.cluster <- c(merge.cluster,paste0("LE",j))
} else {
merge.cluster <- c(merge.cluster,"LE")
}}}}
if (exists("OE")==TRUE){
for (j in 1:length(OE)){
if (OE[j]==i){
if (length(OE)>1){
merge.cluster <- c(merge.cluster,paste0("OE",j))
} else {
merge.cluster <- c(merge.cluster,"OE")
}}}}
if (exists("Endo")==TRUE){
for (j in 1:length(Endo)){
if (Endo[j]==i){
if (length(Endo)>1){
merge.cluster <- c(merge.cluster,paste0("Endo",j))
} else {
merge.cluster <- c(merge.cluster,"Endo")
}}}}
if (exists("SM")==TRUE){
for (j in 1:length(SM)){
if (SM[j]==i){
if (length(SM)>1){
merge.cluster <- c(merge.cluster,paste0("SM",j))
} else {
merge.cluster <- c(merge.cluster,"SM")
}}}}
if (exists("Fib")==TRUE){
for (j in 1:length(Fib)){
if (Fib[j]==i){
if (length(Fib)>1){
merge.cluster <- c(merge.cluster,paste0("Fib",j))
} else {
merge.cluster <- c(merge.cluster,"Fib")
}}}}
if (exists("UKSt")==TRUE){
for (j in 1:length(UKSt)){
if (UKSt[j]==i){
if (length(UKSt)>1){
merge.cluster <- c(merge.cluster,paste0("?St",j))
} else {
merge.cluster <- c(merge.cluster,"?St")
}}}}}
sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster)
clust.order <- c(grep("BE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("LE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("OE",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("Fib",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("SM",as.character(levels(sc10x.Group@ident)),value=TRUE),grep("Endo",as.character(levels(sc10x.Group@ident)),value=TRUE))
clust.order <- c(clust.order,as.character(levels(sc10x.Group@ident))[!(as.character(levels(sc10x.Group@ident)) %in% clust.order)])
sc10x.Group@ident <- factor(sc10x.Group@ident,levels=clust.order)
sc10x.Group <- StashIdent(object=sc10x.Group,save.name="ALLSubPop")
rm(clust.order)
postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_ALLSubPop.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()
rm(Ident)
rm(St)
rm(BE)
rm(LE)
rm(OE)
rm(Endo)
rm(SM)
rm(Fib)
rm(UKSt)
rm(merge.cluster)
#make generic variables of run specific variables
rm(UK.clust)
rm(plot)
rm(St)
rm(qs)
rm(gene.set)
rm(i)
rm(j)
rm(leg)
rm(max.x.rg)
rm(min.x.rg)
rm(max.y.rg)
rm(pop)
rm(Number.Clusters)
Group.List <- c(Group.List,"ALLPop","ALLSubPop")
assign(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st"),sc10x.Group)
rm(sc10x.Group)
opt.stpop <- opt
rm(opt)
#save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files
save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.StPop.Rda"))
rm(list=ls(pattern="^qs."))
save(list=paste0("sc10x.",opt.stpop$g,".cluster",sub.file,".IDepi+st"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.stpop$g,".cluster",sub.file,".IDepi+st.Rda"))
rm(list=paste0("sc10x.",opt.stpop$g,".cluster",sub.file,".IDepi+st"))
rm(sub.folder)
rm(sub.file)
save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData"))
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