口腔医学小站

口腔文献阅读与资源分享

MAnorm2:ChIP-seq数据分组比较软件 – 简书

| 暂无评论

MAnorm2:ChIP-seq数据分组比较软件

32021.05.12 21:06:10字数 1,345阅读 1,880

前言

  ChIP-seq 技术是研究蛋白质与DNA相互作用的重要手段,被广泛用于揭示转录因子和组蛋白修饰在全基因组范围内的结合位点的分布情况。近年来,随着实验技术的发展和测序成本的不断降低,在ChIP-seq样本组之间进行比较分析越来越常见。随着样本的增加,实验人员可以通过增加多个生物学重复控制个体差异造成的影响来提高实验结果的可信度。然而,由于ChIP-seq实验固有的高复杂度和高噪声水平,以及不同比较场景所特有的技术困难,现阶段对多样本ChIP-seq数据进行分组定量比较仍然是一个巨大的计算方法学挑战。
  邵振课题组的方法学论文“MAnorm2 for quantitatively comparing groups of ChIP-seq samples”,报道了其开发的新一代MAnorm2计算模型。该模型沿用了MAnorm的核心假设,通过重构其信号强度变换体系,新发展了以参照样本为基准的多样本并行ChIP-seq信号标准化流程。MAnorm2设计了一个经验贝叶斯框架,利用拟合均值-方差曲线来给单个区域的组内变化水平赋予一个先验分布,并进一步通过平衡先验和后验观测来更准确地估计ChIP-seq信号的组内变化水平,从而提高对组间差异ChIP-seq信号的灵敏度。对该软件的统计实现过程感兴趣的话,可以详细阅读参考文献(PS:本人没有搞明白统计过程…@_@)。

  该模型的应用场景和统计模型具有良好的可扩展性,可用于比较两组之间的比较,有无重复均可;还可以用于同时比较任意多个样本组,并发现其使用效果优于传统的ANOVA方法。说了那么多废话,下面来看一下具体使用过程。

分析流程

  MAnorm2做差异分析时会将所有样本的peak汇总然后划分成无冗余的bin区间,然后考虑每个样本在bin区间的read和是否属于peak开放区,最后根据统计模型来做差异比较。在使用该软件之前,需要用peak、bam文件来准备特定的输入文件。对此,作者很贴心地准备好了转换所用的脚本sam2bed.py、profile_bins.py,详情见github仓库MAnorm2_utils。具体代码如下:

samtools view -O SAM -o sample.sam sample.bam
python sam2bed.py -i sample.sam -o sample.bed
12

第一步bam文件转化为bed文件还是很方便快捷的,转换后格式如下:

chrI    -7      143     E00516:574:H3MHVCCX2:6:1104:2432:2715   60      +
chrI    -23     127     E00516:574:H3MHVCCX2:6:1104:26068:24638 60      +
chrI    -23     127     E00516:574:H3MHVCCX2:6:1104:27397:25605 60      +
chrI    -30     92      E00516:574:H3MHVCCX2:6:1104:15727:34641 60      +
chrI    -60     90      E00516:574:H3MHVCCX2:6:1105:27093:37489 60      +
chrI    -8      142     E00516:574:H3MHVCCX2:6:1106:10805:13281 60      +
chrI    -8      142     E00516:574:H3MHVCCX2:6:1106:10886:13351 60      +
chrI    -20     130     E00516:574:H3MHVCCX2:6:1107:28595:70785 60      +
chrI    -49     101     E00516:574:H3MHVCCX2:6:1108:23734:49531 60      +
chrI    0       150     E00516:574:H3MHVCCX2:6:1110:2077:23372  60      +
chrI    -54     96      E00516:574:H3MHVCCX2:6:1118:8450:47738  48      +
chrI    -12     138     E00516:574:H3MHVCCX2:6:1119:8532:68834  60      +
chrI    -74     76      E00516:574:H3MHVCCX2:6:1120:20060:12683 52      +
12345678910111213

  格式看上去就是正规的bed,但仔细一想发现有不合理的地方,也许大家心里也想到了,那就是bed前三列表示的染色体坐标位置是不可能出现负值的情况,这个转换后起始坐标居然有很多负值,真是unbelievable!这个原因暂且不管了,接着继续往下走:

