Commit 9969a1ea authored by Gervaise H. Henry's avatar Gervaise H. Henry 🤠

Copy VAMC15x analysis to VAMC13x

parent e0484b91
gc()
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(readxl)
source("../r.scripts/sc-TissueMapper.R")
#Create folder structure
setwd("../")
if (!dir.exists("./analysis")){
dir.create("./analysis")
}
if (!dir.exists("./analysis/qc")){
dir.create("./analysis/qc")
}
if (!dir.exists("./analysis/qc/cc")){
dir.create("./analysis/qc/cc")
}
if (!dir.exists("./analysis/tSNE")){
dir.create("./analysis/tSNE")
}
if (!dir.exists("./analysis/tSNE/pre.stress")){
dir.create("./analysis/tSNE/pre.stress")
}
if (!dir.exists("./analysis/pca")){
dir.create("./analysis/pca")
}
if (!dir.exists("./analysis/pca/stress")){
dir.create("./analysis/pca/stress")
}
if (!dir.exists("./analysis/violin")){
dir.create("./analysis/violin")
}
if (!dir.exists("./analysis/violin/stress")){
dir.create("./analysis/violin/stress")
}
if (!dir.exists("./analysis/table")){
dir.create("./analysis/table")
}
if (!dir.exists("./analysis/tSNE/post.stress")){
dir.create("./analysis/tSNE/post.stress")
}
if (!dir.exists("./analysis/cor")){
dir.create("./analysis/cor")
}
if (!dir.exists("./analysis/tSNE/lin")){
dir.create("./analysis/tSNE/lin")
}
if (!dir.exists("./analysis/tSNE/epi")){
dir.create("./analysis/tSNE/epi")
}
if (!dir.exists("./analysis/tSNE/st")){
dir.create("./analysis/tSNE/st")
}
if (!dir.exists("./analysis/tSNE/merge")){
dir.create("./analysis/tSNE/merge")
}
if (!dir.exists("./analysis/pca/ne")){
dir.create("./analysis/pca/ne")
}
if (!dir.exists("./analysis/tSNE/ne")){
dir.create("./analysis/tSNE/ne")
}
if (!dir.exists("./analysis/violin/ne")){
dir.create("./analysis/violin/ne")
}
if (!dir.exists("./analysis/tSNE/FINAL")){
dir.create("./analysis/tSNE/FINAL")
}
if (!dir.exists("./analysis/deg")){
dir.create("./analysis/deg")
}
if (!dir.exists("./analysis/cca")){
dir.create("./analysis/cca")
}
if (!dir.exists("./analysis/diy")){
dir.create("./analysis/diy")
}
if (!dir.exists("./analysis/pseudotime")){
dir.create("./analysis/pseudotime")
}
#Retrieve command-line options
option_list=list(
make_option("--p",action="store",default="DPrF",type='character',help="Project Name"),
make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"),
make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"),
make_option("--hg",action="store",default=3000,type='integer',help="Threshold for cells with maximum genes"),
make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"),
make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"),
make_option("--lx",action="store",default=0.2,type='numeric',help="x low threshold for hvg selection"),
make_option("--hx",action="store",default=5,type='numeric',help="x high threshold for hvg selection"),
make_option("--ly",action="store",default=1,type='numeric',help="y low threshold for hvg selection"),
make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"),
make_option("--cca",action="store",default=50,type='integer',help="Number of CCAs to cacluate"),
make_option("--acca",action="store",default=30,type='integer',help="Number of CCAs to align"),
make_option("--pc",action="store",default=50,type='integer',help="Number of PCs to cacluate"),
make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"),
make_option("--st",action="store",default=TRUE,type='logical',help="Remove stressed cells?"),
make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"),
make_option("--cut.stress",action="store",default=0.9,type='numeric',help="Cutoff for stress score"),
make_option("--res.poststress",action="store",default=0.5,type='numeric',help="Resolution to cluster, post-stress"),
make_option("--cut.ne",action="store",default=0.999,type='numeric',help="Cutoff for NE score")
)
opt=parse_args(OptionParser(option_list=option_list))
rm(option_list)
if (opt$lm==0){opt$lm=-Inf}
sc10x <- scLoad("VAMC013PrRdF")
if (opt$cc==TRUE){
results <- scCellCycle(sc10x)
sc10x <- results[[1]]
genes.s <- results[[2]]
genes.g2m <- results[[3]]
rm(results)
} else {
genes.s=""
genes.g2m=""
}
results <- scQC(sc10x,lg=opt$lg,hg=opt$hg,lm=opt$lm,hm=opt$hm)
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
counts.gene.raw <- results[[3]]
counts.cell.filtered <- results[[4]]
counts.gene.filtered <- results[[5]]
rm(results)
gc()
if (opt$cc==TRUE){
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
} else {
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45)
}
gc()
results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=50,hpc=0.85,file="pre.stress",cca=FALSE)
sc10x <- results[[1]]
genes.hvg.prestress <- results[[2]]
pc.use.prestress <- results[[3]]
rm(results)
sc10x <- scCluster(sc10x,pc.use=pc.use.prestress,res.use=opt$res.prestress,folder="pre.stress",red="pca")
if (opt$st==TRUE){
results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,cut=opt$cut.stress)
sc10x <- results[[1]]
counts.cell.filtered.stress <- results[[2]]
sc10x.Stress <- results[[3]]
rm(results)
results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=50,hpc=0.85,file="post.stress",cca=FALSE)
sc10x <- results[[1]]
genes.hvg.poststress <- results[[2]]
pc.use.poststress <- results[[3]]
rm(results)
sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=opt$res.poststress,folder="post.stress",red="pca")
}
gene.set1 <- read_delim("./genesets/genes.deg.Epi.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Epi"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.St.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.BE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "BE"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.LE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "LE"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE1.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "OE_SCGB"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE2.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "OE_KRT13"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.Endo.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Endo"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.SM.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "SM"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.Fib.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Fib"
gene.set <- c(gene.set,gene.set1)
genes.leu <- read_excel("genesets/40425_2017_215_MOESM1_ESM.xlsx",sheet="S3. Candidate markers")
leu <- as.list(unique(genes.leu[,2]))$Cell
leu.l <- leu[-c(1,3,4,7:9,14,15,17:18,20:21,23:30)]
genes.leu <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu),]
genes.leu.l <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu.l),]
gene.set1 <- split(genes.leu.l[,1], genes.leu.l[,2])
gene.set1 <- lapply(gene.set1,unname)
gene.set1 <- lapply(gene.set1,unlist)
gene.set <- c(gene.set,gene.set1)
rm(gene.set1)
gc()
min.all <- min(table(sc10x@meta.data[,paste0("res",opt$res.poststress)]))
results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=min.all,nm="Pop",folder="lin")
sc10x <- results[[1]]
results.cor.Lin <- results[[2]]
results.clust.Lin.id <- results[[3]]
rm(results)
rm(gene.set)
sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne)
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