前一篇推送已经与大家分享了近期吴柘团队发表在NG上的拟南芥雄配子发生过程染色质特征的文章。
这篇文章涉及的数据有ChIP-seq、RNA-seq、ATAC-seq，分析上整体不是特别复杂，但是对于H3K4me3和H3K27me3两种组蛋白修饰的ChIP-seq分析还是比较经典，而且作者在github上(https://github.com/Blairewen/Arabidopsis-Pollen )提供了详细的代码，所以这次推送我们一起学习作者用到的ChIP-seq分析流程。
一 数据下载
Fig 1
Fig2
示例数据用Spm的H3K4me3和H3K27me3数据（https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-10966），各三个生物学重复（Fig 1 ）。
文章的数据上传到了BioStudies数据库 （Fig 2），我是第一次用这个数据库。我不太清楚这个数据库数据下载的方式。我在这个文件中找到了下载链接，然后用命令行的wget下载。
二 数据的过滤与质控
Fig 3
拿到测序数据我们都会先做这一步，做法有很多种。
作者提供的流程中是用fastp对ChIP-seq数据进行过滤和质控。
Fig 3 中-w 指定运行的线程数， -l 默认值是15，长度小于该值的测序读段被去除，--detect_adapter_for_pe 默认情况下会自动检测单端测序数据的接头，加入此参数以适用于双端数据。
三 比对及比对后过滤
Fig 4
Fig 5
对于DNA的比对，我们常用的软件是bowtie/bowtie2/bwa，但是这篇文章作者用到的是hisat2，为了保持一致，也延续了作者的代码。
Fig 4中 -p 指定线程数；
-x指定参考基因组的index；
--no-temp-splicesite 禁用使用找到的剪接位点；
--summary-file 后加文件名，将比对报告输入到该文件；
Samtools view –q 30 只保留比对质量大于30的比对记录。
用picard去除PCR重复 （Fig 5），注意REMOVE_DUPLICATES一定要写为true，否则只标记，不直接去除。
四 Peak Calling
Fig 6
Fig 7
对组蛋白ChIP-seq数据检测峰值时为了找到组蛋白修饰的富集区，根据组蛋白结合特性不同，峰可以分成宽峰(broad peak)和窄峰 (narrow peak)。
H3K9me2、H3K9me3、H3K27me3、H3K36me3、H2AK121ub和H2AK119ub1等修饰常用宽峰，H3K4me3常用窄峰。
Call peak的软件最为经典的就是MACS2，该软件默认是检验窄峰 （Fig 6），所以检验宽峰时需要加入-broad参数（Fig 7）。
--SPMR参数用以做测序深度的标准化。
Fig 8
Fig 9
由于作者时对三个生物学重复分别进行peak calling，所以下一步需要对不同生物学重复的peak进行合并。
对于H3K4me3可直接合并（Fig 8），而H3K27me3经常是一个较宽的区间，所以先合并了相距5kb以内的峰，然后对三个生物学重复进行合并（Fig 9）。
五 Peak交集的显著性
如果我们获得了两组基因id，想判断有无交集，以及交集是否显著。很简单，对两组基因画一个韦恩图然后做一个超几何分布检验或者费歇尔精确检验即可。但是如果对于两组区间，似乎麻烦了些，我之前一直没有找到合适的方法检验两者的交集是否显著。这篇文章作者也面临这个问题，并且提出了一个比较简单可行的方法。
Fig 10
Fig 11
还是代入目前的例子，我们想知道在Spm中的H3K27me3峰和H3K4me3峰是否存在显著的重叠，通过韦恩图可以看到，两者共有7843个峰，H3K27me3特有2918个峰，H3K4me3特有2998个峰，而H3K27me3共有2918+7843=10761个峰，H3K4me3共有2998+7843=10841个峰 （Fig 10）。
作者对H3K27me3的10761个峰和H3K4me3的10841个峰分别生成了相同染色体上的随机峰 （Fig 11）。然后统计生成的两组随机峰的交集是否超过7843，这样的计算重复一万次，统计这一万次重复中交集超过7843的次数，然后用这个次数除以10000以表征P值。理论上这个次数越大最后结果越接近1，说明两组的交集富集越不显著，反之越显著。
本文使用 文章同步助手 同步
网友评论