SingleR.R 8.83 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
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)
Gervaise Henry's avatar
Gervaise Henry committed
17
library(ggplot2)
18 19 20 21

options(bitmapType="cairo")

option_list=list(
22
  make_option("--p",action="store",default="Biopsy",type='character',help="Project Name"),
23 24
  make_option("--s",action="store",default="hu",type='character',help="Species"),
  make_option("--o",action="store",default="pr",type='character',help="Organ")
25 26 27 28 29 30 31 32 33 34
)
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")

35 36 37
if (opt$s == "hu"){
  hpca.se <- HumanPrimaryCellAtlasData()
  lin.se <- hpca.se
38
  leu.se <- hpca.se
39 40
} else if (opt$s == "mu"){
  igd.se <- ImmGenData()
41 42 43 44 45
  #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
46
}
47

48
if (opt$o == "pr" && opt$s == "hu"){
49
  load("./genesets/Pd.shallow.Rda")
50 51 52 53 54 55
  
  try(
    if (as.numeric(substring(sc10x@version,1,1))<3){
      sc10x <- UpdateSeuratObject(sc10x)
    }
  )
Gervaise Henry's avatar
Gervaise Henry committed
56
  scDWSpr.se <- as.SingleCellExperiment(sc10x,assay="RNA")
57
  scDWSpr.se$Lin[scDWSpr.se$Lin == "Unknown"] <- "Leu"
58 59
  levels(scDWSpr.se$ident)[levels(scDWSpr.se$ident)=="OE2"] <- "Hillock"
  levels(scDWSpr.se$ident)[levels(scDWSpr.se$ident)=="OE1"] <- "Club"
60 61
  rm(sc10x)
}
Gervaise Henry's avatar
Gervaise Henry committed
62 63 64 65

load("./analysis/sc10x.raw.rda")

sc10x.se <- as.SingleCellExperiment(sc10x)
66 67
common <- intersect(rownames(sc10x.se),rownames(lin.se))
lin.se <- lin.se[common,]
Gervaise Henry's avatar
Gervaise Henry committed
68 69 70
sc10x.se <- sc10x.se[common,]
rm(common)

71 72
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)]
73 74
#singler.lin <- SingleR(sc10x.se,ref=lin.se,method="single",labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
#sc10x$lin <- singler.lin$labels
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
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"),
  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"
91
#DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
Gervaise Henry's avatar
Gervaise Henry committed
92 93

Idents(sc10x) <- "lin"
94 95
sc10x.epi <- subset(sc10x,idents="Epi")
sc10x.fmst <- subset(sc10x,idents="FMSt")
Gervaise Henry's avatar
Gervaise Henry committed
96
sc10x.st <- subset(sc10x,idents=setdiff(levels(sc10x),"Epi"))
97 98
sc10x.leu <- subset(sc10x,idents="Leu")

99
if (opt$o == "pr" && opt$s == "hu"){
Gervaise Henry's avatar
Gervaise Henry committed
100
  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi,assay="RNA")
101 102 103 104 105 106
  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)
  
Gervaise Henry's avatar
Gervaise Henry committed
107
  sc10x.se.fmst <- as.SingleCellExperiment(sc10x.fmst,assay="RNA")
108 109 110 111 112 113 114 115 116 117 118 119 120 121
  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"
}
Gervaise Henry's avatar
Gervaise Henry committed
122

Gervaise Henry's avatar
Gervaise Henry committed
123
sc10x.se.leu <- as.SingleCellExperiment(sc10x.leu,assay="RNA")
124 125 126
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,]
Gervaise Henry's avatar
Gervaise Henry committed
127 128 129
sc10x.se.leu <- sc10x.se.leu[common,]
rm(common)

130 131
singler.leu <- SingleR(sc10x.se.leu,ref=leu.se,labels=leu.se$label.main,BPPARAM=MulticoreParam(workers=10))
sc10x.leu$leu <- singler.leu$labels
Gervaise Henry's avatar
Gervaise Henry committed
132
Idents(sc10x.leu) <- "leu"
133
#DimPlot(sc10x.leu,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
Gervaise Henry's avatar
Gervaise Henry committed
134 135

sc10x$pops <- sc10x$lin
136
sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
137 138 139
sc10x.epi$pops <- sc10x.epi$lin
sc10x.fmst$pops <- sc10x.fmst$lin
sc10x.st$pops <- sc10x.st$lin
140
sc10x.leu$pops <- sc10x.leu$leu
Gervaise Henry's avatar
Gervaise Henry committed
141
if (opt$o == "pr" && opt$s == "hu"){
142 143 144 145
  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
Gervaise Henry's avatar
Gervaise Henry committed
146 147
  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
  sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
148 149
  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
}
Gervaise Henry's avatar
Gervaise Henry committed
150 151 152
#DimPlot(sc10x,group.by="pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

sc10x$leu <- "non-Leu"
153 154
sc10x$leu[names(sc10x.leu$leu)] <- sc10x.leu$leu
#DimPlot(sc10x,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
155

156 157
res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))

158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
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")
})
189

Gervaise Henry's avatar
Gervaise Henry committed
190 191 192
if (!dir.exists(paste0("./analysis/vis/singler"))){
  dir.create(paste0("./analysis/vis/singler"))
}
Gervaise Henry's avatar
Gervaise Henry committed
193

Gervaise Henry's avatar
Gervaise Henry committed
194 195 196 197
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()
Gervaise Henry's avatar
Gervaise Henry committed
198 199

Idents(sc10x) <- "pops"
200 201 202 203
Idents(sc10x.epi) <- "pops"
Idents(sc10x.fmst) <- "pops"
Idents(sc10x.st) <- "pops"
Idents(sc10x.leu) <- "leu"
Gervaise Henry's avatar
Gervaise Henry committed
204 205 206 207 208 209

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()

Gervaise Henry's avatar
Gervaise Henry committed
210 211 212 213 214 215
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()
}
216

Gervaise Henry's avatar
Gervaise Henry committed
217
rm(list=ls(pattern="^sc10x.se"))
Gervaise Henry's avatar
Gervaise Henry committed
218

219 220 221 222 223
#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")
224

Gervaise Henry's avatar
Gervaise Henry committed
225 226
save(list=ls(pattern="^singler"),file='./analysis/singler_objects.RData')
save(list=ls(pattern="^sc10x"),file='./analysis/sc10x.id.rda')
227 228
save(list=ls(pattern="^sc10x.epi"),file='./analysis/sc10x.epi.id.rda')
save(list=ls(pattern="^sc10x.st"),file='./analysis/sc10x.st.id.rda')
Gervaise Henry's avatar
Gervaise Henry committed
229
#save.image(file="./analysis/sc10x.id.RData")