Commit 14ae840e authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Finalize stress removal before align

parent e2110a4e
......@@ -64,26 +64,47 @@ counts.cell.filtered.3gene <- results[[4]]
counts.gene.filtered.3gene <- results[[5]]
rm(results)
#results <- scQC(sc10x,sp=opt$s,feature="percent.ribo")
#sc10x <- results[[1]]
#counts.cell.filtered.4ribo <- results[[4]]
#counts.gene.filtered.4ribo <- results[[5]]
#rm(results)
counts.cell.filtered <- counts.cell.filtered.3gene
counts.gene.filtered <- counts.gene.filtered.3gene
# results <- scCellCycle(sc10x,sub=FALSE,sp="hu")
# sc10x <- results[[1]]
# genes.s <- results[[2]]
# genes.g2m <- results[[3]]
# rm(results)
sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))]
pc.use.prestress <- list()
for (i in names(sc10x)){
sc10x.temp <- sc10x[[i]]
sc10x.temp <- SCTransform(sc10x.temp,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
results <- scPC(sc10x.temp,pc=1000,hpc=0.9,file=paste0(i,".pre.stress"),print="2",assay="SCT")
sc10x.temp <- results[[1]]
pc.use.prestress.temp <- results[[2]]
rm(results)
sc10x.temp <- scCluster(sc10x.temp,res=0.5,red="pca",dim=pc.use.prestress.temp,print="2",folder=paste0(i,".pre.stress"),assay="SCT")
sc10x[i] <- sc10x.temp
pc.use.prestress[i] <- pc.use.prestress.temp
}
rm(sc10x.temp)
if (opt$s == "hu"){
genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress <- genes.stress[1]
colnames(genes.stress) <- "scDWS.Stress"
anchor <- c("EGR1","FOS","JUN")
} else if (opt$s == "mu"){
genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/van.den.Brink1.txt",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress <- genes.stress[1]
colnames(genes.stress) <- "van.den.Brink1.Stress"
anchor <- c("Egr1","Fos","Jun")
}
results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt="triangle",anchor=anchor)
sc10x.preStress <- results[[1]]
sc10x <- results[[2]]
rm(results)
counts.cell.destress <- lapply(sc10x,ncol)
rm(anchor)
merges <- NULL
for (i in names(counts.cell.filtered[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))])){
if (counts.cell.filtered[[i]]<750){
for (i in names(counts.cell.destress[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))])){
if (counts.cell.destress[[i]]<750){
merges <- c(merges,i)
}
}
......@@ -94,146 +115,40 @@ if (length(merges)>1){
sc10x$merge <- sc10x.merge
rm(sc10x.merge)
} else if (length(merges)==1){
merges <- c(merges,names(counts.cell.filtered[counts.cell.filtered == min(sapply(counts.cell.filtered[setdiff(as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1])),merges)],min))]))
sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
merges <- c(merges,names(counts.cell.destress[counts.cell.destress == min(sapply(counts.cell.destress[setdiff(as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1])),merges)],min))]))
sc10x.merge <- merge(x=sc10x.med[[merges[1]]],y=sc10x.med[merges[-1]])
sc10x[merges] <- NULL
sc10x$merge <- sc10x.merge
}
rm(merges)
gc()
if (length(sc10x)>1){
sc10x <- scAlign(sc10x)
sc10x@project.name <- project.name
sc10x$ALL <- "ALL"
}
gc()
results <- scPC(sc10x,pc=1000,hpc=0.9,file="pre.stress",print="2")
results <- scPC(sc10x,pc=1000,hpc=0.9,file="ALL",print="2")
sc10x <- results[[1]]
pc.use.prestress <- results[[2]]
rm(results)
sc10x <- scCluster(sc10x,res=0.5,red="pca",dim=pc.use.prestress,print="2",folder="pre.stress")
if (opt$s == "hu"){
genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress <- genes.stress[1]
colnames(genes.stress) <- "scDWS.Stress"
anchor <- c("EGR1","FOS","JUN")
} else if (opt$s == "mu"){
genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/van.den.Brink1.txt",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress <- genes.stress[1]
colnames(genes.stress) <- "van.den.Brink1.Stress"
anchor <- c("Egr1","Fos","Jun")
}
results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt="minerror",anchor=anchor)
#sc10x.preStress <- results[[1]]
sc10x.high <- results[[2]]
rm(results)
counts.cell.destress.high <- as.list(table(sc10x.high$samples))
results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt="triangle",anchor=anchor)
#sc10x.preStress <- results[[1]]
sc10x.med <- results[[2]]
rm(results)
counts.cell.destress.med <- as.list(table(sc10x.med$samples))
results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt="renyi",anchor=anchor)
#sc10x.preStress <- results[[1]]
sc10x.low <- results[[2]]
rm(results)
counts.cell.destress.low <- as.list(table(sc10x.low$samples))
rm(anchor)
#if (opt$p == "muPr" || opt$p == "muPrUr"){
# hto <- as.character(unlist(sc10x.groups[sc10x.groups$HTO==1,1]))
# sc10x.hto.data <- Read10X(data.dir=paste0("./analysis/DATA/10x/",hto,"/filtered_feature_bc_matrix/"))
# sc10x.hto <- CreateSeuratObject(counts=sc10x.hto.data$`Gene Expression`,project=hto)
# sc10x.hto[["HTO"]] <- CreateAssayObject(counts=sc10x.hto.data$`Antibody Capture`)
# rm(sc10x.hto.data)
# sc10x.hto <- NormalizeData(object=sc10x.hto,assay="HTO",normalization.method="CLR")
# sc10x.hto <- HTODemux(sc10x.hto,assay="HTO",positive.quantile=.99)
# sc10x.hto <- subset(sc10x.hto,cells=substr(names(sc10x$samples)[sc10x$samples == paste0(hto,"_GEX")],start=1,stop=16))
# # HTOHeatmap(sc10x.hto)
# # plot1 <- FeatureScatter(sc10x.hto,feature1 = "hto_Anterior-Prostate",feature2 = "hto_Ventral-Prostate",group.by = "HTO_maxID")
# # plot2 <- FeatureScatter(sc10x.hto,feature1 = "hto_Anterior-Prostate",feature2 = "hto_Dorsolateral-Prostate",group.by = "HTO_maxID")
# # plot3 <- FeatureScatter(sc10x.hto,feature1 = "hto_Dorsolateral-Prostate",feature2 = "hto_Ventral-Prostate",group.by = "HTO_maxID")
# # grid.arrange(plot1,plot2,plot3,ncol=1)
# Idents(sc10x) <- "samples"
# Idents(sc10x,cells=paste0(names(sc10x.hto$HTO_maxID)[sc10x.hto$HTO_maxID=="Anterior-Prostate"],"_4")) <- paste0(hto,"_AP")
# Idents(sc10x,cells=paste0(names(sc10x.hto$HTO_maxID)[sc10x.hto$HTO_maxID=="Dorsolateral-Prostate"],"_4")) <- paste0(hto,"_DLP")
# Idents(sc10x,cells=paste0(names(sc10x.hto$HTO_maxID)[sc10x.hto$HTO_maxID=="Ventral-Prostate"],"_4")) <- paste0(hto,"_VP")
# sc10x$HTO_maxID <- Idents(sc10x)
# sc10x$samples_HTO <- Idents(sc10x)
# Idents(sc10x,cells=paste0(names(sc10x.hto$hash.ID)[sc10x.hto$hash.ID=="Anterior-Prostate"],"_4")) <- paste0(hto,"_AP")
# Idents(sc10x,cells=paste0(names(sc10x.hto$hash.ID)[sc10x.hto$hash.ID=="Dorsolateral-Prostate"],"_4")) <- paste0(hto,"_DLP")
# Idents(sc10x,cells=paste0(names(sc10x.hto$hash.ID)[sc10x.hto$hash.ID=="Ventral-Prostate"],"_4")) <- paste0(hto,"_VP")
# Idents(sc10x,cells=paste0(names(sc10x.hto$hash.ID)[sc10x.hto$hash.ID=="Negative"],"_4")) <- paste0(hto,"_negative")
# Idents(sc10x,cells=paste0(names(sc10x.hto$hash.ID)[sc10x.hto$hash.ID=="Doublet"],"_4")) <- paste0(hto,"_doublet")
# sc10x$hash.ID <- Idents(sc10x)
#}
results <- scPC(sc10x.high,pc=1000,hpc=0.9,file="ALL.high",print="2")
sc10x.high <- results[[1]]
pc.use.poststress.high <- results[[2]]
rm(results)
results <- scPC(sc10x.med,pc=1000,hpc=0.9,file="ALL.med",print="2")
sc10x.med <- results[[1]]
pc.use.poststress.med <- results[[2]]
rm(results)
results <- scPC(sc10x.low,pc=1000,hpc=0.9,file="ALL.low",print="2")
sc10x.low <- results[[1]]
pc.use.poststress.low <- results[[2]]
pc.use.poststress <- results[[2]]
rm(results)
res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
sc10x.high <- scCluster(sc10x.high,res=res,red="pca",dim=pc.use.poststress.high,print="2",folder="ALL.high")
sc10x.med <- scCluster(sc10x.med,res=res,red="pca",dim=pc.use.poststress.med,print="2",folder="ALL.med")
sc10x.low <- scCluster(sc10x.low,res=res,red="pca",dim=pc.use.poststress.low,print="2",folder="ALL.low")
#DimPlot(sc10x,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
#DefaultAssay(sc10x) <- "integrated"
#sc10x <- FindClusters(sc10x,resolution=0.1,verbose=FALSE,algorithm=1,group.singletons=FALSE)
#if (opt$p == "muPr" || opt$p == "muPrUr"){
# rm(sc10x.hto.data)
# rm(sc10x.hto)
# Idents(sc10x) <- "samples"
# sc10x.tmp <- subset(sc10x,ident=paste0(hto,"_GEX"))
# postscript(paste0("./analysis/vis/ALL/HTO_HTO_maxID.eps"))
# DimPlot(sc10x.tmp,reduction="umap",group.by="HTO_maxID")
# dev.off()
# postscript(paste0("./analysis/vis/ALL/HTO_hash.ID.eps"))
# DimPlot(sc10x.tmp,reduction="umap",group.by="hash.ID")
# dev.off()
#}sceasy
sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="2",folder="ALL")
DefaultAssay(object=sc10x) <- "SCT"
#scShinyOutput(sc10x,anal="raw")
sc10x@meta.data <- sc10x@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.high@meta.data <- sc10x.high@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.med@meta.data <- sc10x.med@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.low@meta.data <- sc10x.low@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
saveRDS(sc10x,paste0("./analysis/",project.name,"_raw.rds"))
saveRDS(sc10x.high,paste0("./analysis/",project.name,"_raw.high.rds"))
saveRDS(sc10x.med,paste0("./analysis/",project.name,"_raw.med.rds"))
saveRDS(sc10x.low,paste0("./analysis/",project.name,"_raw.low.rds"))
#save(sc10x,file=paste0("./analysis/sc10x.raw.rda"))
save.image(file="./analysis/sc10x.raw.RData")
library(sceasy)
library(reticulate)
use_condaenv('sceasy')
convertFormat(sc10x.high,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.high.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.med,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.med.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.low,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.low.h5ad"),assay="SCT",main_layer="scale.data")
#saveRDS(sc10x,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",project.name,"_raw.rds"))
convertFormat(sc10x,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.h5ad"),assay="SCT",main_layer="scale.data")
saveRDS(sc10x,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",project.name,"_raw.rds"))
......@@ -493,7 +493,7 @@ scCellCycle <- function(sc10x,sub=FALSE,sp="hu"){
}
scPC <- function(sc10x,pc=50,hpc=0.9,file="pre.stress",print="tsne"){
scPC <- function(sc10x,pc=50,hpc=0.9,file="pre.stress",print="tsne",assay="integrated"){
#Scale Seurat object & calculate # of PCs to use
#Inputs:
......@@ -508,7 +508,7 @@ scPC <- function(sc10x,pc=50,hpc=0.9,file="pre.stress",print="tsne"){
#Run PCA
Idents(object=sc10x) <- "ALL"
sc10x <- RunPCA(object=sc10x,npcs=pc,verbose=FALSE,assay="integrated")
sc10x <- RunPCA(object=sc10x,npcs=pc,verbose=FALSE,assay=assay)
#Calculate PCs to use
pc.use <- sc10x[["pca"]]@stdev^2
......@@ -535,7 +535,7 @@ scAlign <- function(sc10x.l){
#sc10x.l[[i]] <- NormalizeData(sc10x.l[[i]],verbose=FALSE)
gc()
#sc10x.l[[i]] <- ScaleData(sc10x.l[[i]],vars.to.regress=c("nFeature_RNA","percent.mito"),verbose = FALSE)
sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,assay="RNA")
gc()
#sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
}
......@@ -556,14 +556,14 @@ scAlign <- function(sc10x.l){
#gc()
gc()
sc10x <- SCTransform(sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,return.only.var.genes=FALSE,assay="RNA")
sc10x <- SCTransform(sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,return.only.var.genes=FALSE,assay="RNA")
gc()
return(sc10x)
}
scCluster <- function(sc10x,res=0.1,red="pca",dim,print="tsne",folder=FALSE){
scCluster <- function(sc10x,res=0.1,red="pca",dim,print="tsne",folder=FALSE,assay="integrated"){
#Cluster Seurat object and produce visualizations
#Inputs:
......@@ -588,11 +588,11 @@ scCluster <- function(sc10x,res=0.1,red="pca",dim,print="tsne",folder=FALSE){
}
DefaultAssay(sc10x) <- "integrated"
DefaultAssay(sc10x) <- assay
#Calculste Vis
sc10x <- RunTSNE(sc10x,dims=1:dim,reduction="pca",assay="integrated")
sc10x <- RunUMAP(sc10x,dims=1:dim,reduction="pca",assay="integrated")
sc10x <- RunTSNE(sc10x,dims=1:dim,reduction="pca",assay=assay)
sc10x <- RunUMAP(sc10x,dims=1:dim,reduction="pca",assay=assay)
sc10x <- FindNeighbors(sc10x,reduction=red,verbose=FALSE)
......@@ -662,7 +662,7 @@ scCluster <- function(sc10x,res=0.1,red="pca",dim,print="tsne",folder=FALSE){
}
scScore <- function(sc10x,score,geneset,cut.pt=0.9,anchor=FALSE){
scScore <- function(sc10x.l,score,geneset,cut.pt=0.9,anchor=FALSE){
#Runs custom PCA analysis for stress ID
#Inputs:
......@@ -684,87 +684,97 @@ scScore <- function(sc10x,score,geneset,cut.pt=0.9,anchor=FALSE){
dir.create(paste0("./analysis/vis/",score))
}
#Score geneset
sc10x <- AddModuleScore(object=sc10x,features=geneset,name=score,assay="SCT")
Idents(object=sc10x) <- paste0(score,"1")
#CDF
cdf <- ecdf(as.numeric(levels(sc10x)))
if (cut.pt == "renyi"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="RenyiEntropy"]
} else if (cut.pt == "triangle"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
sc10x.l.negative <- list()
sc10x.l.positive <- list()
for (i in names(sc10x.l)){
sc10x <- sc10x.l[[i]]
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="Triangle"]
} else if (cut.pt == "minerror"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
#Score geneset
sc10x <- AddModuleScore(object=sc10x,features=geneset,name=score,assay="SCT")
Idents(object=sc10x) <- paste0(score,"1")
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="MinErrorI"]
} else {
cut.x <- quantile(cdf,probs=cut.pt)
cut.x <- unname(cut.x)
}
postscript(paste0("./analysis/score_id/",score,"/CDF_",score,".eps"))
plot(cdf,main=paste0("Cumulative Distribution of ",score," Score"),xlab=paste0(score," Score"),ylab="CDF")
abline(v=cut.x,col="red")
dev.off()
#KDE
postscript(paste0("./analysis/score_id/",score,"/Histo_",score,".eps"))
plot(ggplot(data.frame(Score=as.numeric(levels(sc10x))),aes(x=Score))+geom_histogram(bins=100,aes(y=..density..))+geom_density()+geom_vline(xintercept=cut.x,size=1,color="red")+ggtitle(paste0(score," Score"))+cowplot::theme_cowplot())
dev.off()
Idents(object=sc10x) <- "ALL"
predicate <- paste0(score,"1 >= ",cut.x)
Idents(object=sc10x, cells = rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]] >= cut.x]) <- score
sc10x[[score]] <- Idents(object=sc10x)
Idents(sc10x) <- score
sc10x.negative <- subset(x=sc10x,idents="ALL")
sc10x.positive <- subset(x=sc10x,idents=score)
#Generate plots
postscript(paste0("./analysis/vis/",score,"/3Vis_",score,".eps"))
plot1 <- DimPlot(sc10x,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("ALL")
plot2 <- DimPlot(sc10x.negative,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot3 <- DimPlot(sc10x.positive,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot4 <- DimPlot(sc10x,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Negative")
plot5 <- DimPlot(sc10x.negative,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot6 <- DimPlot(sc10x.positive,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot7 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Positive")
plot8 <- DimPlot(sc10x.negative,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot9 <- DimPlot(sc10x.positive,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
grid.arrange(plot1,plot2,plot3,plot4,plot5,plot6,plot7,plot8,plot9,ncol=3)
dev.off()
#Generate violin plot of gene exvpression
if (anchor!=FALSE){
postscript(paste0("./analysis/score_id/",score,"/Violin_",score,".eps"))
plot <- VlnPlot(object=sc10x,features=anchor,pt.size=0.1,assay="SCT")
plot(plot)
#CDF
cdf <- ecdf(as.numeric(levels(sc10x)))
if (cut.pt == "renyi"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="RenyiEntropy"]
} else if (cut.pt == "triangle"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="Triangle"]
} else if (cut.pt == "minerror"){
h <- hist(data.frame(sc10x[[paste0(score,"1")]])[,paste0(score,"1")],breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
cells.remove <- rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]][,1] < cutoff.temp]
sc10x.temp <- subset(sc10x,cells=cells.remove,invert=TRUE)
thresh <- list()
thresh[["all"]] <- scThresh(list(all=sc10x.temp),feature=paste0(score,"1"),sub=score)
cut.x <- thresh$all$all$threshold[thresh$all$all$method=="MinErrorI"]
} else {
cut.x <- quantile(cdf,probs=cut.pt)
cut.x <- unname(cut.x)
}
postscript(paste0("./analysis/score_id/",score,"/CDF_",i,".",score,".eps"))
plot(cdf,main=paste0("Cumulative Distribution of ",score," Score"),xlab=paste0(score," Score"),ylab="CDF")
abline(v=cut.x,col="red")
dev.off()
#KDE
postscript(paste0("./analysis/score_id/",score,"/Histo_",i,".",score,".eps"))
plot(ggplot(data.frame(Score=as.numeric(levels(sc10x))),aes(x=Score))+geom_histogram(bins=100,aes(y=..density..))+geom_density()+geom_vline(xintercept=cut.x,size=1,color="red")+ggtitle(paste0(score," Score"))+cowplot::theme_cowplot())
dev.off()
Idents(object=sc10x) <- "ALL"
predicate <- paste0(score,"1 >= ",cut.x)
Idents(object=sc10x, cells = rownames(sc10x[[paste0(score,"1")]])[sc10x[[paste0(score,"1")]] >= cut.x]) <- score
sc10x[[score]] <- Idents(object=sc10x)
Idents(sc10x) <- score
sc10x.negative <- subset(x=sc10x,idents="ALL")
sc10x.positive <- subset(x=sc10x,idents=score)
#Generate plots
postscript(paste0("./analysis/vis/",score,"/3Vis_",i,".",score,".eps"))
plot1 <- DimPlot(sc10x,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("ALL")
plot2 <- DimPlot(sc10x.negative,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot3 <- DimPlot(sc10x.positive,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot4 <- DimPlot(sc10x,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Negative")
plot5 <- DimPlot(sc10x.negative,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot6 <- DimPlot(sc10x.positive,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot7 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Positive")
plot8 <- DimPlot(sc10x.negative,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
plot9 <- DimPlot(sc10x.positive,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
grid.arrange(plot1,plot2,plot3,plot4,plot5,plot6,plot7,plot8,plot9,ncol=3)
dev.off()
#Generate violin plot of gene exvpression
if (anchor!=FALSE){
postscript(paste0("./analysis/score_id/",score,"/Violin_",i,".",score,".eps"))
plot <- VlnPlot(object=sc10x,features=anchor,pt.size=0.1,assay="SCT")
plot(plot)
dev.off()
}
sc10x.l[[i]] <- sc10x
sc10x.l.negative[[i]] <- sc10x.negative
sc10x.l.positive[[i]] <- sc10x.positive
}
results <- list(
sc10x <- sc10x,
sc10x.negative <- sc10x.negative,
sc10x.positive <- sc10x.positive
sc10x <- sc10x.l,
sc10x.negative <- sc10x.l.negative,
sc10x.positive <- sc10x.l.positive
)
return(results)
}
......
Markdown is supported
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