scRNA分析|一(尽)文(力)解决你的单细胞火山图问题

2022-11-03 15:32:39 浏览数 (1)

当有了聚类结果(cluster)或注释结果(celltype)后就可以 找不同cluster/celltype间,不同样本间 或者 不同分组间的差异,为后面的 机制探索 或者 样本间/组间异质性研究 提供一些帮助。

但是在对差异结果进行可视化的时候,你是否经常会遇到以下情况,

1)单细胞数据的火山图 中间经常会缺一块 ?Y轴最高处会堆积很多点?

2)除bulk常见的火山图外,文献中是否还有什么其他的展示形式?

3)可以同时展示单细胞各个cluster /celltype的差异情况吗?

一 差异marker分析

Seurat可以通过FindMarkers函数 和 FindAllMarkers函数寻找不同cluster的差异表达基因。

1 FindMarkers函数
1.1 特定cluster/celltype 与其他的差异

是one-others的差异分析方法,由ident.1来制定cluster/celltype ,然今后计算ident.1与其他所有的差异。本例就是Epi与其余的所有celltype找差异marker。

代码语言:javascript复制
library(Seurat) 
library(tidyverse)
#载入数据
load("sce.anno.RData")
#设置颜色
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",  
         "#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
         "#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
         "#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
# find all markers of cluster Epi
cluster1.markers <- FindMarkers(sce2, ident.1 = "Epi", min.pct = 0.25)
head(cluster1.markers, n = 5)
代码语言:javascript复制
#       p_val avg_log2FC pct.1 pct.2 p_val_adj
# MXRA8     0 -1.5057649 0.024 0.385         0
# MFAP2     0 -0.8893903 0.020 0.279         0
# C1QA      0 -3.0473042 0.048 0.333         0
# C1QC      0 -2.5871212 0.022 0.273         0
# C1QB      0 -3.0683349 0.043 0.316         0
1.2 指定cluster/celltype 之间的差异

示例为计算Epi 和 T B 细胞之间的差异

代码语言:javascript复制
cluster2.markers <- FindMarkers(sce2, ident.1 = "Epi", ident.2 = c("T", "B"), min.pct = 0.25)
head(cluster2.markers, n = 5)

FindMarkers函数常用参数:

min.pct: 基因在两个cluster/celltype中任何一个被检测到的百分比阈值,默认0.1 。 logfc.threshold :默认为0.25 thresh.test:设定基因在两个cluster/celltype中基因差异表达量。 test.use :差异marker的检验方式,默认为“wilcox”,可选“bimod”,“roc”,“t”,“negbinom”,“LR”,“MAST”,“DESeq2”等; only.pos:是否仅返回阳性标记(默认为FALSE) 更多参数可以使用 ??FindMarkers 查看。

2 FindAllMarkers 函数

每个cluster/celltype 分别与其他所有cluster/celltype 进行比较,展示各个cluster/celltype 的前2个基因

代码语言:javascript复制
sce2.markers <- FindAllMarkers(sce2, only.pos = TRUE, min.pct = 0.25, 
                               logfc.threshold = 0.25)
sce2.markers %>%
    group_by(cluster) %>%
    top_n(n = 2, wt = avg_log2FC)
代码语言:javascript复制
      p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
1 1.88e-117       1.08 0.913 0.588 2.57e-113 0       LDHB    
2 5.01e- 85       1.34 0.437 0.108 6.88e- 81 0       CCR7    
3 0               5.57 0.996 0.215 0         1       S100A9

注意FindAllMarkers比FindMarkers的结果多了cluster列,结果为该cluster与其他所有cluster相比的差异基因。

