Commit e2110a4e authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Mod fitering/align anchors/pca numbers

parent 4ce5c19e
gc()
library(methods)
# library(methods)
library(optparse)
library(Seurat)
library(readr)
library(readxl)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
# library(readxl)
# library(fBasics)
# library(pastecs)
# library(qusage)
# library(RColorBrewer)
# library(monocle)
library(dplyr)
library(viridis)
library(gridExtra)
library(SingleR)
library(sctransform)
# library(SingleR)
# library(sctransform)
library(autothresholdr)
library(ggplot2)
options(bitmapType="cairo")
#options(bitmapType="cairo")
options(future.globals.maxSize= 1000000000000)
option_list=list(
make_option("--p",action="store",default="huPr_Pd",type='character',help="Project Name"),
make_option("--p",action="store",default="huPr_PdPb",type='character',help="Project Name"),
make_option("--s",action="store",default="hu",type='character',help="Species")
)
opt=parse_args(OptionParser(option_list=option_list))
......@@ -107,7 +107,7 @@ if (length(sc10x)>1){
sc10x$ALL <- "ALL"
}
results <- scPC(sc10x,pc=100,hpc=0.9,file="pre.stress",print="2")
results <- scPC(sc10x,pc=1000,hpc=0.9,file="pre.stress",print="2")
sc10x <- results[[1]]
pc.use.prestress <- results[[2]]
rm(results)
......@@ -126,11 +126,24 @@ if (opt$s == "hu"){
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 <- results[[2]]
#sc10x.preStress <- results[[1]]
sc10x.low <- results[[2]]
rm(results)
counts.cell.destress <- as.list(table(sc10x$samples))
counts.cell.destress.low <- as.list(table(sc10x.low$samples))
rm(anchor)
#if (opt$p == "muPr" || opt$p == "muPrUr"){
......@@ -161,13 +174,31 @@ rm(anchor)
# sc10x$hash.ID <- Idents(sc10x)
#}
results <- scPC(sc10x,pc=100,hpc=0.9,file="post.stress",print="2")
sc10x <- results[[1]]
pc.use.poststress <- results[[2]]
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]]
rm(results)
res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="2",folder="ALL")
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)
......@@ -180,20 +211,29 @@ sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="2",folde
# postscript(paste0("./analysis/vis/ALL/HTO_hash.ID.eps"))
# DimPlot(sc10x.tmp,reduction="umap",group.by="hash.ID")
# dev.off()
#}
#}sceasy
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.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"))
saveRDS(sc10x,paste0("./analysis/",project.name,"_raw.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,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"))
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"))
......@@ -540,7 +540,7 @@ scAlign <- function(sc10x.l){
#sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
}
sc10x.features <- SelectIntegrationFeatures(object.list=sc10x.l,nfeatures=3000)
sc10x.features <- SelectIntegrationFeatures(object.list=sc10x.l,nfeatures=5000)
sc10x.l <- PrepSCTIntegration(object.list=sc10x.l,anchor.features=sc10x.features,verbose=FALSE)
sc10x.l <- lapply(sc10x.l,FUN=function(x) { RunPCA(x,features=sc10x.features,verbose=FALSE) })
......@@ -691,15 +691,32 @@ scScore <- function(sc10x,score,geneset,cut.pt=0.9,anchor=FALSE){
#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)
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)
......
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