Commit 7c902417 authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Significant edits

parent 11e5bd18
......@@ -8,10 +8,10 @@
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load python/3.6.4-anaconda
source activate umap
module load python/3.7.x-anaconda
conda activate seurat
module load gdal
module load cairo/1.14.8
module load R/3.6.1-gccmkl
module load R/4.0.2-gccmkl
module load hdf5_18/1.8.17
Rscript ../r.scripts/sc-TissueMapper_RAW.R --p "$1" --s "$2"
......@@ -49,83 +49,84 @@ counts.cell.filtered.2mito <- results[[4]]
counts.gene.filtered.2mito <- results[[5]]
rm(results)
results <- scQC(sc10x,sp=opt$s,feature="percent.ribo")
sc10x <- results[[1]]
counts.cell.filtered.3ribo <- results[[4]]
counts.gene.filtered.3ribo <- results[[5]]
rm(results)
#results <- scQC(sc10x,sp=opt$s,feature="percent.ribo")
#sc10x <- results[[1]]
#counts.cell.filtered.3ribo <- results[[4]]
#counts.gene.filtered.3ribo <- results[[5]]
#rm(results)
results <- scQC(sc10x,sp=opt$s,feature="nCount_RNA")
sc10x <- results[[1]]
counts.cell.filtered.4gene <- results[[4]]
counts.gene.filtered.4gene <- results[[5]]
rm(results)
#results <- scQC(sc10x,sp=opt$s,feature="nCount_RNA")
#sc10x <- results[[1]]
#counts.cell.filtered.3gene <- results[[4]]
#counts.gene.filtered.3gene <- results[[5]]
#rm(results)
counts.cell.filtered <- counts.cell.filtered.4gene
counts.gene.filtered <- counts.gene.filtered.4gene
counts.cell.filtered <- counts.cell.filtered.2mito
counts.gene.filtered <- counts.gene.filtered.2mito
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")
if (ncol(sc10x.temp) > 100) {
pc.calc <- 100
} else if (ncol(sc10x.temp) <= 100) {
pc.calc <- ncol(sc10x.temp)-1
}
results <- scPC(sc10x.temp,pc=pc.calc,hpc=0.9,file=paste0(i,".pre.stress"),print="umap",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="umap",folder=paste0(i,".pre.stress"),assay="SCT")
sc10x[i] <- sc10x.temp
pc.use.prestress[i] <- pc.use.prestress.temp
rm(sc10x.temp)
rm(pc.calc)
rm(pc.use.prestress.temp)
}
rm(i)
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)
# 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")
# if (ncol(sc10x.temp) > 100) {
# pc.calc <- 100
# } else if (ncol(sc10x.temp) <= 100) {
# pc.calc <- ncol(sc10x.temp)-1
# }
# results <- scPC(sc10x.temp,pc=pc.calc,hpc=0.9,file=paste0(i,".pre.stress"),print="umap",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="umap",folder=paste0("preStress/",i),assay="SCT")
# sc10x[i] <- sc10x.temp
# pc.use.prestress[i] <- pc.use.prestress.temp
# rm(sc10x.temp)
# rm(pc.calc)
# rm(pc.use.prestress.temp)
# }
# rm(i)
merges <- NULL
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)
# 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)
# }
# }
# rm(i)
# if (length(merges)>1){
# sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
# sc10x[merges] <- NULL
# sc10x$merge <- sc10x.merge
# rm(sc10x.merge)
# } else if (length(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[[merges[1]]],y=sc10x[merges[-1]])
# sc10x[merges] <- NULL
# sc10x$merge <- sc10x.merge
# }
# rm(merges)
merges <- NULL
for (i in unique(sc10x.groups$Patient[sc10x.groups$Keep==1])){
for (j in sc10x.groups$Samples[sc10x.groups$Keep==1]){
if (sc10x.groups[sc10x.groups$Samples==j,]$Patient==i){
merges <- c(merges,j)
}
}
}
rm(i)
if (length(merges)>1){
if (length(merges)>1){
sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
} else {
sc10x.merge <- sc10x[[merges[1]]]
}
sc10x[merges] <- NULL
sc10x$merge <- sc10x.merge
sc10x[i] <- sc10x.merge
merges <- NULL
rm(sc10x.merge)
} else if (length(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[[merges[1]]],y=sc10x[merges[-1]])
sc10x[merges] <- NULL
sc10x$merge <- sc10x.merge
}
rm(merges)
rm(i)
rm(j)
gc()
if (length(sc10x)>1){
......@@ -152,6 +153,27 @@ 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="umap",folder="ALL")
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(list(ALL=sc10x),score="Stress",geneset=as.list(genes.stress),cut.pt="renyi",anchor=anchor)
sc10x.preStress <- results[[1]]
sc10x <- results[[2]]
#sc10x <- results[[1]]$ALL
rm(results)
#counts.cell.destress <- lapply(sc10x,ncol)
counts.cell.destress <- ncol(sc10x)
rm(anchor)
DefaultAssay(object=sc10x) <- "SCT"
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")]
......
......@@ -546,7 +546,8 @@ 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]],variable.features.n=5000,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,assay="RNA")
#sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],variable.features.n=5000,return.only.var.genes=FALSE,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,assay="RNA")
sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],variable.features.n=10000,return.only.var.genes=FALSE,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
gc()
#sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
}
......@@ -554,20 +555,30 @@ scAlign <- function(sc10x.l){
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,npcs=500 ,verbose=FALSE) })
# pc.calc <- min(unlist(lapply(sc10x.l,FUN=function(x) { ncol(x) })))
# if (pc.calc > 1000) {
# pc.calc <- 1000
# } else if (pc.calc > 500) {
# pc.calc <- 500
# } else if (pc.calc > 100) {
# pc.calc <- ncol(sc10x)-1
# }
sc10x.l <- lapply(sc10x.l,FUN=function(x) { RunPCA(x,features=sc10x.features,npcs=50,verbose=FALSE) })
sc10x.anchors <- FindIntegrationAnchors(object.list=sc10x.l,normalization.method="SCT",anchor.features=sc10x.features,verbose=FALSE,reduction="rpca",dims=1:500)
sc10x.anchors <- FindIntegrationAnchors(object.list=sc10x.l,normalization.method="SCT",anchor.features=sc10x.features,verbose=FALSE,reduction="rpca",dims=1:50)
#sc10x.anchors <- FindIntegrationAnchors(object.list=sc10x.l,normalization.method="SCT",anchor.features=sc10x.features,verbose=FALSE)
sc10x <- IntegrateData(anchorset=sc10x.anchors,normalization.method="SCT",verbose=FALSE)
#sc10x <- FindIntegrationAnchors(object.list=sc10x.l,dims=1:30,scale=FALSE)
#sc10x <- IntegrateData(anchorset=sc10x,dims=1:30)
#sc10x <- FindIntegrationAnchors(object.list=sc10x.l,dims=1:50,scale=FALSE)
#sc10x <- IntegrateData(anchorset=sc10x,dims=1:50)
#gc()
#sc10x <- ScaleData(object=sc10x,do.par=TRUE,num.cores=45,verbose=FALSE,assay="integrated")
#gc()
gc()
sc10x <- ScaleData(object=sc10x,do.par=TRUE,num.cores=45,verbose=FALSE,assay="integrated")
gc()
gc()
sc10x <- SCTransform(sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),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")
sc10x <- SCTransform(sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,return.only.var.genes=FALSE,assay="RNA")
gc()
return(sc10x)
......@@ -593,7 +604,7 @@ scCluster <- function(sc10x,res=0.1,red="pca",dim,print="umap",folder=FALSE,assa
sub <- ""
} else if (print != "0") {
if (!dir.exists(paste0("./analysis/vis/",folder))){
dir.create(paste0("./analysis/vis/",folder))
dir.create(paste0("./analysis/vis/",folder), recursive=TRUE)
}
sub <- paste0(folder,"/")
}
......
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