Capture Hi-C 是一種檢測(cè)至少在一端涉及感興趣的 DNA 區(qū)域，如基因啟動(dòng)子染色體相互作用的有效方法。 Chicdiff，是一個(gè)用于捕獲 Hi-C 數(shù)據(jù)中差分相互作用的魯棒檢測(cè)的 R 軟件包。Chicdiff 增強(qiáng)了計(jì)數(shù)數(shù)據(jù)的最先進(jìn)的差異檢驗(yàn)方法，具有定制的標(biāo)準(zhǔn)化和多種檢驗(yàn)程序，可以解釋 Capture Hi-C 的特定統(tǒng)計(jì)特性。作者用Chicdiff 在人單核細(xì)胞和 CD4 + T 細(xì)胞中發(fā)表的promoter-Capture Hi-C 數(shù)據(jù)，鑒定了多種細(xì)胞類(lèi)型特異性相互作用，并證實(shí)了啟動(dòng)子相互作用和基因表達(dá)之間的總體正相關(guān)性。
install.packages("devtools")
library(devtools)
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="Chicdiff")
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="ChicdiffData", force=T)
Chicdiff 是一種在捕獲 Hi-C 數(shù)據(jù)中檢測(cè)統(tǒng)計(jì)顯著差異相互作用的方法。這個(gè)文檔將引導(dǎo)你通過(guò)一個(gè)典型的Chicdiff analysis.。
簡(jiǎn)而言之，Chicdiff 將 DESeq2中實(shí)現(xiàn)的計(jì)數(shù)數(shù)據(jù)的適度差異檢驗(yàn)與用于信號(hào)歸一化和 p 值加權(quán)的 CHi-C 特定程序相結(jié)合。為了提高威力，Chicdiff 還將每個(gè)誘餌的每個(gè)顯著相互作用區(qū)域周?chē)膸讉€(gè)片段的讀數(shù)結(jié)合起來(lái)。更具體地說(shuō)，Chicdiff 使用來(lái)自Chicago的背景估計(jì)來(lái)構(gòu)建用于 DESeq2差異檢驗(yàn)的定制縮放矩陣，來(lái)自 DESeq2的 Wald test檢驗(yàn)產(chǎn)生 p 值提交給加權(quán)過(guò)程，如下所述。
由于 CHi-C 數(shù)據(jù)中的信噪比和效應(yīng)大小與距離有很強(qiáng)的相關(guān)性，Chicdiff 進(jìn)行非均勻的多重檢驗(yàn)校正，因此較長(zhǎng)距離的較弱信號(hào)得到更有力的校正。它使用 IHW 軟件包以交互距離作為協(xié)變量來(lái)學(xué)習(xí) p 值上的權(quán)重，使得被拒絕的零假設(shè)的總分?jǐn)?shù)最大化。由于提交用于 Chicdiff 測(cè)試的一組相互作用通常由于選擇偏差而具有不均勻的 Chicdiff p 值分布，因此從整個(gè) Capture Hi-C 數(shù)據(jù)集采樣的“訓(xùn)練集”的相互作用被用于權(quán)重訓(xùn)練。通過(guò)這種方式學(xué)到的權(quán)重然后被應(yīng)用到檢驗(yàn)集中的 p 值，并且得到的加權(quán) p 值被調(diào)整用于多個(gè)檢驗(yàn)并報(bào)告給用戶(hù)。
來(lái)自 ChicaffData 的示例數(shù)據(jù)集使用約0.6 Gb RAM，在標(biāo)準(zhǔn)筆記本電腦上處理需要幾分鐘。在每種條件下，兩個(gè)全基因組 CHi-C 數(shù)據(jù)的生物學(xué)重復(fù)的典型工作需要約30-60分鐘，并使用約10-15 Gb RAM。具有更多重復(fù)的高覆蓋率數(shù)據(jù)集將需要更長(zhǎng)的時(shí)間來(lái)處理和使用更多 RAM。
該軟件包只包含來(lái)自?xún)蓚€(gè)重復(fù)的單核細(xì)胞和初始 CD4 + T 細(xì)胞的19號(hào)染色體數(shù)據(jù)(相比之下，文中分別公布了三個(gè)和四個(gè)重復(fù))
library(Chicdiff)
library(ChicdiffData)
輸入文件
1.Two restriction map information files (“design files”)
- Restriction map file (.rmap)
限制映射文件(.Rmap)-包含限制性片段坐標(biāo)的底層文件。默認(rèn)情況下，有4列: chr、 start、 end、 FragmentID - Bait map file (.baitmap)
誘餌映射文件(。Baitmap)-包含帶誘餌的限制性片段的坐標(biāo)及其相關(guān)注釋的底層文件。默認(rèn)情況下，有5列: chr、 start、 end、 fragmentID、 baitAnnote。文件中指定的區(qū)域(包括它們的片段 ID)必須是。Rmap 文件。BaitAnnote 是一個(gè)文本字段，僅用于對(duì)輸出和繪圖進(jìn)行注釋
dataPath <- system.file("extdata", package="ChicdiffData")
testDesignDir <- file.path(dataPath, "designDir")
dir(testDesignDir)
# [1] "chr19_GRCh37_HindIII.baitmap" "chr19_GRCh37_HindIII.rmap"
- One or more peakfile(s) defining interactions of interest.
chicagoTools/makePeakMatrix.R
peakFiles <- file.path(dataPath,"peakMatrix_cd4_mono_unmerged.txt") ##makePeakMatrix.R生成
makePeakMatrix.R 腳本需要兩個(gè)位置參數(shù)：
<names-file>：這是一個(gè)制表符分隔的文件，包含樣本名稱(chēng)（第一列）和指向輸入 Rds 文件的完整路徑（第二列）。這些輸入 Rds 文件應(yīng)該是 CHiCAGO 分析的輸出結(jié)果。<output-prefix>：這是輸出文件名的前綴，可以包含文件夾路徑。
此外，腳本還接受許多可選參數(shù)，例如 --cutoff（可能用于設(shè)置峰值閾值）和 --scorecol（可能用于指定在輸入文件中表示分?jǐn)?shù)或強(qiáng)度的列）。
一個(gè)使用這個(gè)腳本的命令可能看起來(lái)像這樣：
Rscript makePeakMatrix.R --cutoff 5 names_and_paths.tsv output_prefix
在這個(gè)命令中，names_and_paths.tsv 是一個(gè)包含樣本名稱(chēng)和對(duì)應(yīng) Rds 文件路徑的文件，output_prefix 是輸出文件名的前綴。--cutoff 5 是一個(gè)可選參數(shù)，設(shè)置峰值閾值為 5。
- Count data files in Chicago input format
.chinput 文件，chicagoTools/bam2chicago.sh轉(zhuǎn)換bam文件而成，
testDataPath_CD4 <- file.path(dataPath, "CD4")
testDataPath_Mono <- file.path(dataPath, "monocytes")
dir(testDataPath_CD4)
dir(testDataPath_Mono)
countData <- list(
CD4 = c(NCD4_22 = file.path(testDataPath_CD4, "unitTest_CD41.chinput"),
NCD4_23 = file.path(testDataPath_CD4, "unitTest_CD42.chinput")
),
Mono = c(Mon_2 = file.path(testDataPath_Mono, "unitTest_Mono2.chinput"),
Mon_3 = file.path(testDataPath_Mono, "unitTest_Mono3.chinput")
)
)
4 Chicago output files for each separate replicate (saved as either Rds or .Rda)
testDataPath_rda <- system.file("data", package="ChicdiffData")
dir(testDataPath_rda)
## [1] "unitTest_CD41.RDa" "unitTest_CD42.RDa" "unitTest_Mono2.RDa"
## [4] "unitTest_Mono3.RDa"
chicagoData <- list(
CD4 = c(NCD4_22 = file.path(testDataPath_rda, "unitTest_CD41.RDa"),
NCD4_23 = file.path(testDataPath_rda, "unitTest_CD42.RDa")
),
Mono = c(Mon_2 = file.path(testDataPath_rda, "unitTest_Mono2.RDa"),
Mon_3 = file.path(testDataPath_rda, "unitTest_Mono3.RDa")
)
如果Peakfile(s) 文件是在示例級(jí)別生成的(與每次重復(fù)相反) ，則不應(yīng)該命名每個(gè)向量中的元素
主要流程
我們?cè)跍y(cè)試數(shù)據(jù)上運(yùn)行 Chicdiff，首先，我們使用 setChicdiffExperiment()設(shè)置 Chicdiff 實(shí)驗(yàn)設(shè)置:
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test")
注意，除了返回設(shè)置列表之外，該函數(shù)還會(huì)將其保存為 Rds 文件 * * * _ sets?？芍貜?fù)使用的 Rds。
可選: 您可以直接在命令行中修改設(shè)置(參見(jiàn)？默認(rèn)設(shè)置)。例如，我們可以將運(yùn)行模式切換到并行模式(這會(huì)加快運(yùn)行時(shí)速度，但增加內(nèi)存需求) :
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test", settings = list(parallel=TRUE))
或者，可以在設(shè)置文件中提供自定義參數(shù):
settingsFile = file.path(dataPath,"SettingsFile.txt")
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test", settingsFile = settingsFile)
最后，我們使用chicdiffPipeline()運(yùn)行流程:
output <- chicdiffPipeline(chicdiff.settings)
流水線(xiàn)將基于 RUtended 選項(xiàng)定義擴(kuò)展區(qū)域(默認(rèn)情況下，從軟化相互作用的另一端的每個(gè)相互作用峰每個(gè)方向5個(gè)片段) ，執(zhí)行標(biāo)準(zhǔn)化，差異檢驗(yàn)，基于“測(cè)試集”的權(quán)重估計(jì)和 p 值加權(quán)和多重測(cè)試校正。
對(duì)于每個(gè)bait-region 相互作用，輸出數(shù)據(jù)表列出標(biāo)準(zhǔn)化相互作用讀取計(jì)數(shù)的 log2FoldChange，以及原始的，加權(quán)的和加權(quán)調(diào)整的 p 值(分別為 pvalue，weighted_pvalue 和weighted_padj)。
head(output)
## group baseMean log2FoldChange lfcSE stat pvalue padj
## 1: 1 98.04145 0.4654394 0.3430257 1.3568645 0.17482427 0.4287443
## 2: 1 99.90272 0.5800085 0.3265890 1.7759582 0.07573981 0.2708560
## 3: 1 58.45174 0.4928256 0.4360129 1.1303006 0.25834959 0.5231848
## 4: 1 61.89843 0.5010415 0.4122627 1.2153453 0.22423442 0.4872660
## 5: 1 60.04236 0.5523034 0.4108116 1.3444201 0.17881258 0.4340686
## 6: 1 35.74176 0.3311124 0.4737406 0.6989319 0.48459458 0.7169037
## baitID maxOE minOE regionID OEchr OEstart OEend baitchr baitstart
## 1: 320526 320524 320517 100 19 436227 524426 19 530387
## 2: 320526 320536 320528 101 19 542386 622492 19 530387
## 3: 320526 320544 320534 102 19 594189 749359 19 530387
## 4: 320526 320545 320535 103 19 596117 753972 19 530387
## 5: 320526 320546 320536 104 19 598270 763809 19 530387
## 6: 320526 320550 320540 105 19 638626 864727 19 530387
## baitend avDist uniform shuff avgLogDist avWeights weight
## 1: 539467 -53060.62 0.009665534 0.05048031 10.87919 7.148253 2.8532
## 2: 539467 45032.11 0.430116629 0.44131718 10.71513 7.148253 2.8532
## 3: 539467 111436.45 0.092592717 0.06689947 11.62121 7.148253 2.8532
## 4: 539467 125665.00 0.303875251 0.02621149 11.74137 7.148253 2.8532
## 5: 539467 140364.82 0.067529660 0.19320485 11.85200 7.148253 2.8532
## 6: 539467 214648.36 0.163497301 0.10922867 12.27676 7.148253 2.8532
## weighted_pvalue weighted_padj
## 1: 0.06127305 0.2951175
## 2: 0.02654557 0.1621728
## 3: 0.09054731 0.3850313
## 4: 0.07859050 0.3491772
## 5: 0.06267089 0.2997665
## 6: 0.16984248 0.5832029
除了返回輸出數(shù)據(jù)表之外，chicdiffPipeline還將其保存為 Rds文件如 _ result.Rds.它還將為 Rds 文件 * _ count put 中的每個(gè)交互保存最終的 count 表以及_countput.Rds. .可以使用 saveAuxFiles = TRUE 設(shè)置生成更多的輸出文件-請(qǐng)參見(jiàn)？詳細(xì)信息請(qǐng)參閱?chicdiffpipeline。
2.3 Output diagnostic plots
chicdiffPipeline()生成了幾個(gè)保存在工作目錄中的圖。這些為 chr19數(shù)據(jù)集生成的圖是作為 ChicffData 包的一部分便存有的。請(qǐng)注意，對(duì)于全基因組數(shù)據(jù)，下面描述的趨勢(shì)應(yīng)該更加明顯。
resultsPath <- file.path(dataPath, "CD4_Mono_results")
IHWdecisionBoundaryPlot <- png::readPNG(file.path(resultsPath,"test_IHWdecisionBoundaryPlot.png"))
grid::grid.raster(IHWdecisionBoundaryPlot)
這張圖顯示了作為協(xié)變量函數(shù)的未加權(quán) p 值的隱含決策邊界。邊界對(duì)協(xié)變量(距離)的依賴(lài)性表明，這個(gè)協(xié)變量是信息量大的。低 p 值的趨勢(shì)是在低距離富集。
Estimated weights
IHWweightPlot <- png::readPNG(file.path(resultsPath,"test_IHWweightPlot.png"))
grid::grid.raster(IHWweightPlot)
從圖中可以看出，權(quán)重與距離(stratum)有關(guān)。正如預(yù)期的那樣，在數(shù)據(jù)的不同隨機(jī)子集(folds)上計(jì)算的權(quán)重函數(shù)的行為是相似的。對(duì)于這個(gè) chr19數(shù)據(jù)，跨越較長(zhǎng)距離的交互會(huì)受到懲罰，賦予以較低的權(quán)重，而跨越較短距離的交互則被優(yōu)先考慮。
差異的分析的可視化
進(jìn)行兩種比較條件的可視化，可以使用 plotdiffBaits()函數(shù)繪制相互作用的原始count與它們與相應(yīng)bait碎片的線(xiàn)性距離的圖像。對(duì)應(yīng)于重要相互作用的其他端的擴(kuò)展區(qū)域表示為根據(jù)其調(diào)整后的加權(quán)值進(jìn)行顏色編碼的區(qū)間。
Chicdiff 管道自動(dòng)調(diào)用 plotdefBaits ()為從前100個(gè)包含與最低加權(quán) p 值相互作用的誘餌中選擇的4個(gè)隨機(jī)誘餌生成配置文件。顯示了來(lái)自誘餌的1Mb 內(nèi)的相互作用，并且根據(jù)分別在每個(gè)條件中檢測(cè)到的相應(yīng)相互作用的 Chicago score(score>5: red; 3)對(duì)圖上的數(shù)據(jù)點(diǎn)進(jìn)行顏色編碼。
To use plotDiffBaits() outside of the pipeline to plot the profiles of baits of interest, it needs to be provided with the output data table (saved in test_results.Rds by chicdiffPipeline, assuming that outprefix="test" by default), count table (saved in test_countput.Rds), baitmap file (from Chicago design directory) and a vector of baitIDs of interest (see ?plotDiffBaits for the description of additional parameters).
outputRds <- file.path(resultsPath, "test_results.Rds")
countputRds <- file.path(resultsPath, "test_countput.Rds")
bmapRds <- file.path(testDesignDir, "chr19_GRCh37_HindIII.baitmap")
baits <- c(327182, 330614,330598, 327700)
plotDiffBaits(outputRds, countputRds, bmapRds, baits)
在差異相互作用區(qū)域內(nèi)對(duì)單個(gè)片段進(jìn)行優(yōu)先排序
除非 RUexpad 設(shè)置設(shè)置為零，否則 Chicdiff 工作在誘餌和池區(qū)域之間的交互級(jí)別，包含多個(gè)限制性片段。
為了在單個(gè)片段獲得差異信號(hào)的指示，Chicdiff 提供了函數(shù) getCandidateInteractions ()。對(duì)于落在多個(gè)匯集區(qū)域內(nèi)的碎片，該函數(shù)通過(guò)取最小值(默認(rèn)值)或更正式地使用調(diào)和平均值方法(由方法參數(shù)確定)來(lái)組合它們的微分 p 值，用于在包調(diào)和平均值中實(shí)現(xiàn)的相關(guān)檢驗(yàn)。還可以指定用于篩選輸出的最大 p 值截止值(pvcut)。
為了優(yōu)先考慮差異相互作用區(qū)域內(nèi)推定的“驅(qū)動(dòng)器”相互作用，getCandidateInteractions 允許在條件之間通過(guò)最小 | asinh (Chicago 得分) | (minDeltaAsinhScore)過(guò)濾它們。Asinh 轉(zhuǎn)換有助于壓縮得分的上限，在這個(gè)范圍內(nèi)，他們之間的差異比那些低范圍的差異更難解釋。
outCI <- getCandidateInteractions(chicdiff.settings = chicdiff.settings,
output = output, peakFiles = peakFiles,
pvcut = 0.05, minDeltaAsinhScore = 1)
head(outCI)
