Capture Hi-C 是一种检测至少在一端涉及感兴趣的 DNA 区域，如基因启动子染色体相互作用的有效方法。 Chicdiff，是一个用于捕获 Hi-C 数据中差分相互作用的鲁棒检测的 R 软件包。Chicdiff 增强了计数数据的最先进的差异检验方法，具有定制的标准化和多种检验程序，可以解释 Capture Hi-C 的特定统计特性。作者用Chicdiff 在人单核细胞和 CD4 + T 细胞中发表的promoter-Capture Hi-C 数据，鉴定了多种细胞类型特异性相互作用，并证实了启动子相互作用和基因表达之间的总体正相关性。
install.packages("devtools")
library(devtools)
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="Chicdiff")
install_github("RegulatoryGenomicsGroup/chicdiff", subdir="ChicdiffData", force=T)
Chicdiff 是一种在捕获 Hi-C 数据中检测统计显著差异相互作用的方法。这个文档将引导你通过一个典型的Chicdiff analysis.。
简而言之，Chicdiff 将 DESeq2中实现的计数数据的适度差异检验与用于信号归一化和 p 值加权的 CHi-C 特定程序相结合。为了提高威力，Chicdiff 还将每个诱饵的每个显著相互作用区域周围的几个片段的读数结合起来。更具体地说，Chicdiff 使用来自Chicago的背景估计来构建用于 DESeq2差异检验的定制缩放矩阵，来自 DESeq2的 Wald test检验产生 p 值提交给加权过程，如下所述。
由于 CHi-C 数据中的信噪比和效应大小与距离有很强的相关性，Chicdiff 进行非均匀的多重检验校正，因此较长距离的较弱信号得到更有力的校正。它使用 IHW 软件包以交互距离作为协变量来学习 p 值上的权重，使得被拒绝的零假设的总分数最大化。由于提交用于 Chicdiff 测试的一组相互作用通常由于选择偏差而具有不均匀的 Chicdiff p 值分布，因此从整个 Capture Hi-C 数据集采样的“训练集”的相互作用被用于权重训练。通过这种方式学到的权重然后被应用到检验集中的 p 值，并且得到的加权 p 值被调整用于多个检验并报告给用户。
来自 ChicaffData 的示例数据集使用约0.6 Gb RAM，在标准笔记本电脑上处理需要几分钟。在每种条件下，两个全基因组 CHi-C 数据的生物学重复的典型工作需要约30-60分钟，并使用约10-15 Gb RAM。具有更多重复的高覆盖率数据集将需要更长的时间来处理和使用更多 RAM。
该软件包只包含来自两个重复的单核细胞和初始 CD4 + T 细胞的19号染色体数据(相比之下，文中分别公布了三个和四个重复)
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。文件中指定的区域(包括它们的片段 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")
- Count data files in Chicago input format
.chinput 文件，chicagoTools/bam2chicago.sh转换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) 文件是在示例级别生成的(与每次重复相反) ，则不应该命名每个向量中的元素
主要流程
我们在测试数据上运行 Chicdiff，首先，我们使用 setChicdiffExperiment()设置 Chicdiff 实验设置:
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test")
注意，除了返回设置列表之外，该函数还会将其保存为 Rds 文件 * * * _ sets。可重复使用的 Rds。
可选: 您可以直接在命令行中修改设置(参见？默认设置)。例如，我们可以将运行模式切换到并行模式(这会加快运行时速度，但增加内存需求) :
chicdiff.settings <- setChicdiffExperiment(designDir = testDesignDir, chicagoData = chicagoData, countData = countData, peakfiles = peakFiles, outprefix="test", settings = list(parallel=TRUE))
或者，可以在设置文件中提供自定义参数:
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 选项定义扩展区域(默认情况下，从软化相互作用的另一端的每个相互作用峰每个方向5个片段) ，执行标准化，差异检验，基于“测试集”的权重估计和 p 值加权和多重测试校正。
对于每个bait-region 相互作用，输出数据表列出标准化相互作用读取计数的 log2FoldChange，以及原始的，加权的和加权调整的 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
除了返回输出数据表之外，chicdiffPipeline还将其保存为 Rds文件如 _ result.Rds.它还将为 Rds 文件 * _ count put 中的每个交互保存最终的 count 表以及_countput.Rds. .可以使用 saveAuxFiles = TRUE 设置生成更多的输出文件-请参见？详细信息请参阅?chicdiffpipeline。
2.3 Output diagnostic plots
chicdiffPipeline()生成了几个保存在工作目录中的图。这些为 chr19数据集生成的图是作为 ChicffData 包的一部分便存有的。请注意，对于全基因组数据，下面描述的趋势应该更加明显。
resultsPath <- file.path(dataPath, "CD4_Mono_results")
IHWdecisionBoundaryPlot <- png::readPNG(file.path(resultsPath,"test_IHWdecisionBoundaryPlot.png"))
grid::grid.raster(IHWdecisionBoundaryPlot)
image.png
这张图显示了作为协变量函数的未加权 p 值的隐含决策边界。边界对协变量(距离)的依赖性表明，这个协变量是信息量大的。低 p 值的趋势是在低距离富集。
Estimated weights
IHWweightPlot <- png::readPNG(file.path(resultsPath,"test_IHWweightPlot.png"))
grid::grid.raster(IHWweightPlot)
image.png
从图中可以看出，权重与距离(stratum)有关。正如预期的那样，在数据的不同随机子集(folds)上计算的权重函数的行为是相似的。对于这个 chr19数据，跨越较长距离的交互会受到惩罚，赋予以较低的权重，而跨越较短距离的交互则被优先考虑。
差异的分析的可视化
进行两种比较条件的可视化，可以使用 plotdiffBaits()函数绘制相互作用的原始count与它们与相应bait碎片的线性距离的图像。对应于重要相互作用的其他端的扩展区域表示为根据其调整后的加权值进行颜色编码的区间。
Chicdiff 管道自动调用 plotdefBaits ()为从前100个包含与最低加权 p 值相互作用的诱饵中选择的4个随机诱饵生成配置文件。显示了来自诱饵的1Mb 内的相互作用，并且根据分别在每个条件中检测到的相应相互作用的 Chicago score(score>5: red; 3)对图上的数据点进行颜色编码。
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)
image.png
在差异相互作用区域内对单个片段进行优先排序
除非 RUexpad 设置设置为零，否则 Chicdiff 工作在诱饵和池区域之间的交互级别，包含多个限制性片段。
为了在单个片段获得差异信号的指示，Chicdiff 提供了函数 getCandidateInteractions ()。对于落在多个汇集区域内的碎片，该函数通过取最小值(默认值)或更正式地使用调和平均值方法(由方法参数确定)来组合它们的微分 p 值，用于在包调和平均值中实现的相关检验。还可以指定用于筛选输出的最大 p 值截止值(pvcut)。
为了优先考虑差异相互作用区域内推定的“驱动器”相互作用，getCandidateInteractions 允许在条件之间通过最小 | asinh (Chicago 得分) | (minDeltaAsinhScore)过滤它们。Asinh 转换有助于压缩得分的上限，在这个范围内，他们之间的差异比那些低范围的差异更难解释。
outCI <- getCandidateInteractions(chicdiff.settings = chicdiff.settings,
output = output, peakFiles = peakFiles,
pvcut = 0.05, minDeltaAsinhScore = 1)
head(outCI)
网友评论