python profile_bins.py --peaks=myc1_peaks.narrowPeak,yaf9d_peaks.narrowPeak --reads=myc1_reads.bed,yaf9d_reads.bed --labs=myc1,yaf9d -n myc1_yaf9d
1

这一步输入peak文件和上一步得到的bed文件,就会生成最终需要的输入文件,格式如下所示:

chrom   start   end     myc1.read_cnt   yaf9d.read_cnt  myc1.occupancy  yaf9d.occupancy
chrI    0       591     412     548     1       1
chrI    1808    2185    176     238     0       1
chrI    11703   12113   120     276     0       1
chrI    45896   46719   521     774     0       1
chrI    56742   57311   340     321     1       0
chrI    60586   61545   921     848     1       0
chrI    62426   63070   690     601     1       0
chrI    68394   68881   504     873     0       1
123456789

MAnorm做差异比较分为三种情况,一是两组间比较且两组都有重复样本;二是两组间比较但组内没有重复样本;三是多组间比较。下面先来看看第一种情况:

library(MAnorm2)
# 自带的测试数据
head(H3K27Ac)
chrom  start    end GM12890_H3K27Ac_1.read_cnt GM12891_H3K27Ac_1.read_cnt
1  chr1  29023  29346                          8                         16
2  chr1 712983 715873                        440                        498
3  chr1 740056 741095                         81                         54
4  chr1 751252 753001                          2                        139
5  chr1 760716 764177                        234                        329
6  chr1 800556 801985                          4                         26
  GM12891_H3K27Ac_2.read_cnt GM12892_H3K27Ac_1.read_cnt
1                          9                         23
2                        477                        439
3                         39                         44
4                         84                         11
5                        350                        326
6                         59                         16
  GM12892_H3K27Ac_2.read_cnt GM12890_H3K27Ac_1.occupancy
1                         22                           0
2                        508                           1
3                         40                           1
4                         11                           0
5                        376                           1
6                         19                           0
  GM12891_H3K27Ac_1.occupancy GM12891_H3K27Ac_2.occupancy
1                           0                           0
2                           1                           1
3                           0                           0
4                           1                           0
5                           1                           1
6                           0                           1
  GM12892_H3K27Ac_1.occupancy GM12892_H3K27Ac_2.occupancy
1                           1                           1
2                           1                           1
3                           1                           0
4                           0                           0
5                           1                           1
6                           0                           0

# 标准化
norm <- normalize(H3K27Ac, count = 5:6, occupancy = 10:11)
norm <- normalize(norm, count = 7:8, occupancy = 12:13)
conds <- list(GM12891 = bioCond(norm[5:6], norm[10:11], name = "GM12891"), GM12892 = bioCond(norm[7:8], norm[12:13], name = "GM12892"))
conds <- normBioCond(conds)

# 拟合mean-variance曲线
conds <- fitMeanVarCurve(conds, method = "parametric", occupy.only = TRUE)
# 统计检验差异区域
res <- diffTest(conds[[1]], conds[[2]])
head(res)
GM12891.mean GM12892.mean       Mval   Mval.se     Mval.t         pval
1     2.962750     4.378686  1.4159352 0.6952966  2.0364479 4.170540e-02
2     8.599608     8.842037  0.2424297 0.1687241  1.4368406 1.507633e-01
3     4.978799     5.282421  0.3036218 0.4299372  0.7062003 4.800636e-01
4     6.286417     3.357143 -2.9292743 0.4749251 -6.1678654 6.921801e-10
5     8.043483     8.401049  0.3575665 0.1852377  1.9303114 5.356827e-02
6     4.742082     4.015697 -0.7263855 0.5490536 -1.3229775 1.858429e-01
          padj
1 7.478628e-02
2 2.234701e-01
3 5.746783e-01
4 5.306019e-09
5 9.286274e-02
6 2.663536e-01
12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364

接着来看第二种情况怎么处理,与情况一大体相似,具体代码如下:

library(MAnorm2)

# H3K27Ac为自带测试数据
norm <- normalize(H3K27Ac, c(5, 7), c(10, 12))
conds <- list(GM12891 = bioCond(norm[5], norm[10], name = "GM12891"), GM12892 = bioCond(norm[7], norm[12], name = "GM12892"))

