Capture Hi-C 是一種檢測至少在一端涉及感興趣的 DNA 區(qū)域姊舵，如基因啟動子染色體相互作用的有效方法晰绎。 Chicdiff，是一個用于捕獲 Hi-C 數(shù)據(jù)中差分相互作用的魯棒檢測的 R 軟件包括丁。Chicdiff 增強了計數(shù)數(shù)據(jù)的最先進的差異檢驗方法荞下，具有定制的標準化和多種檢驗程序，可以解釋 Capture Hi-C 的特定統(tǒng)計特性躏将。作者用Chicdiff 在人單核細胞和 CD4 + T 細胞中發(fā)表的promoter-Capture Hi-C 數(shù)據(jù)锄弱，鑒定了多種細胞類型特異性相互作用，并證實了啟動子相互作用和基因表達之間的總體正相關性祸憋。
install.packages("devtools")
library(devtools)
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="Chicdiff")
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="ChicdiffData", force=T)
Chicdiff 是一種在捕獲 Hi-C 數(shù)據(jù)中檢測統(tǒng)計顯著差異相互作用的方法会宪。這個文檔將引導你通過一個典型的Chicdiff analysis.。
簡而言之蚯窥，Chicdiff 將 DESeq2中實現(xiàn)的計數(shù)數(shù)據(jù)的適度差異檢驗與用于信號歸一化和 p 值加權(quán)的 CHi-C 特定程序相結(jié)合掸鹅。為了提高威力，Chicdiff 還將每個誘餌的每個顯著相互作用區(qū)域周圍的幾個片段的讀數(shù)結(jié)合起來拦赠。更具體地說巍沙，Chicdiff 使用來自Chicago的背景估計來構(gòu)建用于 DESeq2差異檢驗的定制縮放矩陣，來自 DESeq2的 Wald test檢驗產(chǎn)生 p 值提交給加權(quán)過程荷鼠，如下所述句携。
由于 CHi-C 數(shù)據(jù)中的信噪比和效應大小與距離有很強的相關性，Chicdiff 進行非均勻的多重檢驗校正允乐，因此較長距離的較弱信號得到更有力的校正矮嫉。它使用 IHW 軟件包以交互距離作為協(xié)變量來學習 p 值上的權(quán)重，使得被拒絕的零假設的總分數(shù)最大化牍疏。由于提交用于 Chicdiff 測試的一組相互作用通常由于選擇偏差而具有不均勻的 Chicdiff p 值分布蠢笋，因此從整個 Capture Hi-C 數(shù)據(jù)集采樣的“訓練集”的相互作用被用于權(quán)重訓練。通過這種方式學到的權(quán)重然后被應用到檢驗集中的 p 值鳞陨，并且得到的加權(quán) p 值被調(diào)整用于多個檢驗并報告給用戶昨寞。
來自 ChicaffData 的示例數(shù)據(jù)集使用約0.6 Gb RAM，在標準筆記本電腦上處理需要幾分鐘厦滤。在每種條件下援岩，兩個全基因組 CHi-C 數(shù)據(jù)的生物學重復的典型工作需要約30-60分鐘，并使用約10-15 Gb RAM掏导。具有更多重復的高覆蓋率數(shù)據(jù)集將需要更長的時間來處理和使用更多 RAM享怀。
該軟件包只包含來自兩個重復的單核細胞和初始 CD4 + T 細胞的19號染色體數(shù)據(jù)(相比之下，文中分別公布了三個和四個重復)
library(Chicdiff)
library(ChicdiffData)
輸入文件
1.Two restriction map information files (“design files”)
- Restriction map file (.rmap)
限制映射文件(.Rmap)-包含限制性片段坐標的底層文件碘菜。默認情況下凹蜈，有4列: chr限寞、 start、 end仰坦、 FragmentID - Bait map file (.baitmap)
誘餌映射文件(履植。Baitmap)-包含帶誘餌的限制性片段的坐標及其相關注釋的底層文件。默認情況下悄晃，有5列: chr玫霎、 start、 end妈橄、 fragmentID庶近、 baitAnnote。文件中指定的區(qū)域(包括它們的片段 ID)必須是眷蚓。Rmap 文件鼻种。BaitAnnote 是一個文本字段，僅用于對輸出和繪圖進行注釋
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 腳本需要兩個位置參數(shù)：
<names-file>：這是一個制表符分隔的文件沙热，包含樣本名稱（第一列）和指向輸入 Rds 文件的完整路徑（第二列）叉钥。這些輸入 Rds 文件應該是 CHiCAGO 分析的輸出結(jié)果。<output-prefix>：這是輸出文件名的前綴篙贸，可以包含文件夾路徑投队。
此外，腳本還接受許多可選參數(shù)爵川，例如 --cutoff（可能用于設置峰值閾值）和 --scorecol（可能用于指定在輸入文件中表示分數(shù)或強度的列）敷鸦。
一個使用這個腳本的命令可能看起來像這樣：
Rscript makePeakMatrix.R --cutoff 5 names_and_paths.tsv output_prefix
在這個命令中，names_and_paths.tsv 是一個包含樣本名稱和對應 Rds 文件路徑的文件寝贡，output_prefix 是輸出文件名的前綴扒披。--cutoff 5 是一個可選參數(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) 文件是在示例級別生成的(與每次重復相反) 洞焙，則不應該命名每個向量中的元素
主要流程
我們在測試數(shù)據(jù)上運行 Chicdiff，首先拯啦，我們使用 setChicdiffExperiment()設置 Chicdiff 實驗設置:
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test")
注意澡匪，除了返回設置列表之外，該函數(shù)還會將其保存為 Rds 文件 * * * _ sets褒链⊙淝椋可重復使用的 Rds。
可選: 您可以直接在命令行中修改設置(參見甫匹？默認設置)甸鸟。例如惦费，我們可以將運行模式切換到并行模式(這會加快運行時速度，但增加內(nèi)存需求) :
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test", settings = list(parallel=TRUE))
或者抢韭，可以在設置文件中提供自定義參數(shù):
settingsFile = file.path(dataPath,"SettingsFile.txt")
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test", settingsFile = settingsFile)
最后薪贫，我們使用chicdiffPipeline()運行流程:
output <- chicdiffPipeline(chicdiff.settings)
流水線將基于 RUtended 選項定義擴展區(qū)域(默認情況下，從軟化相互作用的另一端的每個相互作用峰每個方向5個片段) 刻恭，執(zhí)行標準化瞧省，差異檢驗，基于“測試集”的權(quán)重估計和 p 值加權(quán)和多重測試校正鳍贾。
對于每個bait-region 相互作用鞍匾，輸出數(shù)據(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 中的每個交互保存最終的 count 表以及_countput.Rds. .可以使用 saveAuxFiles = TRUE 設置生成更多的輸出文件-請參見梳码？詳細信息請參閱?chicdiffpipeline。
2.3 Output diagnostic plots
chicdiffPipeline()生成了幾個保存在工作目錄中的圖伍掀。這些為 chr19數(shù)據(jù)集生成的圖是作為 ChicffData 包的一部分便存有的掰茶。請注意，對于全基因組數(shù)據(jù)蜜笤，下面描述的趨勢應該更加明顯濒蒋。
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 值的隱含決策邊界。邊界對協(xié)變量(距離)的依賴性表明把兔，這個協(xié)變量是信息量大的沪伙。低 p 值的趨勢是在低距離富集。
Estimated weights
IHWweightPlot <- png::readPNG(file.path(resultsPath,"test_IHWweightPlot.png"))
grid::grid.raster(IHWweightPlot)
從圖中可以看出县好，權(quán)重與距離(stratum)有關围橡。正如預期的那樣，在數(shù)據(jù)的不同隨機子集(folds)上計算的權(quán)重函數(shù)的行為是相似的缕贡。對于這個 chr19數(shù)據(jù)翁授，跨越較長距離的交互會受到懲罰，賦予以較低的權(quán)重晾咪，而跨越較短距離的交互則被優(yōu)先考慮收擦。
差異的分析的可視化
進行兩種比較條件的可視化，可以使用 plotdiffBaits()函數(shù)繪制相互作用的原始count與它們與相應bait碎片的線性距離的圖像谍倦。對應于重要相互作用的其他端的擴展區(qū)域表示為根據(jù)其調(diào)整后的加權(quán)值進行顏色編碼的區(qū)間塞赂。
Chicdiff 管道自動調(diào)用 plotdefBaits ()為從前100個包含與最低加權(quán) p 值相互作用的誘餌中選擇的4個隨機誘餌生成配置文件。顯示了來自誘餌的1Mb 內(nèi)的相互作用昼蛀，并且根據(jù)分別在每個條件中檢測到的相應相互作用的 Chicago score(score>5: red; 3)對圖上的數(shù)據(jù)點進行顏色編碼宴猾。
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)對單個片段進行優(yōu)先排序
除非 RUexpad 設置設置為零圆存，否則 Chicdiff 工作在誘餌和池區(qū)域之間的交互級別，包含多個限制性片段仇哆。
為了在單個片段獲得差異信號的指示辽剧，Chicdiff 提供了函數(shù) getCandidateInteractions ()。對于落在多個匯集區(qū)域內(nèi)的碎片税产，該函數(shù)通過取最小值(默認值)或更正式地使用調(diào)和平均值方法(由方法參數(shù)確定)來組合它們的微分 p 值怕轿，用于在包調(diào)和平均值中實現(xiàn)的相關檢驗。還可以指定用于篩選輸出的最大 p 值截止值(pvcut)辟拷。
為了優(yōu)先考慮差異相互作用區(qū)域內(nèi)推定的“驅(qū)動器”相互作用撞羽，getCandidateInteractions 允許在條件之間通過最小 | asinh (Chicago 得分) | (minDeltaAsinhScore)過濾它們。Asinh 轉(zhuǎn)換有助于壓縮得分的上限衫冻，在這個范圍內(nèi)诀紊，他們之間的差異比那些低范圍的差異更難解釋。
outCI <- getCandidateInteractions(chicdiff.settings = chicdiff.settings,
output = output, peakFiles = peakFiles,
pvcut = 0.05, minDeltaAsinhScore = 1)
head(outCI)
