前言
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
第一步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 +
格式看上去就是正规的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
这一步输入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
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
接着来看第二种情况怎么处理，与情况一大体相似，具体代码如下：
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
最后看一下多组之间的比较，据说比直接使用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
三种比较方法用起来都挺方便的，其实该软件还支持一些其他情况如Combining Replicates and Using Local Regression，这里就不介绍了，感兴趣的可以自己去查阅文献和资料。不过，有一点还是要提醒一下，MAnorm2不是直接比较peak区域的差异，而是先把peak划分为bin，划分bin大小视情况而定，默认组蛋白推荐2000转录因子1000，然后再进行比较，所以最后得到的是差异的bin区域。
最后
ChIP-seq的差异分析，目前并没有统一的标准，MAnorm2在MAnorm的核心假设基础上，通过重构其信号强度变换体系，发展了以参照样本为基准的多样本并行ChIP-seq信号标准化流程，考虑了不同样本间的系统误差，从而提高对组间差异信号的灵敏度，使其展现处更优越的性能。对于ChIP-seq差异分析大家有什么看法？
