凭什么你分成两组就应该有全局差异

在生信技能树的教程:《你确定你的差异基因找对了吗?》, 提到过,必须要对你的转录水平的全局表达矩阵做好质量控制,最好是看到标准3张图

  • 左边的热图,说明我们实验的两个分组,normal和npc的很多基因表达量是有明显差异的
  • 中间的PCA图,说明我们的normal和npc两个分组非常明显的差异
  • 右边的层次聚类也是如此,说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来,我们是警告大家如果强行进行后续差异分析是有风险的。但并不意味着这样就没办法进行后续分析,我在教程:PCA都分不开的两个组强行找差异是为何提到过无数的这样的例子!

而且呢,本来你自己实验设计好分成两组就不一定是有全局差异,比如文章:ATAC-seq reveals alterations in open chromatin in pancreatic islets from subjects with type 2 diabetes. Sci Rep 2019 May 23;9(1):7785. PMID: 31123324,其配套数据集在:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129383,实验材料是Human islets of 9 non-diabetic donors and 6 donors diagnosed with T2D,也就是说是15个样品的ATAC-seq数据。

首先看15个样品的ATAC-seq数据全局差异,如下所示 :

 

如果采用我们转录组授课提到的3张图标准,可以看到这个组内差异和组间差异其实是混淆的。糖尿病患者和正常人勉强还是可以区分开来,但是混杂因素有点多,而且可以肯定这个混杂因素并不是性别,有可能是纳入的病人的年龄或者其它基础疾病,总之呢仅仅是一个糖尿病无法彻底把样品分成两组!

但是这并不影响研究者进行后续差异分析; (b) Volkano plot of human islet ATAC-seq data of donors with type 2 diabetes versus non-diabetic donors analysed with R Diffbind and edgeR packages

  • Interestingly, we found 1,078 differential ATAC-seq peaks in T2D versus control islets.

可以看到,符合要求的在糖尿病患者里面跟正常人相比有统计学显著的差异的ATAC-seq peaks仍然是有一千多个!

火山图有点丑:

丑爆了的火山图

作者在糖尿病的病人组和正常对照组里面各自挑选了一个信号值改变的peak进行IGV的载入bw文件的可视化,如下所示:

信号值上下调的IGV可视化

可以看到,这个时候研究者很鸡贼哦,他们并没有去可视化这些peaks在全部的15个样品的ATAC-seq数据的差异,而是选择了4个糖尿病患者和4个正常对照,确实就很容易看到信号值的上下调!

如果是转录组数据差异分析,基本上看我六年前的表达芯片的公共数据库挖掘系列推文即可;

但这个文章做的是 ATAC-seq数据,差异分析首先就不一样,拿到了1,078 ATAC-seq peaks 虽然说可以对应到具体的基因,然后进行生物学功能数据库注释,但是比较有特色的分析应该是peaks的注释。

比如研究者,就针对在《糖尿病的病人组》里面上调的1,078 ATAC-seq peaks 进行HOMER的注释:

 

上游分析这里我略过了,感兴趣可以去看教程ATAC-seq项目的标准分析仅收费1600, 有一个2020综述《From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis》值得看:

  • MACS2 进行 Peak calleing
  • csaw 进行差异 Peak 分析
  • MEME suite 进行 motif 检测和富集
  • ChIPseeker 进行注释和可视化
  • HMMRATAC 进行核小体检测
  • HINT-ATAC 进行足迹分析

文末友情推荐

(0)

相关推荐