sc-TissueMapper_RAW.R 6.18 KB
Newer Older
1
2
3
4
5
gc()
library(optparse)
library(Seurat)
library(readr)
library(dplyr)
Gervaise Henry's avatar
Gervaise Henry committed
6
7
library(autothresholdr)
library(RColorBrewer)
8
9
library(viridis)
library(gridExtra)
10
library(ggplot2)
11

Gervaise Henry's avatar
Gervaise Henry committed
12
options(future.globals.maxSize= 1000000000000)
13

14
option_list=list(
15
  make_option("--p",action="store",default="huPr_PdPb",type='character',help="Project Name"),
Gervaise Henry's avatar
Gervaise Henry committed
16
  make_option("--s",action="store",default="hu",type='character',help="Species"),
Gervaise Henry's avatar
Gervaise Henry committed
17
  make_option("--m",action="store",default="Patient",type='character',help="Group to merge on")
18
19
20
21
22
)
opt=parse_args(OptionParser(option_list=option_list))
rm(option_list)


23
24
25
26
27
setwd("../")

source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")

28
project.name=opt$p
Gervaise Henry's avatar
Gervaise Henry committed
29

30
31
scFolders()

Gervaise Henry's avatar
Gervaise Henry committed
32
results <- scLoad(p=project.name,cellranger=4,aggr=FALSE,ncell=0,nfeat=0,seurat=TRUE)
33
34
35
sc10x <- results[[1]]
sc10x.groups <- results[[2]]
rm(results)
36

Gervaise Henry's avatar
Gervaise Henry committed
37
38
sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$ALL==1,1]))]

Gervaise Henry's avatar
Gervaise Henry committed
39
results <- scQC(sc10x,sp=opt$s,feature=c("nFeature_RNA"))
40
41
42
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
counts.gene.raw <- results[[3]]
Gervaise Henry's avatar
Gervaise Henry committed
43
44
counts.cell.filtered.1umi <- results[[4]]
counts.gene.filtered.1umi <- results[[5]]
45
46
rm(results)

Gervaise Henry's avatar
Gervaise Henry committed
47
48
49
50
51
results <- scQC(sc10x,sp=opt$s,feature="percent.mito")
sc10x <- results[[1]]
counts.cell.filtered.2mito <- results[[4]]
counts.gene.filtered.2mito <- results[[5]]
rm(results)
52

Gervaise Henry's avatar
Gervaise Henry committed
53
54
55
56
57
#results <- scQC(sc10x,sp=opt$s,feature="percent.ribo")
#sc10x <- results[[1]]
#counts.cell.filtered.3ribo <- results[[4]]
#counts.gene.filtered.3ribo <- results[[5]]
#rm(results)
Gervaise Henry's avatar
Gervaise Henry committed
58

Gervaise Henry's avatar
Gervaise Henry committed
59
60
61
62
63
#results <- scQC(sc10x,sp=opt$s,feature="nCount_RNA")
#sc10x <- results[[1]]
#counts.cell.filtered.3gene <- results[[4]]
#counts.gene.filtered.3gene <- results[[5]]
#rm(results)
64

Gervaise Henry's avatar
Gervaise Henry committed
65
66
counts.cell.filtered <- counts.cell.filtered.2mito
counts.gene.filtered <- counts.gene.filtered.2mito
67

Gervaise Henry's avatar
Gervaise Henry committed
68
sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))]
69

Gervaise Henry's avatar
Gervaise Henry committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
# pc.use.prestress <- list()
# for (i in names(sc10x)){
#   sc10x.temp <- sc10x[[i]]
#   sc10x.temp <- SCTransform(sc10x.temp,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
#   if (ncol(sc10x.temp) > 100) {
#     pc.calc <- 100
#   } else if (ncol(sc10x.temp) <= 100) {
#     pc.calc <- ncol(sc10x.temp)-1
#   }
#   results <- scPC(sc10x.temp,pc=pc.calc,hpc=0.9,file=paste0(i,".pre.stress"),print="umap",assay="SCT")
#   sc10x.temp <- results[[1]]
#   pc.use.prestress.temp <- results[[2]]
#   rm(results)
#   sc10x.temp <- scCluster(sc10x.temp,res=0.5,red="pca",dim=pc.use.prestress.temp,print="umap",folder=paste0("preStress/",i),assay="SCT")
#   sc10x[i] <- sc10x.temp
#   pc.use.prestress[i] <- pc.use.prestress.temp
#   rm(sc10x.temp)
#   rm(pc.calc)
#   rm(pc.use.prestress.temp)
# }
# rm(i)
91

