前言
??ChIP-seq 技術是研究蛋白質(zhì)與DNA相互作用的重要手段雁歌，被廣泛用于揭示轉(zhuǎn)錄因子和組蛋白修飾在全基因組范圍內(nèi)的結合位點的分布情況台妆。近年來，隨著實驗技術的發(fā)展和測序成本的不斷降低钳降，在ChIP-seq樣本組之間進行比較分析越來越常見陷虎。隨著樣本的增加菩掏，實驗人員可以通過增加多個生物學重復控制個體差異造成的影響來提高實驗結果的可信度魂角。然而，由于ChIP-seq實驗固有的高復雜度和高噪聲水平智绸，以及不同比較場景所特有的技術困難野揪，現(xiàn)階段對多樣本ChIP-seq數(shù)據(jù)進行分組定量比較仍然是一個巨大的計算方法學挑戰(zhàn)。
??邵振課題組的方法學論文“MAnorm2 for quantitatively comparing groups of ChIP-seq samples”瞧栗，報道了其開發(fā)的新一代MAnorm2計算模型斯稳。該模型沿用了MAnorm的核心假設，通過重構其信號強度變換體系迹恐，新發(fā)展了以參照樣本為基準的多樣本并行ChIP-seq信號標準化流程挣惰。MAnorm2設計了一個經(jīng)驗貝葉斯框架，利用擬合均值-方差曲線來給單個區(qū)域的組內(nèi)變化水平賦予一個先驗分布殴边，并進一步通過平衡先驗和后驗觀測來更準確地估計ChIP-seq信號的組內(nèi)變化水平憎茂，從而提高對組間差異ChIP-seq信號的靈敏度。對該軟件的統(tǒng)計實現(xiàn)過程感興趣的話锤岸，可以詳細閱讀參考文獻（PS：本人沒有搞明白統(tǒng)計過程...@_@）竖幔。
??該模型的應用場景和統(tǒng)計模型具有良好的可擴展性，可用于比較兩組之間的比較是偷，有無重復均可拳氢；還可以用于同時比較任意多個樣本組，并發(fā)現(xiàn)其使用效果優(yōu)于傳統(tǒng)的ANOVA方法晓猛。說了那么多廢話饿幅，下面來看一下具體使用過程凡辱。
分析流程
??MAnorm2做差異分析時會將所有樣本的peak匯總然后劃分成無冗余的bin區(qū)間戒职，然后考慮每個樣本在bin區(qū)間的read和是否屬于peak開放區(qū)，最后根據(jù)統(tǒng)計模型來做差異比較透乾。在使用該軟件之前洪燥，需要用peak、bam文件來準備特定的輸入文件乳乌。對此捧韵，作者很貼心地準備好了轉(zhuǎn)換所用的腳本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文件轉(zhuǎn)化為bed文件還是很方便快捷的，轉(zhuǎn)換后格式如下：
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 +
??格式看上去就是正規(guī)的bed，但仔細一想發(fā)現(xiàn)有不合理的地方芒篷，也許大家心里也想到了搜变，那就是bed前三列表示的染色體坐標位置是不可能出現(xiàn)負值的情況，這個轉(zhuǎn)換后起始坐標居然有很多負值针炉，真是unbelievable挠他！這個原因暫且不管了，接著繼續(xù)往下走：
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做差異比較分為三種情況，一是兩組間比較且兩組都有重復樣本镰烧；二是兩組間比較但組內(nèi)沒有重復樣本拢军；三是多組間比較。下面先來看看第一種情況:
library(MAnorm2)
# 自帶的測試數(shù)據(jù)
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)
# 統(tǒng)計檢驗差異區(qū)域
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為自帶測試數(shù)據(jù)
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
最后看一下多組之間的比較，據(jù)說比直接使用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區(qū)域的差異鸿竖，而是先把peak劃分為bin，劃分bin大小視情況而定铸敏，默認組蛋白推薦2000轉(zhuǎn)錄因子1000缚忧，然后再進行比較，所以最后得到的是差異的bin區(qū)域杈笔。
最后
??ChIP-seq的差異分析闪水，目前并沒有統(tǒng)一的標準，MAnorm2在MAnorm的核心假設基礎上蒙具，通過重構其信號強度變換體系球榆，發(fā)展了以參照樣本為基準的多樣本并行ChIP-seq信號標準化流程，考慮了不同樣本間的系統(tǒng)誤差禁筏，從而提高對組間差異信號的靈敏度持钉，使其展現(xiàn)處更優(yōu)越的性能。對于ChIP-seq差異分析大家有什么看法篱昔？
