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

Update cca with new Seurat recommendations

parent a6e6dfb9
#!/bin/bash
#SBATCH --job-name sc10x.raw
#SBATCH -p 256GB,256GBv1,384GB
#SBATCH -N 1
#SBATCH -t 7-0:0:0
#SBATCH -o job_%j.out
#SBATCH -e job_%j.out
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load python/3.6.4-anaconda
source activate umap
module load R/3.5.1-gccmkl
module load hdf5_18/1.8.17
Rscript ../r.scripts/SingleR.R
gc()
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(readxl)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(gridExtra)
library(SingleR)
library(sctransform)
library(autothresholdr)
options(bitmapType="cairo")
option_list=list(
make_option("--p",action="store",default="huBl",type='character',help="Project Name"),
make_option("--s",action="store",default="hu",type='character',help="Species")
)
opt=parse_args(OptionParser(option_list=option_list))
rm(option_list)
setwd("../")
source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")
load("./analysis/sc10x.raw.RData")
load("./genesets/scDWS.250.rda")
load("./genesets/cibersort.rda")
singler = CreateBigSingleRObject(GetAssayData(sc10x,assay="SCT"), annot = NULL, "PdPbPc_deep", min.genes = 0,
singler = CreateBigSingleRObject(GetAssayData(sc10x,assay="SCT"), annot = NULL, "PdPbPc", min.genes = 0,
technology = "10X", species = "Human", citation = "",
ref.list = list(scDWS.250), normalize.gene.length = F, variable.genes = "de",
fine.tune = T, do.signatures = T, clusters = NULL, do.main.types = F,
reduce.file.size = T, numCores = 10)
reduce.file.size = T, numCores = 15)
# cells.epi <-names(singler$singler[[1]]$SingleR.single.main$labels[singler$singler[[1]]$SingleR.single.main$labels[,1]=="Epi",])
# cells.fmst <-names(singler$singler[[1]]$SingleR.single.main$labels[singler$singler[[1]]$SingleR.single.main$labels[,1]=="FMSt",])
......@@ -94,11 +129,11 @@ for (i in c("epi","fmst","st","leu")){
dev.off()
}
singler.leu = CreateSinglerObject(GetAssayData(sc10x.leu,assay="SCT"), annot = NULL, "PdPbPc_deep", min.genes = 0,
singler.leu = CreateSinglerObject(GetAssayData(sc10x.leu,assay="SCT"), annot = NULL, "PdPbPc", min.genes = 0,
technology = "10X", species = "Human", citation = "",
ref.list = list(cibersort), normalize.gene.length = F, variable.genes = "de",
fine.tune = T, do.signatures = T, clusters = NULL, do.main.types = T,
reduce.file.size = T, numCores = SingleR.numCores)
reduce.file.size = T, numCores = 15)
sc10x.leu$leu.lin <- unlist(singler.leu$singler[[1]]$SingleR.single.main$labels[,])
sc10x.leu$leu <- unlist(singler.leu$singler[[1]]$SingleR.single$labels[,])
......@@ -132,4 +167,4 @@ dev.off()
save(list=ls(pattern="^singler"),file='./analysis/singler_objects.RData')
save(list=ls(pattern="^sc10x"),file='./analysis/sc10x.id.rda')
save.image(file="./analysis/sc10x.id.RData")
\ No newline at end of file
save.image(file="./analysis/sc10x.id.RData")
......@@ -17,9 +17,10 @@ library(sctransform)
library(autothresholdr)
options(bitmapType="cairo")
options(future.globals.maxSize= 10000000000)
option_list=list(
make_option("--p",action="store",default="huBl",type='character',help="Project Name"),
make_option("--p",action="store",default="PdPbPc",type='character',help="Project Name"),
make_option("--s",action="store",default="hu",type='character',help="Species")
)
opt=parse_args(OptionParser(option_list=option_list))
......
......@@ -151,7 +151,7 @@ scQC <- function(sc10x,sp="hu"){
for (i in names(sc10x)){
sc10x.temp <- sc10x[[i]]
sc10x.temp[["percent.mito"]] <- PercentageFeatureSet(object=sc10x.temp,pattern=mito.pattern)
sc10x.temp <- subset(sc10x.temp,cell=names(which(is.na(sc10x.temp$percent.mito))),invert=TRUE)
#sc10x.temp <- subset(sc10x.temp,cell=names(which(is.na(sc10x.temp$percent.mito))),invert=TRUE)
sc10x[i] <- sc10x.temp
}
......@@ -170,7 +170,7 @@ scQC <- function(sc10x,sp="hu"){
h[[i]] <- hist(data.frame(sc10x[[j]][[i]])$percent.mito,breaks=1000,plot=FALSE)
cutoff.temp <- mean(c(h[[i]]$mids[which.max(h[[i]]$counts)],h[[i]]$mids[-which.max(h[[i]]$counts)][which.max(h[[i]]$counts[-which.max(h[[i]]$counts)])]))
cells.remove[[j]] <- c(cells.remove[[j]],rownames(sc10x[[j]][["percent.mito"]])[sc10x[[j]][[i]][,1] > cutoff.temp])
sc10x.temp[[j]] <- subset(sc10x[[j]],cells=WhichCells(sc10x[[j]],cells=cells.remove[[j]],invert=TRUE))
sc10x.temp[[j]] <- subset(sc10x[[j]],cells=setdiff(colnames(sc10x[[j]]),cells.remove[[j]]))
thresh.l[[i]] <- scThresh(sc10x.temp,feature=i,sub="lower")
cutoff.l.mito[[j]] <- thresh.l[[i]][[j]]$threshold[thresh.l[[i]][[j]]$method=="Triangle"]
}
......@@ -241,7 +241,7 @@ scQC <- function(sc10x,sp="hu"){
for (i in names(sc10x)){
counts.cell.raw[i] <- ncol(GetAssayData(object=sc10x[[i]],slot="counts"))
counts.gene.raw[i] <- nrow(GetAssayData(object=sc10x[[i]],slot="counts"))
sc10x.sub[[i]] <- subset(sc10x[[i]],cells=WhichCells(sc10x[[i]],cells=cells.remove[[i]],invert=TRUE))
sc10x.sub[[i]] <- subset(sc10x[[i]],cells=setdiff(colnames(sc10x[[i]]),cells.remove[[i]]))
counts.cell.filtered[i] <- ncol(GetAssayData(object=sc10x.sub[[i]],slot="counts"))
counts.gene.filtered[i] <- nrow(GetAssayData(object=sc10x.sub[[i]],slot="counts"))
}
......@@ -455,21 +455,29 @@ scPC <- function(sc10x,pc=50,hpc=0.9,file="pre.stress",print="tsne"){
scCCA <- function(sc10x.l){
for (i in 1:length(sc10x.l)){
sc10x.l[[i]] <- NormalizeData(sc10x.l[[i]],verbose=FALSE)
#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]] <- 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")
gc()
sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
#sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
}
sc10x <- FindIntegrationAnchors(object.list=sc10x.l,dims=1:30,scale=FALSE)
sc10x <- IntegrateData(anchorset=sc10x,dims=1:30)
sc10x.features <- SelectIntegrationFeatures(object.list=sc10x.l,nfeatures=3000)
sc10x.l <- PrepSCTIntegration(object.list=sc10x.l,anchor.features=sc10x.features,verbose=FALSE)
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)
gc()
sc10x <- ScaleData(object=sc10x,do.par=TRUE,num.cores=45,verbose=FALSE,assay="integrated")
gc()
#sc10x <- FindIntegrationAnchors(object.list=sc10x.l,dims=1:30,scale=FALSE)
#sc10x <- IntegrateData(anchorset=sc10x,dims=1:30)
#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"),verbose=FALSE,return.only.var.genes=FALSE,assay="RNA")
gc()
return(sc10x)
}
......
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