92
merges <- NULL
Gervaise Henry's avatar
Gervaise Henry committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# for (i in names(counts.cell.destress[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))])){
#   if (counts.cell.destress[[i]]<750){
#     merges <- c(merges,i)
#   }
# }
# rm(i)
# if (length(merges)>1){
#   sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
#   sc10x[merges] <- NULL
#   sc10x$merge <- sc10x.merge
#   rm(sc10x.merge)
# } else if (length(merges)==1){
#   merges <- c(merges,names(counts.cell.destress[counts.cell.destress == min(sapply(counts.cell.destress[setdiff(as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1])),merges)],min))]))
#   sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
#   sc10x[merges] <- NULL
#   sc10x$merge <- sc10x.merge
# }
# rm(merges)
merges <- NULL
for (i in unique(sc10x.groups$Patient[sc10x.groups$Keep==1])){
  for (j in sc10x.groups$Samples[sc10x.groups$Keep==1]){
Gervaise Henry's avatar
Gervaise Henry committed
114
    if (sc10x.groups[sc10x.groups$Samples==j,][opt$merge]==i){
Gervaise Henry's avatar
Gervaise Henry committed
115
116
      merges <- c(merges,j)
    }
117
  }
Gervaise Henry's avatar
Gervaise Henry committed
118
  if (length(merges)>1){
Gervaise Henry's avatar
Gervaise Henry committed
119
  sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
Gervaise Henry's avatar
Gervaise Henry committed
120
121
122
  } else {
    sc10x.merge <- sc10x[[merges[1]]]
  }
Gervaise Henry's avatar
Gervaise Henry committed
123
  sc10x[merges] <- NULL
Gervaise Henry's avatar
Gervaise Henry committed
124
125
  sc10x[i] <- sc10x.merge
  merges <- NULL
Gervaise Henry's avatar
Gervaise Henry committed
126
127
  rm(sc10x.merge)
}
128
rm(merges)
Gervaise Henry's avatar
Gervaise Henry committed
129
130
rm(i)
rm(j)
131

132
gc()
133
if (length(sc10x)>1){
Gervaise Henry's avatar
Gervaise Henry committed
134
  sc10x <- scAlign(sc10x)
135
136
137
  sc10x@project.name <- project.name
  sc10x$ALL <- "ALL"
}
138
gc()
139

140
141
142
143
144
145
146
if (ncol(sc10x) > 1000) {
  pc.calc <- 1000
} else if (ncol(sc10x) > 500) {
  pc.calc <- 500
} else if (ncol(sc10x) > 100) {
  pc.calc <- ncol(sc10x)-1
}
Gervaise Henry's avatar
Gervaise Henry committed
147
results <- scPC(sc10x,pc=pc.calc,hpc=0.9,file="ALL",print="umap")
148
sc10x <- results[[1]]
149
pc.use.poststress <- results[[2]]
150
rm(results)
151
rm(pc.calc)
152
153

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

Gervaise Henry's avatar
Gervaise Henry committed
155
sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="umap",folder="ALL")
Gervaise Henry's avatar
Gervaise Henry committed
156

Gervaise Henry's avatar
Gervaise Henry committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
if (opt$s == "hu"){
  genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
  genes.stress <- genes.stress[1]
  colnames(genes.stress) <- "scDWS.Stress"
  anchor <- c("EGR1","FOS","JUN")
} else if (opt$s == "mu"){
  genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/van.den.Brink1.txt",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
  genes.stress <- genes.stress[1]
  colnames(genes.stress) <- "van.den.Brink1.Stress"
  anchor <- c("Egr1","Fos","Jun")
}

results <- scScore(list(ALL=sc10x),score="Stress",geneset=as.list(genes.stress),cut.pt="renyi",anchor=anchor)
sc10x.preStress <- results[[1]]
sc10x <- results[[2]]
#sc10x <- results[[1]]$ALL
rm(results)
#counts.cell.destress <- lapply(sc10x,ncol)
counts.cell.destress <- ncol(sc10x)
rm(anchor)

Gervaise Henry's avatar
Gervaise Henry committed
178
DefaultAssay(object=sc10x) <- "SCT"
179

180
sc10x@meta.data <- sc10x@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
181

182
saveRDS(sc10x,paste0("./analysis/",project.name,"_raw.rds"))
183

Gervaise Henry's avatar
Gervaise Henry committed
184
185
186
library(sceasy)
library(reticulate)
use_condaenv('sceasy')
187
188
convertFormat(sc10x,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.h5ad"),assay="SCT",main_layer="scale.data")
saveRDS(sc10x,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",project.name,"_raw.rds"))
Gervaise Henry's avatar
Gervaise Henry committed
189
190
191

rm(list=ls(pattern="sc10x"))
save.image(file="./analysis/sc10x.raw.RData")