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 }