绘制火山图&热图

上一篇文章中,我们已经对基因进行了差异分析,接下来我们根据结果中的FDR值和FC值筛选出上调基因和下调基因(上调基因:基因转录成mRNA时受到正向调控,促进表达;下调基因:转录成mRNA时受到抑制,表达量减少),并绘制成火山图与热图。

所用工具:R语言;
所需要包:ggplot2、pheatmap。

第一部分:火山图

首先,加载所需的包并导入数据:

library(ggplot2)
diff_stat <- read.csv("F:/公众号/图文素材/绘制火山图&热图/data.csv", header = TRUE, row.names = 1)

在这里插入图片描述
其次,筛选上调趋势数据和下调趋势数据,对于Fold Change值和p值阈值的选择,还需在实际的分析中视情况而定,本文以|log2FC| ≥2以及FDR p-value < 0.05作为差异OTUs的判断依据:

diff_stat[which(diff_stat$FDR < 0.05 & diff_stat$logFC >= 2),'diff'] <- 'up' #上调趋势筛选
diff_stat[which(diff_stat$FDR < 0.05 & diff_stat$logFC <= -2),'diff'] <- 'dowm' #下调趋势筛选
diff_stat[!(diff_stat$diff %in% c('up', 'dowm')),'diff'] <- 'no'

最后,我们根据判断依据,将OTUs划分为“富集”(up)、“下降”(down)以及“无差异”(no)三种水平。然后,在作图时根据预先划分的OTUs差异水平对点分别着色。火山图实质上就是一种散点图,ggplot2作为一个非常好用的作图R包,我们直接用ggplot2进行绘制:

p1 <- ggplot(diff_stat, aes(x = logFC, y = -log10(FDR))) +
  geom_point(aes(color = diff), size = 0.5) +
  scale_colour_manual(limits = c('up', 'dowm', 'no'), values = c('blue', 'red', 'gray40'), labels = c('Enriched OTUs', 'Depleted OTUs', 'No diff OTUs')) +
  labs(x = 'log2 Fold Change', y = '-log10 FDR p-value')

在这里插入图片描述
我们可以对图进行美化,修改背景颜色、添加分界线、调整标签位置:

p1 <- p1 +
  theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent')) +
  geom_vline(xintercept = c(-2, 2), color = 'gray', linetype = 2, size = 0.5) + 
  geom_hline(yintercept = -log10(0.05), color = 'gray', linetype = 2, size = 0.5) +
  theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent'), legend.position = c(0.2, 0.9))

在这里插入图片描述

知识笔记

差异分析是一个典型的多重假设检验过程,对于多重假设检验,单次检验中差异显著基因的假阳性率(p-value较小)可能会较大,而q-value和FDR值较常见的BH校正方法得到的FDR值而言,改进了其对假阳性估计的保守性。
即q-value相比于p-value更加严格,当差异基因结果较少时,可以退而求其次看p-value。Fold ChangeFold Change表示实验组比上对照组的差异表达倍数,一般表达相差2倍以上是有意义的,放宽要求1.5倍或者1.2倍也可以接受。

第二部分:热图

首先,加载所需的包并导入数据:

library(pheatmap)
sign.gene <- read.csv("F:/公众号/图文素材/绘制火山图&热图/data.csv", header = T , row.names = 1)

其次,筛选数据:

sign.gene.FDR <- sign.gene$FDR < 0.05
sign.gene.fc <- abs(sign.gene$logFC) > 2
sign.gene.all <- sign.gene.FDR & sign.gene.fc
sign.gene.real <- sign.gene[sign.gene.all, ]

如果样本中存在缺失值(例如:NA),我们可以用na.omit()进行删除:

sign.gene.real<-na.omit(sign.gene.real)

最后,绘制热图,并用基因标签代替基因ID作为热图的行标签:

pheatmap(log2(sign.gene.real[,3:85]+1), labels_row = sign.gene$Symbol)

在这里插入图片描述

也可以根据各自的需求进行美化:

pheatmap(log2(sign.gene.real[,3:85]+1),
         labels_row = sign.gene$Symbol,
         main="Heatmap",
         color = colorRampPalette(c("blue","white","red"))(256)
)

在这里插入图片描述
说明:
Color参数中的256是指色阶值,也可以理解为色阶分辨率,数值越大,热图上颜色越丰富,一般设置为256。

知识笔记

热图又称为聚类图,可以衡量样本或基因之间表达的相似性。
如本文所示的热图中,横坐标代表样本聚类,一列代表一个样本,聚类基于样本间基因表达的相似性,样本间基因表达越接近,靠的越近,以此类推。
纵坐标代表基因聚类,一行代表一个基因,聚类基于基因在样本中表达的相似性,基因在样本中表达越接近,靠的越近,以此类推。
色阶代表基因表达丰度,越红代表上调得越明显,越蓝代表下调得越明显。

总结

从热图上可以看出,所筛选出的差异基因并没有很好的区分出突变组和未突变组,所以在基因的筛选上,或者差异分析模型的选择上,需进行进一步的调整。

参考文章

https://blog.csdn.net/u012325865/article/details/87344725
http://blog.sciencenet.cn/blog-3406804-1188483.html
https://www.jianshu.com/p/f9040ca31f46

获取代码

以下是我的个人公众号,该篇的数据及代码可在公众号中回复“火山图&热图”即可获得,谢谢大家支持。
在这里插入图片描述

  • 9
    点赞
  • 74
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值