conds$blind <- bioCond(norm[c(5, 7)], norm[c(10, 12)], occupy.num = 2, name = "blind")
conds <- fitMeanVarCurve(conds, method = "parametric", occupy.only = TRUE, init.coef = c(0.1, 10))
res <- diffTest(conds[[1]], conds[[2]])
head(res)
GM12891.mean GM12892.mean       Mval   Mval.se      Mval.t       pval
1     3.423438     4.554589  1.1311505 1.3268737  0.85249295 0.44690929
2     8.622954     8.779719  0.1567654 0.5014371  0.31263231 0.77181360
3     5.246252     5.475733  0.2294818 0.8931912  0.25692353 0.81123771
4     6.680080     3.523562 -3.1565184 0.9576459 -3.29612253 0.03509689
5     7.991326     8.350939  0.3596130 0.5272771  0.68201903 0.53654681
6     4.146230     4.044394 -0.1018357 1.2841904 -0.07929955 0.94100180
       padj
1 0.8219632
2 0.9459544
3 0.9568767
4 0.5210109
5 0.8637627
6 0.9864314
123456789101112131415161718192021222324

最后看一下多组之间的比较,据说比直接使用ANOVA分析的结果要好,下面来看一下具体代码:

norm <- normalize(H3K27Ac, count = 4, occupancy = 9)
norm <- normalize(norm, count = 5:6, occupancy = 10:11)
norm <- normalize(norm, count = 7:8, occupancy = 12:13)

conds <- list(GM12890 = bioCond(norm[4], norm[9], name = "GM12890"), GM12891 = bioCond(norm[5:6], norm[10:11], name = "GM12891"), GM12892 = bioCond(norm[7:8], norm[12:13], name = "GM12892"))

conds <- normBioCond(conds)
conds <- fitMeanVarCurve(conds, method = "parametric", occupy.only = TRUE)
res <- aovBioCond(conds)
head(res)
conds.mean between.ms   within.ms  prior.var posterior.var      mod.f
1   3.877391 1.24798909 0.154992959 0.42340216    0.42340216  2.9475265
2   8.838753 0.02103485 0.003559212 0.02236859    0.02236859  0.9403745
3   5.848512 0.22555106 0.073885521 0.11474682    0.11474682  1.9656410
4   3.979706 9.02918943 0.124069261 0.39503584    0.39503584 22.8566336
5   8.233196 0.15408184 0.004609297 0.02930508    0.02930508  5.2578535
6   3.997513 2.98081654 0.272270951 0.39030110    0.39030110  7.6372230
          pval         padj
1 5.246933e-02 6.504417e-02
2 3.904816e-01 4.199157e-01
3 1.400661e-01 1.623538e-01
4 1.184378e-10 4.159661e-10
5 5.206469e-03 7.439914e-03
6 4.821655e-04 7.897353e-04
123456789101112131415161718192021222324

  三种比较方法用起来都挺方便的,其实该软件还支持一些其他情况如Combining Replicates and Using Local Regression,这里就不介绍了,感兴趣的可以自己去查阅文献和资料。不过,有一点还是要提醒一下,MAnorm2不是直接比较peak区域的差异,而是先把peak划分为bin,划分bin大小视情况而定,默认组蛋白推荐2000转录因子1000,然后再进行比较,所以最后得到的是差异的bin区域。

最后

  ChIP-seq的差异分析,目前并没有统一的标准,MAnorm2在MAnorm的核心假设基础上,通过重构其信号强度变换体系,发展了以参照样本为基准的多样本并行ChIP-seq信号标准化流程,考虑了不同样本间的系统误差,从而提高对组间差异信号的灵敏度,使其展现处更优越的性能。对于ChIP-seq差异分析大家有什么看法?

参考

多样本ChIP-seq数据分组定量比较的MAnorm2计算模型

12人点赞

更多精彩内容,就在简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

生信店小二记录和分享生信相关的内容。工欲善其事,必先利其器。『鹅厂号:生信杂记』
总资产78共写了7.7W字获得683个赞共481个粉丝

发表评论

*为必填字段!