rm(list = ls())
library(Seurat)
library(dplyr)
library(patchwork)
library(ggplot2)
load( "sce.Rdata")
scRNA = sce
scRNA@meta.data$celltype = Idents(scRNA)
ctys = levels(scRNA)
ctys
## [1] "naive B" "CD8 T" "Naive CD4 T" "plasma B" "CD14+ Mono"
## [6] "endothelial" "Fibroblasts" "NK" "DC"
scRNA.markers <- FindAllMarkers(scRNA, min.pct = 0.25,
logfc.threshold = 0.25)
head(scRNA.markers)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## CD79A 0.000000e+00 5.356611 0.915 0.060 0.000000e+00 naive B CD79A
## BANK1 0.000000e+00 6.778892 0.853 0.026 0.000000e+00 naive B BANK1
## MS4A1 0.000000e+00 5.677155 0.848 0.050 0.000000e+00 naive B MS4A1
## HLA-DRA 1.482079e-306 2.775984 0.994 0.336 3.060197e-302 naive B HLA-DRA
## HLA-DQB1 7.341732e-295 2.622558 0.926 0.207 1.515921e-290 naive B HLA-DQB1
## HLA-DQA1 5.404715e-292 2.727477 0.865 0.138 1.115966e-287 naive B HLA-DQA1
colnames(scRNA.markers)[6] = "celltype"
k = scRNA.markers$p_val_adj<0.05;table(k)
## k
## FALSE TRUE
## 3960 8946
scRNA.markers = scRNA.markers[k,]
#上下调
scRNA.markers$label <- ifelse(scRNA.markers$avg_log2FC<0,"sigDown","sigUp")
topgene <- scRNA.markers %>%
group_by(celltype) %>%
top_n(n = 10, wt = avg_log2FC) %>%
bind_rows(group_by(scRNA.markers, celltype) %>%
top_n(n = 10, wt = -avg_log2FC))
head(topgene)
## # A tibble: 6 × 8
## # Groups: celltype [1]
## p_val avg_log2FC pct.1 pct.2 p_val_adj celltype gene label
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> <chr>
## 1 0 6.78 0.853 0.026 0 naive B BANK1 sigUp
## 2 3.13e-250 6.98 0.498 0.011 6.47e-246 naive B LINC00926 sigUp
## 3 2.59e-207 6.02 0.423 0.01 5.34e-203 naive B CD24 sigUp
## 4 1.42e-204 6.90 0.411 0.008 2.92e-200 naive B LINC02397 sigUp
## 5 1.15e-196 7.13 0.403 0.01 2.38e-192 naive B IGHD sigUp
## 6 5.77e-166 7.47 0.335 0.006 1.19e-161 naive B PAX5 sigUp
#根据log2FC范围确定背景柱长度：
dfbar = scRNA.markers %>%
group_by(celltype) %>%
summarise(low = round(min(avg_log2FC)-0.5),
up = round(max(avg_log2FC)+0.5))
#绘制背景柱和散点图：
p1 <- ggplot()+
geom_col(aes(x = celltype ,y = low),dfbar,
fill = "#dcdcdc",alpha = 0.6)+
geom_col(aes(x = celltype ,y = up),dfbar,
fill = "#dcdcdc",alpha = 0.6)+
geom_jitter(aes(x = celltype, y = avg_log2FC, color = label),scRNA.markers,
width =0.4,size = 1)+
scale_color_manual(values = c("#0077c0","#c72d2e"))+
theme_classic()
p1
#X轴的色块标签：
library(RColorBrewer)
mycol <- colorRampPalette(rev(brewer.pal(n = 7, name ="Set1")))(length(ctys))
p2 <- p1 +
geom_tile(aes(x = ctys,y = 0),
height = 0.5,fill = mycol, show.legend = F)+
geom_text(aes(x= ctys, y = 0, label = ctys),
size = 3,fontface = "bold")
p2
library(ggrepel)
#给每个Cluster差异表达top基因加上标签,调整细节：
p3 <- p2 +
geom_text_repel(aes(x = celltype,y = avg_log2FC,label = gene),
topgene,size = 3 )+
labs(x = "CellType",y = "Average log2FoldChange",
title = "Differential expression genes")+
theme(
plot.title = element_text(size = 14,color = "black",face = "bold"),
axis.title = element_text(size = 12,color = "black",face = "bold"),
axis.line.y = element_line(color = "black",linewidth = 0.8),
axis.line.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid = element_blank(),
legend.position = c(0.98,0.96),
legend.background = element_blank(),
legend.title = element_blank(),
legend.direction = "vertical",
legend.justification = c(1,0),
legend.text = element_text(size = 12)
)+
guides(color = guide_legend(override.aes = list(size = 4)))
p3