FindAllMarkers 与 FindMarkers 参数类似,但是需要注意only.pos的默认值,only.pos:仅返回阳性标记(默认为TRUE

二 ,可视化-火山图

1,常规方式绘制火山图
代码语言:javascript复制
#自定义阈值
log2FC = 1
padj = 0.05 

cluster1.markers$threshold="ns";
cluster1.markers[which(cluster1.markers$avg_log2FC  > log2FC & cluster1.markers$p_val_adj <padj),]$threshold="up";
cluster1.markers[which(cluster1.markers$avg_log2FC  < (-log2FC) & cluster1.markers$p_val_adj < padj),]$threshold="down";
cluster1.markers$threshold=factor(cluster1.markers$threshold, levels=c('down','ns','up'))

ggplot(data=cluster1.markers, aes(x=avg_log2FC, y=-log10(p_val_adj), color=threshold))  
  geom_point(alpha=0.8, size=0.8)  
  geom_vline(xintercept = c(-log2FC, log2FC), linetype=2, color="grey") 
  geom_hline(yintercept = -log10(padj), linetype=2, color="grey") 
  #labs(title= ifelse(""==title, "", paste("DEG:", title))) 
  xlab(bquote(Log[2]*FoldChange)) 
  ylab(bquote(-Log[10]*italic(P.adj)) ) 
  theme_classic(base_size = 14)  
  scale_color_manual('',labels=c(paste0("down(",table(cluster1.markers$threshold)[[1]],')'),'ns',
                                 paste0("up(",table(cluster1.markers$threshold)[[3]],')' )),
                     values=c("blue", "grey","red" ) ) 
  guides(color=guide_legend(override.aes = list(size=3, alpha=1)))

更多颜色,参数的设置以及如何添加感兴趣的基因详见ggplot2-plotly|让你的火山图“活”过来 和 火山图|给你geneList,帮我标到火山图上。

回答开始提出的问题1:

(1)因为单细胞自身区别于bulk数据的特异性,大概率会出现很多P值为0或者无限接近于0的基因,绘制常规火山图就会出现Y轴顶出现很多点。

(2)中间的空白是因为FindMarkers的默认参数含还有很多阈值,很多基因被过滤掉了。可以通过调整min.pct,logfc.threshold 和thresh.test等参数得到更多的基因,自行尝试比较。

2, 单细胞火山图展示方式

在生信技能树 两个不同单细胞亚群差异分析,何必一定要做火山图 中提到NC文章 《CD177 modulates the function and homeostasis of tumor-infiltrating regulatory T cells》有以下呈现方式,横坐标使用cluster/celltype间的pct差值,纵坐标使用log2FC ,P值可以使用颜色来表示。

文献阈值参数:Genes labeled have logfold change > 1, Δ Percentage Difference > 20% and adjusted P-value <0.05.

使用ggplot2 点图绘制方式

代码语言:javascript复制
library(ggrepel)
cluster1.markers <- cluster1.markers %>%
  mutate(Difference = pct.1 - pct.2) %>% 
  rownames_to_column("gene")
  
ggplot(cluster1.markers, aes(x=Difference, y=avg_log2FC, color = threshold))   
  geom_point(size=0.5)   
  scale_color_manual(values=c( "blue","grey","red") )   
  geom_label_repel(data=subset(cluster1.markers, avg_log2FC >= 1 & Difference >= 0.2 & p_val_adj <= 0.05), 
                   aes(label=gene),  #添加label
                   color="black", #设置label中标签的颜色
                   segment.colour = "black",#设置label框的颜色
                   label.padding = 0.1, 
                   #max.overlaps = 200,
                   segment.size = 0.3,  #框的大小
                   size=4) 
  geom_label_repel(data=subset(cluster1.markers, avg_log2FC <= -1 & Difference <= -0.2 & p_val_adj <= 0.05), 
                   aes(label=gene), label.padding = 0.1, 
                   color="black",
                   segment.colour = "black",
                   segment.size = 0.3, size=4) 
  geom_vline(xintercept = 0.0,linetype=2) 
  geom_hline(yintercept = 0,linetype=2) 
  theme_classic()

这样就解决了传统火山图 上限一行点的问题,然后可以根据ggplot2的一些参数进行细节的调整和美化。

此外单细胞还有一个特性就是有很多celltype ,那需要像NC文章那样每种cluster绘制一张图吗?是否有方式可以同时展示所有cluster的差异?

3 所有cluster热图

如文献A Spatiotemporal Organ-Wide Gene Expression and Cell Atlas of the Developing Human HeartsciPainter中的Fig2-D所示 ,提出了同时展示所有cluster差异基因的方式

这里介绍一种简单实现的方式,推荐老俊俊封装的scRNAtoolVis-R包,可以一行代码简单出图 jjVolcano 一行代码绘制单细胞火山图,可以去大佬的github star一下哈

支持和鼓励:

代码语言:javascript复制
#install.packages('devtools')
#devtools::install_github('junjunlab/scRNAtoolVis')
library(scRNAtoolVis)
jjVolcano(diffData = sce_marker,
          log2FC.cutoff = 0.25, 
          size  = 3.5, #设置点的大小
          fontface = 'italic', #设置字体形式
          #aesCol = c('purple','orange'), #设置点的颜色
          tile.col = colour, #设置cluster的颜色
          #col.type = "adjustP", #设置矫正方式
          topGeneN = 5 #设置展示topN的基因
         )

还可以通过polar函数调整形状

代码语言:javascript复制
jjVolcano(diffData = sce_marker,
          tile.col = colour,
          size  = 3.5,
          topGeneN = 5,
          fontface = 'italic',
          polar = T)  
  ylim(-8,10)

更多参数可以去github或者使用??jjVolcano查看 。

◆ ◆ ◆ ◆ ◆

精心整理(含图PLUS版)|R语言生信分析,可视化(R统计,ggplot2绘图,生信图形可视化汇总)

0 人点赞