Merge branch 'develop' into 'master'

Develop

See merge request !6

bash.scripts/sc10x.RUN.sh

@@ -8,37 +8,9 @@
 #SBATCH --mail-type ALL
 #SBATCH --mail-user gervaise.henry@utsouthwestern.edu
 
-git pull origin Seurat3.0
-rm ../analysis/*.rda
-rm ../analysis/*.RData
-rm -r ../analysis/cor/
-rm -r ../analysis/qc/
-rm -r ../analysis/score_id/
-rm -r ../analysis/shiny/
-rm -r ../analysis/vis/
-
-module load python/3.6.4-anaconda 
-source activate umap
-module load R/3.5.1-gccmkl
+module load gdal
+module load cairo/1.14.8
+module load R/4.0.2-gccmkl 
 module load hdf5_18/1.8.17
-Rscript ../r.scripts/sc-TissueMapper_RUN.R --p "$1" --s "$2"
-
-module unload R/3.5.1-gccmkl
-module load R/3.6.1-gccmkl
-if [[ "$1" == "huPr_PdPgb" ]] || [[ "$1" == "huPr_PdPb" ]] || [[ "$1" == "huBl" ]]
-then
-	Rscript ../r.scripts/SingleR.R --p "$1" --s "$2" --o "$3"
-elif [[ "$1" == "muPrUr" ]]
-then
-	Rscript ../r.scripts/huPr_muPr.R --p "$1" --r "$3"
-fi
-
-module unload R/3.5.1-gccmkl
-module load R/3.6.1-gccmkl
-if [[ "$1" == "huPr_PdPgb" ]]
-then
-	Rscript ../r.scripts/diy_PdPgb.R
-elif [[ "$1" == "muPrUr" ]]
-then
-	Rscript ../r.scripts/diy_muPrUr_Epi.R
-fi
+Rscript ../r.scripts/sc-TissueMapper_RAW.R --p "$1" --s "$2" --m "$3"
+Rscript ../r.scripts/sc-TissueMapper_ID.R --p "$1" --s "$2" --o "$4"

bash.scripts/sc10x.id.sh

 #!/bin/bash
 #SBATCH --job-name sc10x.id
-#SBATCH -p 128GB,256GB,256GBv1,384GB
+#SBATCH -p 256GB,256GBv1,384GB
 #SBATCH -N 1
 #SBATCH -t 7-0:0:0
 #SBATCH -o job_%j.out
@@ -8,9 +8,8 @@
 #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.6.1-gccmkl
+module load gdal
+module load cairo/1.14.8
+module load R/4.0.2-gccmkl 
 module load hdf5_18/1.8.17
-Rscript ../r.scripts/SingleR.R --p "$1" --s "$2" --o "$3"
+Rscript ../r.scripts/sc-TissueMapper_ID.R --p "$1" --s "$2" --o "$4"

bash.scripts/sc10x.raw.sh

@@ -8,9 +8,8 @@
 #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 gdal
+module load cairo/1.14.8
+module load R/4.0.2-gccmkl 
 module load hdf5_18/1.8.17
-Rscript ../r.scripts/sc-TissueMapper_RUN.R --p "$1" --s "$2"
+Rscript ../r.scripts/sc-TissueMapper_RAW.R --p "$1" --s "$2" --m "$3"

bash.scripts/sc_LinkData.sh

 #!/bin/bash
+# Input: project name (eg. huPr_Pd)
 
 mkdir ../analysis/DATA/10x
 
 for i in `cat ../analysis/DATA/${1}-demultiplex.csv`; do
-if [[ ${i} = *Samples* ]]
+if [[ ${i} != *Samples* ]]
 then
-  continue
-else
   sample=`echo ${i} | cut -f1 -d ','`
   mkdir ../analysis/DATA/10x/${sample}
-  ln -s /work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/DATA/"${sample}"/outs/* ../analysis/DATA/10x/${sample}
+  ln -s /work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/DATA/"${sample}"/outs/* ../analysis/DATA/10x/${sample}
 fi
 done

genesets/genes.deg.Stress.csv

-"","p_val","avg_logFC","pct.1","pct.2","p_val_adj"
-"JUN",0,2.78514474389602,0.953,0.609,0
-"FOS",0,2.6180888667142,0.843,0.459,0
-"EGR1",0,2.27744499658095,0.827,0.313,0
-"ATF3",0,2.24777061964067,0.877,0.331,0
-"JUNB",0,1.97088620438955,0.899,0.663,0
-"GADD45B",0,1.92104125732309,0.72,0.423,0
-"IER2",0,1.88242169361821,0.806,0.625,0
-"ZFP36",0,1.86471330757214,0.93,0.655,0
-"DNAJB1",0,1.81462475826555,0.91,0.467,0
-"RHOB",0,1.70021664456101,0.72,0.179,0
-"NR4A1",0,1.63532584305745,0.727,0.208,0
-"UBC",0,1.6303655140661,0.991,0.936,0

genesets/van.den.Brink1.txt

-van.den.Brink1
-Fos
-Hspa1a
-Jun
-Fosb
-Junb
-Egr1
-Hspa1b
-Ubc
-Zfp36
-Hspb1
-Hsp90aa1
-Mt2
-Dnajb1
-Btg2
-Nr4a1
-Cebpd
-Hspa8
-Mt1
-Ier2
-Dnaja1
-Socs3
-Atf3
-Jund
-Cebpb
-Id3
-Ppp1r15a
-Hspe1
-Cxcl1
-Dusp1
-Hsp90ab1
-Nfkbia
-Hsph1

misc/Pd-demultiplex.csv

-Samples,ALL,Pz,Tz,D17,D27,D35
-D17PrPzF_Via,1,1,0,1,0,0
-D17PrTzF_Via,1,0,1,1,0,0
-D27PrTzF_Via,1,0,1,0,1,0
-D27PrPzF_Via,1,1,0,0,1,0
-D35PrPzF_Via,1,1,0,0,0,1
-D35PrTzF_Via,1,0,1,0,0,1

misc/huPr_Pd-demultiplex.csv

+Samples,ALL,Keep,Zone
+D17PrPzF_Via,1,1,Pz
+D17PrTzF_Via,1,1,Tz
+D27PrPzF_Via,1,1,Pz
+D27PrTzF_Via,1,1,Tz
+D35PrPzF_Via,1,1,Pz
+D35PrTzF_Via,1,1,Tz

r.scripts/SingleR.R

-gc()
-library(methods)
-library(optparse)
-library(Seurat)
-library(readr)
-library(readxl)
-library(fBasics)
-library(pastecs)
-library(qusage)
-library(RColorBrewer)
-library(dplyr)
-library(viridis)
-library(gridExtra)
-library(SingleR)
-library(sctransform)
-library(autothresholdr)
-library(ggplot2)
-
-options(bitmapType="cairo")
-
-option_list=list(
-  make_option("--p",action="store",default="Biopsy",type='character',help="Project Name"),
-  make_option("--s",action="store",default="hu",type='character',help="Species"),
-  make_option("--o",action="store",default="pr",type='character',help="Organ")
-)
-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")
-
-if (opt$s == "hu"){
-  hpca.se <- HumanPrimaryCellAtlasData()
-  lin.se <- hpca.se
-  leu.se <- hpca.se
-} else if (opt$s == "mu"){
-  igd.se <- ImmGenData()
-  #lin.se <- igd.se
-  #leu.se <- igd.se
-  lin.se <- subset(igd.se,select=(as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) >= 920648 & as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) <= 2112477))
-  lin.se <- subset(lin.se,select=!(as.numeric(gsub("_","",substr(colnames(lin.se),4,10))) %in% 1398469:1398471))
-  leu.se <- lin.se
-}
-
-if (opt$o == "pr" && opt$s == "hu"){
-  load("./genesets/Pd.shallow.Rda")
-  
-  try(
-    if (as.numeric(substring(sc10x@version,1,1))<3){
-      sc10x <- UpdateSeuratObject(sc10x)
-    }
-  )
-  scDWSpr.se <- as.SingleCellExperiment(sc10x,assay="RNA")
-  scDWSpr.se$Lin[scDWSpr.se$Lin == "Unknown"] <- "Leu"
-  levels(scDWSpr.se$ident)[levels(scDWSpr.se$ident)=="OE2"] <- "Hillock"
-  levels(scDWSpr.se$ident)[levels(scDWSpr.se$ident)=="OE1"] <- "Club"
-  rm(sc10x)
-}
-
-load("./analysis/sc10x.raw.rda")
-
-sc10x.se <- as.SingleCellExperiment(sc10x)
-common <- intersect(rownames(sc10x.se),rownames(lin.se))
-lin.se <- lin.se[common,]
-sc10x.se <- sc10x.se[common,]
-rm(common)
-
-singler.lin <- SingleR(sc10x.se,ref=lin.se,method="cluster",clusters=sc10x.se$integrated_snn_res.5,labels=lin.se$label.fine,BPPARAM=MulticoreParam(workers=10))
-sc10x$lin <- singler.lin$labels[match(sc10x.se$integrated_snn_res.5,singler.lin@rownames)]
-#singler.lin <- SingleR(sc10x.se,ref=lin.se,method="single",labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
-#sc10x$lin <- singler.lin$labels
-sc10x$lin[sc10x$lin %in% unique(lin.se@colData[lin.se@colData[,"label.main"] %in% c(
-  c("DC","B_cell","Neutrophil","T_cells","Monocyte","Macrophage","NK_cell","Neutrophils","CMP","GMP","MEP","Myelocyte","Pre-B_cell_CD34-","Pro-B_cell_CD34+","Pro-Myelocyte","HSC_-G-CSF","HSC_CD34+"),
-  c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")
-  ),c("label.main","label.fine")])$label.fine] <- "Leu"
-sc10x$lin[sc10x$lin %in% unique(lin.se@colData[lin.se@colData[,"label.main"] %in% c(
-  c("Endothelial_cells","Erythroblast","Platelets"),
-  c("Endothelial cells")
-  ),c("label.main","label.fine")])$label.fine] <- "Endo"
-sc10x$lin[sc10x$lin %in% unique(lin.se@colData[lin.se@colData[,"label.main"] %in% c(
-  c("Smooth_muscle_cells","Fibroblasts","Chondrocytes","Osteoblasts","MSC","Tissue_stem_cells","iPS_cells"),
-  c("Stromal cells","Fibroblasts")
-  ),c("label.main","label.fine")])$label.fine] <- "FMSt"
-sc10x$lin[sc10x$lin %in% unique(lin.se@colData[lin.se@colData[,"label.main"] %in% c(
-  c("Epithelial_cells","Keratinocytes","Neuroepithelial_cell"),
-  c("Epithelial cells")
-  ),c("label.main","label.fine")])$label.fine] <- "Epi"
-#DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-Idents(sc10x) <- "lin"
-sc10x.epi <- subset(sc10x,idents="Epi")
-sc10x.fmst <- subset(sc10x,idents="FMSt")
-sc10x.st <- subset(sc10x,idents=setdiff(levels(sc10x),"Epi"))
-sc10x.leu <- subset(sc10x,idents="Leu")
-
-if (opt$o == "pr" && opt$s == "hu"){
-  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi,assay="RNA")
-  scDWSpr.se.epi <- scDWSpr.se[,scDWSpr.se$Lin=="Epi",]
-  common <- intersect(rownames(sc10x.se.epi),rownames(scDWSpr.se.epi))
-  scDWSpr.se.epi <- scDWSpr.se.epi[common,]
-  sc10x.se.epi <- sc10x.se.epi[common,]
-  rm(common)
-  
-  sc10x.se.fmst <- as.SingleCellExperiment(sc10x.fmst,assay="RNA")
-  scDWSpr.se.fmst <- scDWSpr.se[,scDWSpr.se$Merge_Epi.dws_St.go %in% c("Fib","SM"),]
-  common <- intersect(rownames(sc10x.se.fmst),rownames(scDWSpr.se.fmst))
-  scDWSpr.se.fmst <- scDWSpr.se.fmst[common,]
-  sc10x.se.fmst <- sc10x.se.fmst[common,]
-  rm(common)
-  
-  singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se.epi,labels=scDWSpr.se.epi$ident,BPPARAM=MulticoreParam(workers=10))
-  sc10x.epi$scDWSpr <- singler.epi$labels
-  Idents(sc10x.epi) <- "scDWSpr"
-  
-  singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se.fmst,labels=scDWSpr.se.fmst$ident,BPPARAM=MulticoreParam(workers=10))
-  sc10x.fmst$scDWSpr <- singler.fmst$labels
-  Idents(sc10x.fmst) <- "scDWSpr"
-}
-
-sc10x.se.leu <- as.SingleCellExperiment(sc10x.leu,assay="RNA")
-leu.se <- leu.se[,leu.se$label.main %in% c("DC","B_cell","Neutrophil","T_cells","Monocyte","Macrophage","NK_cell","Neutrophils","CMP","GMP","MEP","Myelocyte","Pre-B_cell_CD34-","Pro-B_cell_CD34+","Pro-Myelocyte","Macrophages","Monocytes","B cells","Eosinophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")]
-common <- intersect(rownames(sc10x.se.leu),rownames(leu.se))
-leu.se <- leu.se[common,]
-sc10x.se.leu <- sc10x.se.leu[common,]
-rm(common)
-
-singler.leu <- SingleR(sc10x.se.leu,ref=leu.se,labels=leu.se$label.main,BPPARAM=MulticoreParam(workers=10))
-sc10x.leu$leu <- singler.leu$labels
-Idents(sc10x.leu) <- "leu"
-#DimPlot(sc10x.leu,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-sc10x$pops <- sc10x$lin
-sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
-sc10x.epi$pops <- sc10x.epi$lin
-sc10x.fmst$pops <- sc10x.fmst$lin
-sc10x.st$pops <- sc10x.st$lin
-sc10x.leu$pops <- sc10x.leu$leu
-if (opt$o == "pr" && opt$s == "hu"){
-  sc10x$pops[names(sc10x.epi$scDWSpr)] <- sc10x.epi$scDWSpr
-  sc10x.epi$pops <- sc10x.epi$scDWSpr
-  sc10x$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-  sc10x.fmst$pops <- sc10x.fmst$scDWSpr
-  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-  sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
-  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-}
-#DimPlot(sc10x,group.by="pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-sc10x$leu <- "non-Leu"
-sc10x$leu[names(sc10x.leu$leu)] <- sc10x.leu$leu
-#DimPlot(sc10x,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
-
-try({
-  results <- scPC(sc10x.epi,pc=100,hpc=0.9,file="epi",print="2")
-  sc10x.epi <- results[[1]]
-  pc.use.epi <- results[[2]]
-  rm(results)
-  sc10x.epi <- scCluster(sc10x.epi,res=res,red="pca",dim=pc.use.epi,print="2",folder="epi")
-})
-
-try({
-  results <- scPC(sc10x.fmst,pc=100,hpc=0.9,file="fmst",print="2")
-  sc10x.fmst <- results[[1]]
-  pc.use.fmst <- results[[2]]
-  rm(results)
-  sc10x.fmst <- scCluster(sc10x.fmst,res=res,red="pca",dim=pc.use.fmst,print="2",folder="fmst")
-})
-
-try({
-  results <- scPC(sc10x.st,pc=100,hpc=0.9,file="st",print="2")
-  sc10x.st <- results[[1]]
-  pc.use.st <- results[[2]]
-  rm(results)
-  sc10x.st <- scCluster(sc10x.st,res=res,red="pca",dim=pc.use.st,print="2",folder="st")
-})
-
-try({
-  results <- scPC(sc10x.leu,pc=100,hpc=0.9,file="leu",print="2")
-  sc10x.leu <- results[[1]]
-  pc.use.leu <- results[[2]]
-  rm(results)
-  sc10x.leu <- scCluster(sc10x.leu,res=res,red="pca",dim=pc.use.leu,print="2",folder="leu")
-})
-
-if (!dir.exists(paste0("./analysis/vis/singler"))){
-  dir.create(paste0("./analysis/vis/singler"))
-}
-
-plot <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-postscript(paste0("./analysis/vis/singler/UMAP_all.lin.eps"))
-print(plot)
-dev.off()
-
-Idents(sc10x) <- "pops"
-Idents(sc10x.epi) <- "pops"
-Idents(sc10x.fmst) <- "pops"
-Idents(sc10x.st) <- "pops"
-Idents(sc10x.leu) <- "leu"
-
-plot <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-postscript(paste0("./analysis/vis/singler/UMAP_all.eps"))
-print(plot)
-dev.off()
-
-for (i in c("epi","fmst","st","leu")){
-  plot <- DimPlot(get(paste0("sc10x.",i)),reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-  postscript(paste0("./analysis/vis/singler/UMAP_",i,".eps"))
-  print(plot)
-  dev.off()
-}
-
-rm(list=ls(pattern="^sc10x.se"))
-
-#scShinyOutput(sc10x,anal="id")
-#scShinyOutput(sc10x.epi,anal="id.epi")
-#scShinyOutput(sc10x.fmst,anal="id.fmst")
-#scShinyOutput(sc10x.st,anal="id.st")
-#scShinyOutput(sc10x.leu,anal="id.leu")
-
-save(list=ls(pattern="^singler"),file='./analysis/singler_objects.RData')
-save(list=ls(pattern="^sc10x"),file='./analysis/sc10x.id.rda')
-save(list=ls(pattern="^sc10x.epi"),file='./analysis/sc10x.epi.id.rda')
-save(list=ls(pattern="^sc10x.st"),file='./analysis/sc10x.st.id.rda')
-#save.image(file="./analysis/sc10x.id.RData")

r.scripts/diy_2021FibPaper_CPDB.R

+library(Seurat)
+library(biomaRt)
+
+options(future.globals.maxSize= 1000000000000)
+
+setwd("../")
+
+sc10x <- readRDS("./analysis/huPr_PdPb_id_all.rds")
+sc10x.fmst <- readRDS("./analysis/huPr_PdPb_id_fmst.rds")
+DefaultAssay(sc10x.fmst) <- "SCT"
+Idents(sc10x.fmst) <- "resolution_0.2"
+sc10x.fmst <- RenameIdents(sc10x.fmst,"6"="Prostate SM")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"2"="Prostate SM")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"1"="Pericytes")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"3"="vSM")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"5"="Interstitial Fibroblasts")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"0"="Interstitial Fibroblasts")
+sc10x.fmst <- RenameIdents(sc10x.fmst,"4"="Peri-epithelial Fibroblasts")
+sc10x.fmst$Populations <- Idents(sc10x.fmst)
+Idents(sc10x.fmst) <- "Populations"
+sc10x$pop[names(sc10x.fmst$Populations)] <- as.character(sc10x.fmst$Populations)
+rm(sc10x.fmst)
+
+as_matrix <- function(mat){
+  
+  tmp <- matrix(data=0L, nrow = mat@Dim[1], ncol = mat@Dim[2])
+  
+  row_pos <- mat@i+1
+  col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])+1
+  val <- mat@x
+  
+  for (i in seq_along(val)){
+    tmp[row_pos[i],col_pos[i]] <- val[i]
+  }
+  
+  row.names(tmp) <- mat@Dimnames[[1]]
+  colnames(tmp) <- mat@Dimnames[[2]]
+  return(tmp)
+}
+
+count_norm <- as.data.frame(as_matrix(GetAssayData(sc10x,assay="SCT",slot="data")))
+count_norm$Symbol <- rownames(count_norm)
+ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl",host="https://sep2019.archive.ensembl.org")
+symbol2ID <- getBM(mart = ensembl,attributes = c("external_gene_name","ensembl_gene_id"),filters="external_gene_name",values=count_norm$Symbol,uniqueRows=TRUE)
+colnames(symbol2ID)[colnames(symbol2ID)=="ensembl_gene_id"] <- "Genes"
+count_norm <- merge(x=count_norm,y=symbol2ID,by.x="Symbol",by.y="external_gene_name",all.x=TRUE)
+count_norm$Symbol <- NULL
+count_norm <- count_norm[,c(which(colnames(count_norm)=="Genes"),which(colnames(count_norm)!="Genes"))]
+write.table(count_norm,'./analysis/cellphonedb_count.txt',sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE)
+
+meta_data <- data.frame(Cell=colnames(sc10x),cell_type=as.character(sc10x$pop))
+write.table(meta_data,'./analysis/cellphonedb_meta.txt',sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE)