1 options(stringsAsFactors = F )
2 library(plyr)
3 library(network)
4 library(tidygraph)
5 library(igraph)
6 library(ggraph)
7 library(scales)
8 library(STRINGdb)
9 library(Seurat)
10 library(progress)
11 library(lattice)
12 #library(tidyverse)
13 library(ggplot2)
14 library(Matrix)
15 #library(Rmisc)
16 library(ggforce)
17 library(rjson)
18 library(cowplot)
19 library(RColorBrewer)
20 library(grid)
21 library(sp)
22
23 #library(readbitmap)
24 library(ggExtra)
25 library(reshape2)
26 library(gridExtra)
27 library(sctransform)
28 library(pheatmap)
29 library(Hmisc)#加载包
30 library(magick)
31 library(imager)
32 library(seewave)
33 library(MASS)
34 library(scales)
35
36 #load CV/PV genes
37 dataPath <- "C:/Gu_lab/Liver/Data/"
38 savePath <- "C:/Gu_lab/Liver/results/"
39 CV_gl <- read.csv(file = paste0(dataPath,"CV markers.csv"))
40 PV_gl <- read.csv(file = paste0(dataPath,"PV markers.csv"))
41 CV_gl <- CV_gl$FeatureName
42 PV_gl <- PV_gl$FeatureName
43 paper_gl <- c("Glul","Cyp2e1","Ass1","Asl","Alb","Cyp2f2","Cyp8b1","Cdh1")
44 #S1_Positive
45
46 library(scCancer)
47 library(cowplot)
48 library(diptest)
49
50 #
51 stat.results <- runScStatistics(dataPath = 'C:/Gu_lab/Liver/Data/S1_positive/outs/',
52 savePath = 'C:/Gu_lab/Liver/Data/S1_positive/outs/',
53 sampleName = 'S1_Positive',
54 genReport = T,
55 species = "mouse",
56 bool.runSoupx = F)
57
58 anno.results <- runScAnnotation(
59 dataPath = paste0( 'C:/Gu_lab/Liver/Data/S1_positive/outs/'),
60 statPath = paste0( 'C:/Gu_lab/Liver/Data/S1_positive/outs/'),
61 savePath = paste0( 'C:/Gu_lab/Liver/Data/S1_positive/outs/'),
62 sampleName = 'S1_Positive',
63 species = "mouse",
64 anno.filter = c("mitochondrial", "ribosome", "dissociation"),
65 nCell.min = 3, bgPercent.max = 1.0,
66 vars.to.regress = c("nUMI", "mito.percent", "ribo.percent"),
67 pc.use = 30,
68 clusterStashName = "default",
69 bool.runDiffExpr = T,
70 bool.runMalignancy = T,
71 genReport = T
72 )
73
74 stat.results <- runScStatistics(dataPath = 'C:/Gu_lab/Liver/Data/S1_Negative/',
75 savePath = 'C:/Gu_lab/Liver/Data/S1_Negative/',
76 sampleName = 'S1_Negative',
77 genReport = T,
78 species = "mouse",
79 bool.runSoupx = F)
80
81 anno.results <- runScAnnotation(
82 dataPath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/'),
83 statPath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/'),
84 savePath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/'),
85 sampleName = 'S1_Negative',
86 species = "mouse",
87 anno.filter = c("mitochondrial", "ribosome", "dissociation"),
88 nCell.min = 3, bgPercent.max = 1.0,
89 vars.to.regress = c("nUMI", "mito.percent", "ribo.percent"),
90 pc.use = 30,
91 clusterStashName = "default",
92 bool.runDiffExpr = T,
93 bool.runMalignancy = T,
94 genReport = T
95 )
96
97 #P1
98 S1P_dataPath <- "C:/Gu_lab/Liver/Data/S1_positive/outs/Nofilter_0.75/"
99 S1P_savePath <- "C:/Gu_lab/Liver/results/S1P/MT_0.7/"
100 S1N_dataPath <- "C:/Gu_lab/Liver/Data/S1_Negative/"
101 S1N_savePath <- "C:/Gu_lab/Liver/results/S1N/"
102 merge_savePath <- "C:/Gu_lab/Liver/results/merge_MNN//"
103 S1_Positive <- readRDS(file = paste0(S1P_dataPath,"expr.RDS"))
104 S1_Positive <- RenameCells(S1_Positive,add.cell.id = "P",for.merge = T )
105 S1_Negative <- readRDS(file = paste0(S1N_dataPath,"expr.RDS"))
106 S1_Negative <- RenameCells(S1_Negative,add.cell.id = "N",for.merge = T )
107
108
109 #MNN
110
111 expr_all <- FindIntegrationAnchors(c(S1_Positive,S1_Negative))
112 expr_all <- IntegrateData(anchorset = expr_all, dims = 1:30)
113 expr_all <- FindVariableFeatures(expr_all)
114 expr_all <- ScaleData(expr_all,features = VariableFeatures(expr_all))
115 #linear reg
116 gene <- intersect(rownames(S1_Negative),rownames(S1_Positive))
117 expr_all_matrix <- cbind(S1_Negative@assays$RNA@counts[gene,],S1_Positive@assays$RNA@counts[gene,])
118 expr_all <- CreateSeuratObject(expr_all_matrix)
119 expr_all <- NormalizeData(expr_all)
120 expr_all <- FindVariableFeatures(expr_all)
121
122 source <- colnames(expr_all)
123 source <- filter_str(source)
124
125 expr_all[["sample.ident"]] <- source
126 expr_all <- ScaleData(object = expr_all,
127 features = union(VariableFeatures(expr_all),union(PV_gl,CV_gl)),
128 vars.to.regress = c("sample.ident"),
129 verbose = F)
130
131 #Harmony
132 expr_all_matrix <- cbind(S1_Negative@assays$RNA@counts[gene,],S1_Positive@assays$RNA@counts[gene,])
133 expr_all <- CreateSeuratObject(expr_all_matrix)
134 expr_all <- RunHarmony(expr_all, "dataset")
135
136
137 #############
138 expr_all <- RunPCA(expr_all)
139 expr_all <- FindNeighbors(expr_all)
140 expr_all <- FindClusters(expr_all,resolution = 0.6)
141 #expr_all_0.6 <- FindClusters(expr_all,resolution = 0.6)
142 ElbowPlot(expr_all)
143 expr_all <- RunUMAP(expr_all,dims = 1:20)
144 DimPlot(expr_all)
145
146 saveRDS(expr_all,file = paste0(merge_savePath,"expr_all_MNN_0.6.RDS"))
147 saveRDS(expr_all,file = paste0(merge_savePath,"expr_filter_MNN_0.6.RDS"))
148
149
150 #delete unused cell and rerun QC
151 expr_all <- SubsetData(expr_all, ident.remove = c("10","9"))
152 expr_all <- NormalizeData(expr_all)
153 expr_all <- FindVariableFeatures(expr_all)
154
155 source <- colnames(expr_all)
156 source <- filter_str(source)
157
158 expr_all[["sample.ident"]] <- source
159 expr_all <- ScaleData(object = expr_all,
160 features = union(VariableFeatures(expr_all),union(PV_gl,CV_gl)),
161 vars.to.regress = c("sample.ident"),
162 verbose = F)
163
164 expr_all <- RunPCA(expr_all)
165 expr_all <- FindNeighbors(expr_all)
166 expr_all <- FindClusters(expr_all,resolution = 0.6)
167 expr_all <- RunUMAP(expr_all,dims = 1:20)
168
169 #############################
170 expr_all_manifast <- data.frame(barcode = colnames(expr_all),
171 source = source,
172 idents_res0.8 = Idents(expr_all),
173 UMAP_x = expr_all@reductions$umap@cell.embeddings[,1],
174 UMAP_y = expr_all@reductions$umap@cell.embeddings[,2],
175 PCA_x = expr_all@reductions$pca@cell.embeddings[,1],
176 PCA_y = expr_all@reductions$pca@cell.embeddings[,2])
177
178 #Plot MT vs RiboP & MT vs Dissoc.Genes
179 #S1_Positive_raw <- Read10X(paste0(S1P_dataPath,"/filtered_feature_bc_matrix"))
180 #S1_Positive_raw[["percent.mt"]] <- PercentageFeatureSet(S1_Positive_raw, pattern = "^Mt-")
181
182
183 #nofilter
184
185
186 anno.results <- runScAnnotation(
187 dataPath = paste0('C:/Gu_lab/Liver/Data/S1_positive/outs/'),
188 statPath = paste0('C:/Gu_lab/Liver/Data/S1_positive/outs/'),
189 savePath = paste0('C:/Gu_lab/Liver/Data/S1_positive/outs/Nofilter_0.75/'),
190 sampleName = 'S1_Positive',
191 species = "mouse",
192 anno.filter = c("mitochondrial", "ribosome", "dissociation"),
193 nCell.min = 3, bgPercent.max = 1.0,
194 vars.to.regress = c("nUMI", "mito.percent", "ribo.percent"),
195 pc.use = 30,
196 clusterStashName = "default",
197 bool.runDiffExpr = T,
198 bool.runMalignancy = T,
199 genReport = T
200 )
201
202
203
204 anno.results <- runScAnnotation(
205 dataPath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/'),
206 statPath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/'),
207 savePath = paste0( 'C:/Gu_lab/Liver/Data/S1_Negative/Nofilter_0.75/'),
208 sampleName = 'S1_Negative',
209 species = "mouse",
210 anno.filter = c("mitochondrial", "ribosome", "dissociation"),
211 nCell.min = 3, bgPercent.max = 1.0,
212 vars.to.regress = c("nUMI", "mito.percent", "ribo.percent"),
213 pc.use = 30,
214 clusterStashName = "default",
215 bool.runDiffExpr = T,
216 bool.runMalignancy = T,
217 genReport = T
218 )
219
220
221 #########################################################
222 expr_all <- readRDS(file = paste0(S1P_dataPath,"/Nofilter_2/expr.RDS"))
223 S1_Positive <- RunUMAP(S1_Positive,dim = 1:20)
224 S1P_manifast <- data.frame(barcode = colnames(S1_Positive),
225 idents_res0.8 = Idents(S1_Positive),
226 UMAP_x = S1_Positive@reductions$umap@cell.embeddings[,1],
227 UMAP_y = S1_Positive@reductions$umap@cell.embeddings[,2],
228 PCA_x = S1_Positive@reductions$pca@cell.embeddings[,1],
229 PCA_y = S1_Positive@reductions$pca@cell.embeddings[,2])
230
231
232 S1_Negative <- RunUMAP(S1_Negative,dim = 1:20)
233 S1N_manifast <- data.frame(barcode = colnames(S1_Negative),
234 idents_res0.8 = Idents(S1_Negative),
235 UMAP_x = S1_Negative@reductions$umap@cell.embeddings[,1],
236 UMAP_y = S1_Negative@reductions$umap@cell.embeddings[,2],
237 PCA_x = S1_Negative@reductions$pca@cell.embeddings[,1],
238 PCA_y = S1_Negative@reductions$pca@cell.embeddings[,2])
239
240
241
242 expr_show <- data.frame(barcode = colnames(expr_all),
243 idents_res0.8 = Idents(expr_all),
244 UMAP_x = expr_all@reductions$umap@cell.embeddings[,1],
245 UMAP_y = expr_all@reductions$umap@cell.embeddings[,2],
246 PCA_x = expr_all@reductions$pca@cell.embeddings[,1],
247 PCA_y = expr_all@reductions$pca@cell.embeddings[,2],
248 MT.percent = expr_all@meta.data$mito.percent*100,
249 Ribo.percent = expr_all@meta.data$ribo.percent*100,
250 Diss.percent = expr_all@meta.data$diss.percent*100,
251 nCount_RNA = expr_all@meta.data$nCount_RNA,
252 nFeature_RNA = expr_all@meta.data$nFeature_RNA)
253
254 Idents <- array(rep("Delete",length(rownames(expr_show))))
255 rownames(Idents) <- expr_show$barcode
256 Idents[S1P_manifast$barcode] <- as.character(S1P_manifast$idents_res0.8)
257 Idents[S1N_manifast$barcode] <- as.character(S1N_manifast$idents_res0.8)
258
259
260 #idents <- data.frame(barcode = expr_show$barcode, idents = Idents)
261 #idents <- merge(idents,S1P_manifast,by.x = "barcode",by.y = "barcode")
262
263 expr_show$idents_res0.8 <- Idents
264
265 expr_show_2 <- data.frame(barcode = colnames(expr_all),
266 idents_res0.8 = Idents(expr_all),
267 UMAP_x = expr_all@reductions$umap@cell.embeddings[,1],
268 UMAP_y = expr_all@reductions$umap@cell.embeddings[,2],
269 PCA_x = expr_all@reductions$pca@cell.embeddings[,1],
270 PCA_y = expr_all@reductions$pca@cell.embeddings[,2],
271 MT.percent = expr_all@meta.data$mito.percent,
272 Ribo.percent = expr_all@meta.data$ribo.percent,
273 Diss.percent = expr_all@meta.data$diss.percent)
274
275 cols = as.array(getDefaultColors(length(table(expr_show$idents_res0.8))))
276 rownames(cols) = names(table(expr_show$idents_res0.8))
277 xpand = 0
278 ypand = 1
279
280 ggplot(expr_show,aes(x =MT.percent,y =Ribo.percent, fill = idents_res0.8)) +
281 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
282 #scale_color_gradientn(colours=pal(5))+
283 #scale_color_brewer("Label") +
284 scale_fill_manual(values = cols, name = "idents")+
285 labs( x = "MT.percent", y = "Ribo.percent") +
286 geom_vline(xintercept = 50, colour = "red", linetype = "dashed") +
287 geom_hline(yintercept = 19.8, colour = "red", linetype = "dashed") +
288 ggplot_config(base.size = 6)
289
290 ggplot(expr_show,aes(x =MT.percent,y = Diss.percent, fill = idents_res0.8)) +
291 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
292 #scale_color_gradientn(colours=pal(5))+
293 #scale_color_brewer("Label") +
294 scale_fill_manual(values = cols, name = "idents")+
295 labs( x = "MT.percent", y = "Diss.percent") +
296 geom_vline(xintercept = 50, colour = "red", linetype = "dashed") +
297 geom_hline(yintercept = 3.7, colour = "red", linetype = "dashed") +
298 ggplot_config(base.size = 6)
299
300
301 ggplot(expr_show,aes(x =MT.percent,y = nCount_RNA, fill = idents_res0.8)) +
302 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
303 #scale_color_gradientn(colours=pal(5))+
304 #scale_color_brewer("Label") +
305 scale_fill_manual(values = cols, name = "idents")+
306 labs( x = "MT.percent", y = "nCount_RNA") +
307 geom_vline(xintercept = 50, colour = "red", linetype = "dashed") +
308 ggplot_config(base.size = 6)
309
310 ggplot(expr_show,aes(x =MT.percent,y = nFeature_RNA, fill = idents_res0.8)) +
311 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
312 #scale_color_gradientn(colours=pal(5))+
313 #scale_color_brewer("Label") +
314 scale_fill_manual(values = cols, name = "idents")+
315 labs( x = "MT.percent", y = "nFeature_RNA") +
316 geom_vline(xintercept = 50, colour = "red", linetype = "dashed") +
317 ggplot_config(base.size = 6)
318
319 #############################################################################
320
321 saveRDS(expr_all,file = paste0(merge_savePath,"expr_0.6.RDS"))
322
323 cols = as.array(getDefaultColors(length(table(Idents(expr_all)))))
324 rownames(cols) = names(table(Idents(expr_all)))
325 xpand = 0
326 ypand = 1
327 #plot(expr_all_manifast$UMAP_x,expr_all_manifast$UMAP_y, pch=19, col= cols[expr_all_manifast$idents_res0.8],xlab = "UMAP_1",ylab = "UMAP_2")
328 p <- ggplot(expr_all_manifast,aes(x =UMAP_x,y =UMAP_y ,fill = idents_res0.8)) +
329 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
330 #scale_color_gradientn(colours=pal(5))+
331 scale_color_brewer("Label")+
332 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
333 coord_equal() +
334 labs(title = "Integration", x = "UMAP1", y = "UMAP2") +
335 scale_fill_manual(values = cols,
336 guide = guide_legend(override.aes = list(size = 3),
337 keywidth = 0.1,
338 keyheight = 0.15,
339 default.unit = "inch"))+
340 #theme_bw()
341 ggplot_config(base.size = 6)
342 #dev.off()
343 ggsave(filename = paste0(merge_savePath,"UMAP_2.png"), p,
344 width = 10, height = 8, dpi = 800)
345
346 p <- ggplot(expr_all_manifast,aes(x =UMAP_x,y =UMAP_y ,fill = source)) +
347 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
348 #scale_color_gradientn(colours=pal(5))+
349 scale_color_brewer("Label")+
350 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
351 coord_equal() +
352 labs(title = "Integration", x = "UMAP1", y = "UMAP2") +
353 scale_fill_manual(values = c("blue","red"),
354 guide = guide_legend(override.aes = list(size = 3),
355 keywidth = 0.1,
356 keyheight = 0.15,
357 default.unit = "inch"))+
358 #theme_bw()
359 ggplot_config(base.size = 6)
360 #dev.off()
361 ggsave(filename = paste0(merge_savePath,"source.png"), p,
362 width = 10, height = 8, dpi = 1000)
363
364 #MarkerPlot
365
366 #Feature Plot
367
368 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),PV_gl)
369 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
370 rel_vst_ct[which(rel_vst_ct>4)] = 4
371
372 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], t(rel_vst_ct))
373 pal <- colorRampPalette(c("mediumseagreen", "lightyellow2", "deeppink"))
374
375 png(filename=paste0(merge_savePath,"PV_marker_plot_2.png"), width=5*400, height=400*(length(genetitle)/5+1))
376 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
377 ggplot(pltdat,aes(x =UMAP_x,y =UMAP_y ,fill = pltdat[,2+x])) +
378 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
379 scale_fill_gradientn(colours=pal(5))+
380 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
381 coord_equal() +
382 labs(title = genetitle[x], x = "UMAP1", y = "UMAP2") +
383 guides(fill = guide_colorbar(title = genetitle[x])) +
384 ggplot_config(base.size = 6)
385 })
386 grid.arrange(grobs = pp, ncol = 4)
387 dev.off()
388
389 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),CV_gl)
390 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
391 rel_vst_ct[which(rel_vst_ct>4)] = 4
392
393 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], t(rel_vst_ct))
394
395 png(filename=paste0(merge_savePath,"CV_marker_plot_2.png"), width=5*400, height=400*(length(genetitle)/5+1))
396 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
397 ggplot(pltdat,aes(x =UMAP_x,y =UMAP_y ,fill = pltdat[,2+x])) +
398 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
399 scale_fill_gradientn(colours=pal(5))+
400 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
401 coord_equal() +
402 labs(title = genetitle[x], x = "UMAP1", y = "UMAP2") +
403 guides(fill = guide_colorbar(title = genetitle[x])) +
404 ggplot_config(base.size = 6)
405 })
406 grid.arrange(grobs = pp, ncol = 4)
407 dev.off()
408
409 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),paper_gl)
410 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
411 rel_vst_ct[which(rel_vst_ct>4)] = 4
412
413 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], t(rel_vst_ct))
414
415 png(filename=paste0(merge_savePath,"paper_marker_plot_2.png"), width=5*400, height=400*(length(genetitle)/5+1))
416 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
417 ggplot(pltdat,aes(x =UMAP_x,y =UMAP_y ,fill = pltdat[,2+x])) +
418 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
419 scale_fill_gradientn(colours=pal(5))+
420 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
421 coord_equal() +
422 labs(title = genetitle[x], x = "UMAP1", y = "UMAP2") +
423 guides(fill = guide_colorbar(title = genetitle[x])) +
424 ggplot_config(base.size = 6)
425 })
426 grid.arrange(grobs = pp, ncol = 4)
427 dev.off()
428
429 boxplot(Cdh1~idents_res0.8, data = expr_all_manifast, xlab = "Cdh1",
430 ylab = "scale data")
431
432 #old seri
433 #rawcount <- SC_obj@assays$RNA@counts[which(rownames( SC_obj@assays$RNA@counts)%in%union(gene_plot,VariableFeatures(SC_obj))),]
434 #count <- as.matrix(rawcount)
435 #vst_ct <- var_stabilize(count)
436
437 gene_plot <- CV_gl
438 gene_plot <- unique(gene_plot)
439 expr_all <- ScaleData(expr_all)
440 vst_ct <- expr_all@assays$integrated@scale.data
441 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
442 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
443 rel_vst_ct <- t(sig_vst_ct)
444 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], rel_vst_ct)
445 colnames(pltdat)[1:2] <- c("x","y")
446 genetitle <- gene_plot
447 png(filename=paste0(merge_savePath,"CV_marker_plot.png"), width=5*400, height=400*(length(gene_plot)/5+1))
448 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
449 pattern_plot2(pltdat, x ,main = T, titlesize = 1.5, title = genetitle[x])
450 })
451 grid.arrange(grobs = pp, ncol = 4)
452 dev.off()
453 gene_plot <- PV_gl
454 gene_plot <- unique(gene_plot)
455 #vst_ct <- expr_all@assays$integrated@scale.data
456 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
457 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
458 rel_vst_ct <- t(sig_vst_ct)
459 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], rel_vst_ct)
460 colnames(pltdat)[1:2] <- c("x","y")
461 genetitle <- gene_plot
462 png(filename=paste0(merge_savePath,"PV_marker_plot.png"), width=5*400, height=400*(length(gene_plot)/5+1))
463 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
464 pattern_plot2(pltdat, x ,main = T, titlesize = 1.5, title = genetitle[x])
465 })
466 grid.arrange(grobs = pp, ncol = 4)
467 dev.off()
468
469 gene_plot <- paper_gl
470 gene_plot <- unique(gene_plot)
471 #vst_ct <- expr_all@assays$integrated@scale.data
472 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
473 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
474 rel_vst_ct <- t(sig_vst_ct)
475 pltdat <- cbind.data.frame(expr_all_manifast[, c("UMAP_x","UMAP_y")], rel_vst_ct)
476 colnames(pltdat)[1:2] <- c("x","y")
477 genetitle <- gene_plot
478 png(filename=paste0(merge_savePath,"Paper_marker_plot.png"), width=5*400, height=400*(length(gene_plot)/5+1))
479 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
480 pattern_plot2(pltdat, x ,main = T, titlesize = 1.5, title = genetitle[x])
481 })
482 grid.arrange(grobs = pp, ncol = 4)
483 dev.off()
484
485 #-------------Markers Plot----------------------
486 Markers_expr <- FindAllMarkers(expr_all, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
487 saveRDS(Markers_expr,file = paste0(merge_savePath,"Markers.RDS"))
488 top5_sc <- Markers_expr %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
489 gl <- top5_sc$gene
490 #expr_all <- ScaleData(expr_all)
491 sc_matrix <- expr_all@assays$integrated@scale.data
492 show_sc_matrix <- cbind(expr_all_manifast,as.numeric(expr_all_manifast$idents_res0.8))
493 colnames(show_sc_matrix) <- c(colnames(expr_all_manifast),"idents_sort")
494 show_sc_matrix <- cbind(show_sc_matrix,t(sc_matrix[intersect(gl,rownames(sc_matrix)),]))
495 show_sc_matrix <- show_sc_matrix[sort(show_sc_matrix$idents_sort,index.return=TRUE)$ix,]
496 show_sc <- t(as.matrix(show_sc_matrix)[,9:length(colnames(show_sc_matrix))])
497 show_sc <- apply(show_sc,2,as.numeric)
498 rownames(show_sc) <- colnames(show_sc_matrix)[9:length(colnames(show_sc_matrix))]
499 show_sc[which(show_sc>4)] <- 4
500 show_sc[which(show_sc<(-3))] <- -3
501 annocol_cs <- data.frame(idents = as.character(show_sc_matrix$idents_sort-1))
502 rownames(annocol_cs) <- (show_sc_matrix$barcode)
503 cols <- array(cols)
504 rownames(cols) <- as.character(0:(length(cols)-1))
505 ann_colors_sc = list(
506 idents = cols
507 )
508 png(filename=paste0(merge_savePath,"Diffgene.png"), width=1200, height=900)
509 pheatmap(show_sc, annotation_col =annocol_cs,cluster_rows =F,cluster_cols =F,show_colnames = F,annotation_colors = ann_colors_sc)
510 dev.off()
511 #
512 pal <- colorRampPalette(c("mediumseagreen", "lightyellow2", "deeppink"))
513 p <- ggplot(expr_all_manifast,aes(x =UMAP_x,y =UMAP_y ,fill = expr_all@meta.data$nCount_RNA)) +
514 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
515 scale_fill_gradientn(colours=pal(5))+
516 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
517 coord_equal() +
518 labs(title = "nCount_RNA", x = "UMAP1", y = "UMAP2") +
519 guides(fill = guide_colorbar(title = "nCount")) +
520 ggplot_config(base.size = 6)
521
522 ggsave(filename = paste0(merge_savePath,"nCount_RNA.png"), p,
523 width = 10, height = 8, dpi = 1000)
524
525 #expr_all[["mito.percent"]] <- PercentageFeatureSet(expr_all, pattern = "^mt-")
526
527 p <- ggplot(expr_all_manifast,aes(x =UMAP_x,y =UMAP_y ,fill = expr_all@meta.data$mito.percent)) +
528 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
529 scale_fill_gradientn(colours=pal(5))+
530 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
531 coord_equal() +
532 labs(title = "mito.percent", x = "UMAP1", y = "UMAP2") +
533 guides(fill = guide_colorbar(title = "mito")) +
534 ggplot_config(base.size = 6)
535
536 ggsave(filename = paste0(merge_savePath,"mito.percent.png"), p,
537 width = 10, height = 8, dpi = 1000)
538
539 ##########
540 PV_gl <- intersect(PV_gl,rownames(sc_matrix))
541 CV_gl <- intersect(CV_gl,rownames(sc_matrix))
542
543 for(i in 1:length(PV_gl)){
544 for(j in 1:length(CV_gl)){
545 gene1 <- PV_gl[i]
546 gene2 <- CV_gl[j]
547 if(gene1!=gene2){
548 gene1_expr <- (sc_matrix[gene1,])
549 gene2_expr <- (sc_matrix[gene2,])
550 sc_id <- expr_all_manifast$idents_res0.8
551 sc_gp_df <- data.frame(gene1_expr = gene1_expr,gene2_expr = gene2_expr,sc_id = sc_id)
552 png(filename=paste0(merge_savePath,"/corrGene/",gene1,"_",gene2,"_sc.png"), width=900, height=900)
553 plot(sc_gp_df$gene1_expr,sc_gp_df$gene2_expr, cex= 1, col="black", pch=21, bg=cols[sc_gp_df$sc_id],xlab = gene1,ylab = gene2)
554 #text(sc_gp_df$gene1_expr,sc_gp_df$gene2_expr, 1:length(sc_gp_df$gene2_expr), cex=0.9)
555 dev.off()
556 }
557 }
558 }
559 #_-----------Evolution
560
561 ##PCA plot
562
563 #plot(expr_all_manifast$UMAP_x,expr_all_manifast$UMAP_y, pch=19, col= cols[expr_all_manifast$idents_res0.8],xlab = "UMAP_1",ylab = "UMAP_2")
564 p <- ggplot(expr_all_manifast,aes(x = PCA_x,y = PCA_y ,fill = idents_res0.8)) +
565 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
566 #scale_color_gradientn(colours=pal(5))+
567 scale_color_brewer("Label")+
568 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
569 coord_equal() +
570 labs(title = "Integration", x = "CD1", y = "CD2") +
571 scale_fill_manual(values = cols,
572 guide = guide_legend(override.aes = list(size = 3),
573 keywidth = 0.1,
574 keyheight = 0.15,
575 default.unit = "inch"))+
576 #theme_bw()
577 ggplot_config(base.size = 6)
578 #dev.off()
579 ggsave(filename = paste0(merge_savePath,"PCA.png"), p,
580 width = 10, height = 10, dpi = 1000)
581
582 p <- ggplot(expr_all_manifast[which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7)),],aes(x =PCA_x,y =PCA_y ,fill = idents_res0.8)) +
583 geom_point(shape = 21, size = 2, stroke = 0.2, color = "lightgrey") +
584 #scale_color_gradientn(colours=pal(5))+
585 scale_color_brewer("Label")+
586 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
587 coord_equal() +
588 labs(title = "Integration", x = "CD1", y = "CD2") +
589 scale_fill_manual(values = cols,
590 guide = guide_legend(override.aes = list(size = 3),
591 keywidth = 0.1,
592 keyheight = 0.15,
593 default.unit = "inch"))+
594 #theme_bw()
595 ggplot_config(base.size = 6)
596
597 ggsave(filename = paste0(merge_savePath,"PCA.png"), p,
598 width = 10, height = 10, dpi = 1000)
599 ############monocel for cluster 0,1,2,3,4,5,6,7,8,9
600 library(monocle)
601 data <- as(as.matrix(expr_all@assays$RNA@counts[,which(expr_all@meta.data$RNA_snn_res.0.8%in%c(0,1,2,3,4,5,6,7,8,9))]), 'sparseMatrix')
602 pd <- new('AnnotatedDataFrame', data = expr_all@meta.data[which(expr_all@meta.data$RNA_snn_res.0.8%in%c(0,1,2,3,4,5,6,7,8,9)),] )
603 fData <- data.frame(gene_short_name = rownames(expr_all@assays$RNA@counts), row.names = rownames(expr_all@assays$RNA@counts))
604 fd <- new('AnnotatedDataFrame', data = fData)
605
606
607 data <- as(as.matrix(expr_all@assays$RNA@counts), 'sparseMatrix')
608 pd <- new('AnnotatedDataFrame', data = expr_all@meta.data )
609
610 #Construct monocle cds
611 monocle_cds <- newCellDataSet(data,
612 phenoData = pd,
613 featureData = fd,
614 lowerDetectionLimit = 0.5,
615 expressionFamily = negbinomial.size())
616 HSMM <- monocle_cds
617 HSMM <- estimateSizeFactors(HSMM)
618 HSMM <- estimateDispersions(HSMM)
619 HSMM <- detectGenes(HSMM, min_expr = 3 )
620 print(head(fData(HSMM)))
621 #expressed_genes <- row.names(subset(fData(HSMM),num_cells_expressed >= 10))
622 expressed_genes <- union(PV_gl,union(CV_gl,paper_gl))
623 HSMM <- setOrderingFilter(HSMM, expressed_genes)
624 plot_ordering_genes(HSMM)
625 HSMM <- reduceDimension(HSMM, max_components = 2,reduction_method = 'DDRTree')
626 HSMM <- orderCells(HSMM)
627
628
629 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),PV_gl)
630 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
631 rel_vst_ct[which(rel_vst_ct>4)] = 4
632 pltdat <- cbind.data.frame(t(HSMM@reducedDimS[,]), t(rel_vst_ct))
633 colnames(pltdat)[1:2] <- c("X","Y")
634
635 png(filename=paste0(merge_savePath,"PV_monocle_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
636 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
637 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
638 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
639 scale_fill_gradientn(colours=pal(5))+
640 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
641 coord_equal() +
642 labs(title = genetitle[x], x = "x", y = "y") +
643 guides(fill = guide_colorbar(title = genetitle[x])) +
644 ggplot_config(base.size = 6)
645 })
646 grid.arrange(grobs = pp, ncol = 4)
647 dev.off()
648
649 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),CV_gl)
650 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
651 rel_vst_ct[which(rel_vst_ct>4)] = 4
652
653 pltdat <- cbind.data.frame(t(HSMM@reducedDimS[,]), t(rel_vst_ct))
654 colnames(pltdat)[1:2] <- c("X","Y")
655
656 png(filename=paste0(merge_savePath,"CV_monocle_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
657 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
658 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
659 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
660 scale_fill_gradientn(colours=pal(5))+
661 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
662 coord_equal() +
663 labs(title = genetitle[x], x = "x", y = "y") +
664 guides(fill = guide_colorbar(title = genetitle[x])) +
665 ggplot_config(base.size = 6)
666 })
667 grid.arrange(grobs = pp, ncol = 4)
668 dev.off()
669
670 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),paper_gl)
671 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
672 rel_vst_ct[which(rel_vst_ct>4)] = 4
673
674 pltdat <- cbind.data.frame(t(HSMM@reducedDimS[,]), t(rel_vst_ct))
675 colnames(pltdat)[1:2] <- c("X","Y")
676
677
678 png(filename=paste0(merge_savePath,"paper_monocle_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
679 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
680 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
681 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
682 scale_fill_gradientn(colours=pal(5))+
683 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
684 coord_equal() +
685 labs(title = genetitle[x], x = "x", y = "y") +
686 guides(fill = guide_colorbar(title = genetitle[x])) +
687 ggplot_config(base.size = 6)
688 })
689 grid.arrange(grobs = pp, ncol = 4)
690 dev.off()
691
692 boxplot(Cdh1~idents_res0.8, data = expr_all_manifast, xlab = "Cdh1",
693 ylab = "scale data")
694
695
696 ##PCA
697 gene_plot <- CV_gl
698 gene_plot <- unique(gene_plot)
699 expr_all <- ScaleData(expr_all)
700 vst_ct <- expr_all@assays$RNA@scale.data
701 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
702 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
703 rel_vst_ct <- t(sig_vst_ct)
704 pltdat <- cbind.data.frame(expr_all_manifast[, c("PCA_x","PCA_y")], rel_vst_ct)
705 colnames(pltdat)[1:2] <- c("x","y")
706 genetitle <- gene_plot
707 png(filename=paste0(merge_savePath,"CV_marker_plot_PCA.png"), width=5*400, height=400*(length(gene_plot)/5+1))
708 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
709 pattern_plot2(pltdat[which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7)),], x ,main = T, titlesize = 1.5, title = genetitle[x])
710 })
711 grid.arrange(grobs = pp, ncol = 4)
712 dev.off()
713 gene_plot <- PV_gl
714 gene_plot <- unique(gene_plot)
715 vst_ct <- expr_all@assays$RNA@scale.data
716 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
717 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
718 rel_vst_ct <- t(sig_vst_ct)
719 pltdat <- cbind.data.frame(expr_all_manifast[, c("PCA_x","PCA_y")], rel_vst_ct)
720 colnames(pltdat)[1:2] <- c("x","y")
721 genetitle <- gene_plot
722 png(filename=paste0(merge_savePath,"PV_marker_plot_PCA.png"), width=5*400, height=400*(length(gene_plot)/5+1))
723 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
724 pattern_plot2(pltdat[which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7)),], x ,main = T, titlesize = 1.5, title = genetitle[x])
725 })
726 grid.arrange(grobs = pp, ncol = 4)
727 dev.off()
728 gene_plot <- paper_gl
729 gene_plot <- unique(gene_plot)
730 expr_all <- ScaleData(expr_all)
731 #vst_ct <- expr_all@assays$RNA@scale.data
732 sig_vst_ct <- vst_ct[intersect(gene_plot,rownames(vst_ct)), ]
733 #rel_vst_ct <- apply(sig_vst_ct, 1, relative_func)
734 rel_vst_ct <- t(sig_vst_ct)
735 pltdat <- cbind.data.frame(expr_all_manifast[, c("PCA_x","PCA_y")], rel_vst_ct)
736 colnames(pltdat)[1:2] <- c("x","y")
737 genetitle <- gene_plot
738 png(filename=paste0(merge_savePath,"paper_gl_marker_plot_PCA.png"), width=5*400, height=400*(length(gene_plot)/5+1))
739 pp <- lapply(1:(ncol(pltdat) - 2), function(x) {
740 pattern_plot2(pltdat[which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7)),], x ,main = T, titlesize = 1.5, title = genetitle[x])
741 })
742 grid.arrange(grobs = pp, ncol = 4)
743 dev.off()
744 ##diffusion map
745 library(destiny)
746
747 res <- DiffusionMap(data = t(as.matrix(expr_all@assays$RNA@data[VariableFeatures(expr_all),which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7))])))
748
749 expr_genelist <- intersect(union(PV_gl,union(CV_gl,paper_gl)),rownames(expr_all))
750 res <- DiffusionMap(data = t(as.matrix(expr_all@assays$RNA@data[expr_genelist,which(expr_all_manifast$idents_res0.8%in%c(0,1,2,3,4,6,7))])))
751
752 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),PV_gl)
753 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
754 rel_vst_ct[which(rel_vst_ct>4)] = 4
755 pltdat <- cbind.data.frame((res@eigenvectors[,1:2]*20), t(rel_vst_ct))
756 colnames(pltdat)[1:2] <- c("X","Y")
757
758 png(filename=paste0(merge_savePath,"PV_DiffusionMap_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
759 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
760 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
761 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
762 scale_fill_gradientn(colours=pal(5))+
763 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
764 coord_equal() +
765 labs(title = genetitle[x], x = "x", y = "y") +
766 guides(fill = guide_colorbar(title = genetitle[x])) +
767 ggplot_config(base.size = 6)
768 })
769 grid.arrange(grobs = pp, ncol = 4)
770 dev.off()
771
772 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),CV_gl)
773 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
774 rel_vst_ct[which(rel_vst_ct>4)] = 4
775
776 pltdat <- pltdat <- cbind.data.frame((res@eigenvectors[,1:2]*20), t(rel_vst_ct))
777 colnames(pltdat)[1:2] <- c("X","Y")
778
779 png(filename=paste0(merge_savePath,"CV_DiffusionMap_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
780 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
781 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
782 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
783 scale_fill_gradientn(colours=pal(5))+
784 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
785 coord_equal() +
786 labs(title = genetitle[x], x = "x", y = "y") +
787 guides(fill = guide_colorbar(title = genetitle[x])) +
788 ggplot_config(base.size = 6)
789 })
790 grid.arrange(grobs = pp, ncol = 4)
791 dev.off()
792
793 genetitle <- intersect(rownames(expr_all@assays$RNA@scale.data),paper_gl)
794 rel_vst_ct <- expr_all@assays$RNA@scale.data[genetitle,]
795 rel_vst_ct[which(rel_vst_ct>4)] = 4
796
797 pltdat <- pltdat <- cbind.data.frame((res@eigenvectors[,1:2]*20), t(rel_vst_ct))
798 colnames(pltdat)[1:2] <- c("X","Y")
799
800 png(filename=paste0(merge_savePath,"paper_DiffusionMap_plot.png"), width=5*400, height=400*(length(genetitle)/5+1))
801 pp <- lapply(1:(ncol(pltdat)-2), function(x) {
802 ggplot(pltdat,aes(x =X,y =Y ,fill = pltdat[,2+x])) +
803 geom_point(shape = 21, size = 1, stroke = 0.2, color = "lightgrey") +
804 scale_fill_gradientn(colours=pal(5))+
805 scale_x_discrete(expand = c(xpand, ypand)) + scale_y_discrete(expand = c(xpand, ypand)) +
806 coord_equal() +
807 labs(title = genetitle[x], x = "x", y = "y") +
808 guides(fill = guide_colorbar(title = genetitle[x])) +
809 ggplot_config(base.size = 6)
810 })
811 grid.arrange(grobs = pp, ncol = 4)
812 dev.off()
813
814
815 #Cell cycle
816
817 cell_cycle_gl <- c("Foxa1", "Foxa2", "Gata4", "Gata6", "vHnf1",
818 "Hex", "Tbx3", "Prox1", "Hnf6/OC-1", "OC-2",
819 "Hnf4",
820 "Hnf6/OC-1",
821 "Hex", "Hnf6/OC-1", "Ptf1a",
822 "Foxa1", "Foxa2", "Hnf6/OC-1", "Hlxb9", "Ptf1a",
823 "Pdx1", "Ptf1a", "Hes1", "Rbpj", "Sox9", "Cpa1",
824 "Ptf1a", "Rpbji",
825 "Hnf6/OC-1",
826 "lsl1", "Ngn3", "NeuroD", "lnsm1",
827 "Mafa", "Pdx1", "Hlxb9", "Pax4", "Pax6", "lsl1", "Nkx2.2", "Nkx6.1",
828 "Foxa2", "Nkx2.2")
829
830 S1_Negative_C <- CellCycleScoring(S1_Negative, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
831
832 # view cell cycle scores and phase assignments
833 head(marrow[[]])
834
835
836
837 #RNA velocity
838 library(Seurat)
839 library(velocyto.R)
840 library(SeuratWrappers)
841
842
843 ##############################################################################################
844 # Functions #
845 ##############################################################################################
846 filter_str <- function(data,strsep = '_',filter_id = 1){
847 filter_barcode <- strsplit(data,strsep)
848 filter_barcode_ <- rep(1,length(filter_barcode))
849 for(i in 1:length(filter_barcode_)){
850 filter_barcode_[i] <- filter_barcode[[i]][filter_id]
851 }
852 return(filter_barcode_)
853 }
854 pattern_plot2 <- function(pltdat, igene, xy = T, main = F, titlesize = 2,
855 pointsize = 1, xpand = 0, ypand = 1, title = NULL) {
856 if (!xy) {
857 xy <- matrix(as.numeric(do.call(rbind, strsplit(as.character(pltdat[,
858 1]), split = "x"))), ncol = 2)
859 rownames(xy) <- as.character(pltdat[, 1])
860 colnames(xy) <- c("x", "y")
861 pd <- cbind.data.frame(xy, pltdat[, 2:ncol(pltdat)])
862 } else {
863 pd <- pltdat
864 }
865
866 # pal <- colorRampPalette(c('seagreen1','orangered')) pal <-
867 # colorRampPalette(c('#00274c','#ffcb05')) pal <-
868 # colorRampPalette(c('deepskyblue','goldenrod1')) pal <-
869 # colorRampPalette(c('deepskyblue','deeppink'))
870 pal <- colorRampPalette(c("mediumseagreen", "lightyellow2", "deeppink"))
871 #pal <- colorRampPalette(c("blue", "lightyellow2", "red"))
872 gpt <- ggplot(pd, aes(x = x, y = y, color = pd[, igene + 2])) + geom_point(size = pointsize) +
873 # scale_color_gradientn(colours=pal(5))+
874 scale_color_gradientn(colours = pal(5)) + scale_x_discrete(expand = c(xpand,
875 ypand)) + scale_y_discrete(expand = c(xpand, ypand)) + coord_equal() +
876 # labs(title = colnames(pd)[igene+2], x = NULL, y = NULL)+
877 theme_bw()
878 if (main) {
879 if (is.null(title)) {
880 title = colnames(pd)[igene + 2]
881 }
882 out = gpt + labs(title = title, x = NULL, y = NULL) + theme(legend.position = "none",
883 plot.title = element_text(hjust = 0.5, size = rel(titlesize)))
884 } else {
885 out = gpt + labs(title = NULL, x = NULL, y = NULL) + theme(legend.position = "none")
886 }
887 return(out)
888 }
889 pattern_plot1 <- function(pltdat, igene, xy = T, main = F, titlesize = 2,
890 pointsize = 1, xpand = 0, ypand = 1, title = NULL) {
891 if (!xy) {
892 xy <- matrix(as.numeric(do.call(rbind, strsplit(as.character(pltdat[,
893 1]), split = "x"))), ncol = 2)
894 rownames(xy) <- as.character(pltdat[, 1])
895 colnames(xy) <- c("x", "y")
896 pd <- cbind.data.frame(xy, pltdat[, 2:ncol(pltdat)])
897 } else {
898 pd <- pltdat
899 }
900
901 # pal <- colorRampPalette(c('seagreen1','orangered')) pal <-
902 # colorRampPalette(c('#00274c','#ffcb05')) pal <-
903 # colorRampPalette(c('deepskyblue','goldenrod1')) pal <-
904 # colorRampPalette(c('deepskyblue','deeppink'))
905 #pal <- colorRampPalette(c("mediumseagreen", "lightyellow2", "deeppink"))
906 pal <- colorRampPalette(c("blue", "lightyellow2", "red"))
907 gpt <- ggplot(pd, aes(x = x, y = y, color = pd[, igene + 2])) + geom_point(size = pointsize) +
908 # scale_color_gradientn(colours=pal(5))+
909 scale_color_gradientn(colours = pal(5)) + scale_x_discrete(expand = c(xpand,
910 ypand)) + scale_y_discrete(expand = c(xpand, ypand)) + coord_equal() +
911 # labs(title = colnames(pd)[igene+2], x = NULL, y = NULL)+
912 theme_bw()
913 if (main) {
914 if (is.null(title)) {
915 title = colnames(pd)[igene + 2]
916 }
917 out = gpt + labs(title = title, x = NULL, y = NULL) + theme(legend.position = "none",
918 plot.title = element_text(hjust = 0.5, size = rel(titlesize)))
919 } else {
920 out = gpt + labs(title = NULL, x = NULL, y = NULL) + theme(legend.position = "none")
921 }
922 return(out)
923 }
924
925 convertMouseGeneList <- function(x){
926 require("biomaRt")
927 human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
928 mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
929
930 genesV2 = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = x , mart = mouse, attributesL = c("hgnc_symbol"), martL = human, uniqueRows=T)
931 #humanx <- unique(genesV2[, 2])
932
933 # Print the first 6 genes found to the screen
934 #print(head(humanx))
935 return(genesV2)
936 }
937
938 getDefaultColors <- function(n = NULL, type = 1){
939 if(type == 1){
940 colors <- c("#cb7c77", "#68d359", "#6a7dc9", "#c9d73d", "#c555cb",
941 "#d7652d", "#7cd5c8", "#c49a3f", "#507d41", "#5d8d9c",
942 "#90353b", "#674c2a", "#1B9E77", "#c5383c", "#0081d1",
943 "#ffd900", "#502e71", "#c8b693", "#aed688", "#f6a97a",
944 "#c6a5cc", "#798234", "#6b42c8", "#cf4c8b", "#666666",
945 "#feb308", "#ff1a1a", "#1aff1a", "#1a1aff", "#ffff1a")
946 }else if(type == 2){
947 if(n <= 8){
948 colors <- c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
949 "#A6D854", "#FFD92F", "#E5C494", "#B3B3B3")
950 }else if(n <= 14){
951 colors <- c("#437BFE", "#FEC643", "#43FE69", "#FE6943", "#C643FE",
952 "#43D9FE", "#B87A3D", "#679966", "#993333", "#7F6699",
953 "#E78AC3", "#333399", "#A6D854", "#E5C494")
954 }
955 else if(n <= 20){
956 colors <- c("#87b3d4", "#d5492f", "#6bd155", "#683ec2", "#c9d754",
957 "#d04dc7", "#81d8ae", "#d34a76", "#607d3a", "#6d76cb",
958 "#ce9d3f", "#81357a", "#d3c3a4", "#3c2f5a", "#b96f49",
959 "#4e857e", "#6e282c", "#d293c8", "#393a2a", "#997579")
960 }else if(n <= 30){
961 colors <- c("#628bac", "#ceda3f", "#7e39c9", "#72d852", "#d849cc",
962 "#5e8f37", "#5956c8", "#cfa53f", "#392766", "#c7da8b",
963 "#8d378c", "#68d9a3", "#dd3e34", "#8ed4d5", "#d84787",
964 "#498770", "#c581d3", "#d27333", "#6680cb", "#83662e",
965 "#cab7da", "#364627", "#d16263", "#2d384d", "#e0b495",
966 "#4b272a", "#919071", "#7b3860", "#843028", "#bb7d91")
967 }else{
968 colors <- c("#982f29", "#5ddb53", "#8b35d6", "#a9e047", "#4836be",
969 "#e0dc33", "#d248d5", "#61a338", "#9765e5", "#69df96",
970 "#7f3095", "#d0d56a", "#371c6b", "#cfa738", "#5066d1",
971 "#e08930", "#6a8bd3", "#da4f1e", "#83e6d6", "#df4341",
972 "#6ebad4", "#e34c75", "#50975f", "#d548a4", "#badb97",
973 "#b377cf", "#899140", "#564d8b", "#ddb67f", "#292344",
974 "#d0cdb8", "#421b28", "#5eae99", "#a03259", "#406024",
975 "#e598d7", "#343b20", "#bbb5d9", "#975223", "#576e8b",
976 "#d97f5e", "#253e44", "#de959b", "#417265", "#712b5b",
977 "#8c6d30", "#a56c95", "#5f3121", "#8f846e", "#8f5b5c")
978 }
979 }else if(type == 3){
980 # colors <- c("#07a2a4", "#9a7fd1", "#588dd5", "#f5994e",
981 # "#c05050", "#59678c", "#c9ab00", "#7eb00a")
982 colors <- c("#c14089", "#6f5553", "#E5C494", "#738f4c", "#bb6240",
983 "#66C2A5", "#2dfd29", "#0c0fdc")
984 }
985 if(!is.null(n)){
986 if(n <= length(colors)){
987 colors <- colors[1:n]
988 }else{
989 step <- 16777200 %/% (n - length(colors)) - 2
990 add.colors <- paste0("#", as.hexmode(seq(from = sample(1:step, 1),
991 by = step, length.out = (n-length(colors)))))
992 colors <- c(colors, add.colors)
993 }
994 }
995 return(colors)
996 }
997
998 ggplot_config <- function(base.size = 8){
999 p <- theme_classic() +
1000 theme(plot.title = element_text(size = 2 * base.size),
1001 axis.title.x = element_text(size = 2 * base.size, vjust = -0.2),
1002 axis.title.y = element_text(size = 2 * base.size, vjust = 0.2),
1003 axis.text.x = element_text(size = 2 * base.size),
1004 axis.text.y = element_text(size = 2 * base.size),
1005 panel.grid.major = element_blank(),
1006 panel.grid.minor = element_blank(),
1007 legend.title = element_text(size = 2 * base.size - 2),
1008 legend.text = element_text(size = 1.5 * base.size)
1009 )
1010 return(p)
1011 }
