huPr_muPr.R 7.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
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(
22 23
  make_option("--p",action="store",default="muPrUr",type='character',help="Project Name"),
  make_option("--r",action="store",default="epi",type='character',help="Reference Subset")
24 25 26 27 28 29 30 31 32 33
)
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")

34 35 36 37 38
igd.se <- ImmGenData()
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
rm(igd.se)
39

Gervaise Henry's avatar
Gervaise Henry committed
40
convert <- read.delim("./genesets/Ensemble.mus-hum.txt")
41 42 43 44

if (opt$r == "epi"){
  load("./genesets/sc10x.epi.id.rda")
  sc10x <- sc10x.epi
45
  rm(sc10x.epi)
46 47 48 49 50 51 52 53 54 55
  Idents(sc10x) <- "pops"
  sc10x <- subset(sc10x,idents=c("NE"),invert=TRUE)
  sc10x <- RenameIdents(sc10x,"Hillock"="Ur")
  sc10x <- RenameIdents(sc10x,"Club"="Ur")
  sc10x$pops <- Idents(sc10x)
  sc10x$pops <- factor(sc10x$pops,levels=c("BE","LE","Ur"))
  Idents(sc10x) <- "pops"
} else if (opt$r == "fmst"){
  load("./genesets/sc10x.fmst.id.rda")
  sc10x <- sc10x.st
56
  rm(sc10x.st)
57
  Idents(sc10x) <- "pops"
58 59
}

60 61 62 63 64 65 66 67 68 69 70 71
try(
  if (as.numeric(substring(sc10x@version,1,1))<3){
    sc10x <- UpdateSeuratObject(sc10x)
  }
)
Idents(sc10x) <- "pops"

scDWSpr.se <- as.SingleCellExperiment(sc10x,assay="RNA")
rm(sc10x)
scDWSpr.se <- scDWSpr.se[rownames(scDWSpr.se) %in% convert$Human.gene.name,]
rownames(scDWSpr.se) <- convert$Gene.name[match(rownames(scDWSpr.se),convert$Human.gene.name)]

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

74
sc10x.se <- as.SingleCellExperiment(sc10x,assay="RNA")
75 76 77 78 79
common <- intersect(rownames(sc10x.se),rownames(lin.se))
lin.se <- lin.se[common,]
sc10x.se <- sc10x.se[common,]
rm(common)

80 81 82 83 84 85 86 87 88
singler.lin <- SingleR(sc10x.se,ref=lin.se,method="cluster",clusters=sc10x.se$integrated_snn_res.0.5,labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
sc10x$lin <- singler.lin$labels[match(sc10x.se$integrated_snn_res.0.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% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")] <- "Leu"
sc10x$lin[sc10x$lin %in% c("Endothelial cells")] <- "Endo"
sc10x$lin[sc10x$lin %in% c("Stromal cells","Fibroblasts")] <- "FMSt"
sc10x$lin[sc10x$lin %in% c("Epithelial cells")] <- "Epi"
#DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
89 90

Idents(sc10x) <- "lin"
91 92
sc10x.epi <- subset(sc10x,idents=c("Epi"))
sc10x.fmst <- subset(sc10x,idents=c("FMSt"))
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110

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")
})

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
if (opt$r == "epi"){
  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi,assay="RNA")
  common <- intersect(rownames(sc10x.se.epi),rownames(scDWSpr.se))
  scDWSpr.se <- scDWSpr.se[common,]
  sc10x.se.epi <- sc10x.se.epi[common,]
  rm(common)
  
  singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se,method="cluster",clusters=sc10x.se.epi$integrated_snn_res.0.1,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
  sc10x.epi$hu_pops <- singler.epi$labels[match(sc10x.se.epi$integrated_snn_res.0.1,singler.epi@rownames)]
  #singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
  #sc10x.epi$hu_pops <- singler.epi$labels
  Idents(sc10x.epi) <- "hu_pops"
  sc10x.epi$hu_pops <- factor(sc10x.epi$hu_pops)
  #DimPlot(sc10x.epi,group.by="hu_pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
  sc10x.fmst$hu_pops <- sc10x.fmst$lin
} else if (opt$r == "fmst"){
  sc10x.se.fmst <- as.SingleCellExperiment(sc10x.fmst,assay="RNA")
  common <- intersect(rownames(sc10x.se.fmst),rownames(scDWSpr.se))
  scDWSpr.fmse <- scDWSpr.se[common,]
  sc10x.se.fmst <- sc10x.se.fmst[common,]
  rm(common)
  
  singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se,method="cluster",clusters=sc10x.se.fmst$integrated_snn_res.0.1,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
  sc10x.fmst$hu_pops <- singler.fmst$labels[match(sc10x.se.fmst$integrated_snn_res.0.1,singler.fmst@rownames)]
  #singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
  #sc10x.fmst$hu_pops <- singler.fmst$labels
  Idents(sc10x.fmst) <- "hu_pops"
  sc10x.fmst$hu_pops <- factor(sc10x.fmst$hu_pops)
  #DimPlot(sc10x.fmst,group.by="hu_pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
  sc10x.epi$hu_pops <- sc10x.epi$lin
}
142 143

Idents(sc10x) <- "lin"
144 145
Idents(sc10x.epi) <- "hu_pops"
Idents(sc10x.fmst) <- "hu_pops"
146 147 148 149 150 151 152 153 154 155

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

156
for (i in c("epi","fmst")){
157 158 159 160 161 162
  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()
}

163 164 165 166 167 168
labs <- singler.lin$labels
labs <- replace(labs,labs=="Epithelial cells","Epi")
labs <- replace(labs,labs %in% c("Stromal cells","Fibroblasts"),"FMSt")
labs <- replace(labs,labs == "Endothelial cells","Endo")
labs <- replace(labs,labs %in% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia"),"Leu")
plot <- plotScoreHeatmap(singler.lin,show.labels=FALSE,annotation_col=data.frame(ID=labs,row.names=rownames(singler.lin)))
169
postscript(paste0("./analysis/vis/singler/ScoreHeatmap_lin.norm.eps"))
170 171
print(plot)
dev.off()
172
rm(labs)
173

174 175 176 177 178 179 180 181 182 183 184
if (opt$r == "epi"){
  plot <- plotScoreHeatmap(singler.epi,show.labels=FALSE,show.pruned=FALSE,annotation_col=data.frame(ID=singler.epi$labels,row.names=rownames(singler.epi)))
  postscript(paste0("./analysis/vis/singler/ScoreHeatmap_epi.norm.eps"))
  print(plot)
  dev.off()
} else if (opt$r == "fmst"){
plot <- plotScoreHeatmap(singler.fmst,show.labels=FALSE,show.pruned=FALSE,annotation_col=data.frame(ID=singler.fmst$labels,row.names=rownames(singler.fmst)))
  postscript(paste0("./analysis/vis/singler/ScoreHeatmap_fmst.norm.eps"))
  print(plot)
  dev.off()
}
185 186 187

rm(list=ls(pattern="^sc10x.se"))

188 189 190 191 192
#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")
193 194 195 196

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