3.1.1数据矩阵的构造

preciseTAD要求用户提供“ground-truth”域边界的基因组坐标，以建立响应向量(Y)．作为一个例子，我们考虑从流行派生的TAD边界箭头TAD调用器，是Aiden实验室开发的Juicebox工具套件的一部分(Durand et al. 2016)．边界,箭头在常染色体上应用5 kb GM12878 Hi-C数据(Rao et al. 2014)（GSE63525)．用于调用TADs的命令行脚本示例箭头，以及结果TAD坐标的前三列。关于实现的更详细的教程箭头可以找到在这里．

箭头-c 1 \ #染色体调用TADs在-r 5000 \ #HiC数据分辨率~/GSE63525_GM12878_insitu_primary。hic文件的位置~/arrowhead_output #存储输出的位置

数据(“arrowhead_gm12878_5kb”)头(arrowhead_gm12878_5kb)

## v1 v2 v3 ## 1 1 49375000 50805000 ## 2 1 16830000 17230000 ## 31 163355000 164860000 ## 4 1 231935000 233400000 ## 5 1 149035000 149430000 ## 6 1 3995000 5505000

方法将此TAD坐标数据转换为惟一边界坐标的grange对象preciseTAD: extractBoundaries ()函数。注意，此函数的输入是一个类似bed的数据帧，其中第二列和第三列是域锚点中点的基因组坐标。这里我们只提取CHR1和CHR22的边界(CHR = c("CHR1"， "CHR22")参数)。我们指定预处理= FALSE因为我们只对所有边界感兴趣，而不是根据长度过滤TADs。最后,我们指定分辨率= 5000来匹配使用的分辨率箭头(尽管这个参数会被忽略预处理= FALSE)．如下图所示，共报告了1901个独特的TAD边界箭头1号染色体和22号染色体。

范围< extractBoundaries(域。mat = arrowhead_gm12878_5kb, filter = FALSE, CHR = c("CHR1"， "CHR22")， resolution = 5000) #查看唯一边界边界

## seqnames ranges string# #    ## 787 chr1 815000-815001 * ## 9196 chr1 890000-890001 * ## 37 chr1 915000-915001 * ## 8923 chr1 1005000-1005001 * ## 275 chr1 1015000-1015001 * ## ... ... ... ...## 8379 chr22 50815000-50815001 * ## 16788 chr22 50935000-50935001 * ## 8382 chr22 50945000-50945001 * ## 16791 chr22 51060000-51060001 * ## 16670 chr22 51235000-51235001 * ## ------- ## seqinfo: 22个序列来自一个未指定的基因组;没有seqlengths

为了确定最能预测TAD边界的基因组特征，preciseTAD需要功能基因组注释数据。它们被用来建立特征空间(\(\textbf{X}=\{X_{1}， X_{2}， \cdots, X_{p} \}\))．细胞类型特异性基因组注释数据可从编码数据门户作为BED文件。一旦你下载了你喜欢的功能基因组注释列表，将它们存储在一个特定的文件位置(即，“。/path/to/BEDfiles”)。然后可以将这些文件转换为GRangesList对象，并用于使用preciseTAD: bedToGRangesList ()函数。的信号参数指BED文件中包含峰值信号强度值的列，用于将元数据分配给相应的grange(仅用于下游绘图)。我们已经提供了一个GRangesList对象，其中包含26个特定于GM12878细胞类型的转录因子结合位点(TFBS)。加载它之后，可以使用以下命令查看转录因子列表。

# path <- "。# tfbsList <- bedToGRangesList(filepath=path, bedList=NULL, bedNames=NULL, pattern = "*. "bed"， signal=4) data("tfbsList") names(tfbsList)

# #[1]“Gm12878-Atf3-Haib”“Gm12878-Cfos-Sydh”“Gm12878-Cmyc-Uta”“Gm12878-Ctcf-Broad”“Gm12878-Egr1-Haib”“Gm12878-Ets1-Haib”“Gm12878-Gabp-Haib”“Gm12878-Jund-Sydh”“Gm12878-Max-Sydh”“Gm12878-Mazab85725-Sydh”“Gm12878-Mef2a-Haib”“Gm12878-Mxi1-Sydh”“Gm12878-P300-Sydh”“Gm12878-Pol2-Haib”“Gm12878-Pu1-Haib”# #[16]“Gm12878-Rad21-Haib”“Gm12878-Rfx5-Sydh”“Gm12878-Sin3a-Sydh”“Gm12878-Six5-Haib”“Gm12878-Smc3-Sydh”“Gm12878-Sp1-Haib”“Gm12878-Srf-Haib”“Gm12878-Taf1-Haib”“Gm12878-Tr4-Sydh”“Gm12878-Yy1-Sydh”“Gm12878-Znf143-Sydh”

tfbsList

# # GRangesList对象长度26:# # $“Gm12878-Atf3-Haib”与1677年# #农庄组织对象范围和1元数据列:# # seqnames链|覆盖范围# # < Rle > < IRanges > < Rle > | <数字> # # [1]chr1 1167429 - 1167429 * | 28.5019 # # [2] chr1 1510215 - 1510215 * | 64.7440 # # [3] chr1 1655730 - 1655730 * | 25.1686 # # [4] chr1 1709843 - 1709843 * | 15.5525 # # [5] chr1 2343938 - 2343938 * 19.6961 |  ## ... ... ... ... . ...## [1673] chrX 118699292-118699462 * | 17.6970 ## [1674] chrX 119947928-119948098 * | 13.8884 ## [1675] chrX 128481457-128481627 * | 26.1075 ## [1676] chrX 132399481-132399651 * | 20.7070 ## [1677] chrX 136520393-136520563 * | 24.4197 ## ------- ## seqinfo:来自未指定基因组的24个序列;没有seqlength ## ##…## <25个元素>

使用“ground-truth”边界和下面的TFBS，我们可以构建用于预测建模的数据矩阵。的preciseTAD: createTADdata ()函数可以创建训练和测试数据。这里，我们指定在染色体1上进行训练，在染色体22上进行测试。此外,我们指定分辨率= 5000构建5kb移位基因组箱(以匹配Hi-C数据分辨率)，featureType = "距离"对于一个\ (log_2(距离+ 1)\)-type特征空间和重采样=“罗斯”在训练数据上应用随机欠采样(RUS)来平衡边界和非边界区域的分类。我们还指定一个种子，以确保执行重采样时的可重复性。结果是一个包含两个数据帧的列表:(1)重采样(如果指定)训练数据，(2)测试数据。

set.seed(123) tadData <- createTADdata(bounds)。GR = bounds，分辨率= 5000，基因组元素。GR = tfbsList, featureType = "distance"， resampling = "rus"， trainCHR = "CHR1"， predictCHR = "CHR22") #训练数据的视图子集tadData[[1]][1:5,1:4]

## y Gm12878-Atf3-Haib Gm12878-Cfos-Sydh Gm12878-Cmyc-Uta ## 5799 No 18.73205 17.73674 14.60901 ## 32930 No 21.23327 15.32696 15.03635 ## 5767 Yes 18.07127 17.99360 13.22415 ## 8504 Yes 19.24515 19.26026 18.31016 ## 39486 No 21.23396 17.28550 19.09578

查看tadData[[1]]$y是否平衡

## ##否是## 1595 1595

查看测试数据tadData[[2]][1:5,1:4]

## y Gm12878-Atf3-Haib Gm12878-Cfos-Sydh Gm12878-Cmyc-Uta ## 1 No 24.08978 24.1100 24.12327 ## 2 No 24.08937 24.1096 24.12287 ## 3 No 24.08897 24.1092 24.12248 ## 4 No 24.08857 24.1088 24.12208 ## 5 No 24.08816 24.1084 24.12169

3.1.2使用递归特征消除进行特征选择

现在我们可以实现机器学习算法来预测边界区域。在这里，我们选择随机森林算法进行二进制分类。preciseTAD提供了执行递归特征消除(RFE)的功能，作为特征减少的一种形式，通过使用preciseTAD: TADrfe ()函数的包装器rfe功能脱字符号包库恩(2012)．preciseTAD: TADrfe ()使用5倍交叉验证，在数据中从2到最大数量的特征的最佳子集上实现随机森林模型。我们指定精度作为性能度量。下面是一个示例。

set.seed(123) rfe_res <- TADrfe(trainData = tadData[[1]]， tuneParams = list(ntree = 500, nodesize = 1)， cvFolds = 5, cvMetric = "Accuracy"， verbose = TRUE)

查看RFE性能rfe_res[[1]]

##变量MCC ROC Sens Spec Pos Pred Value Neg Pred Value Accuracy Kappa MCCSD ROCSD SensSD SpecSD Pos Pred Value esd Neg Pred Value Neg Pred Value Accuracy KappaSD ## 12 0.3491082 0.7420849 0.6526646 0.6952978 0.6822854 0.6667258 0.6739812 0.3479624 0.02627631 0.010098192 0.009508308 0.02993689 0.019873672 0.013160545 0.02632109 ## 24 0.4613042 0.8026680 0.6984326 0.7611285 0.7455349 0.7163087 0.7297806 0.4595611 0.009974699 0.015099178 0.02275698 0.015584018 0.0079491990.008769029 0.01753806 ## 38 0.5076806 0.8229440 0.7122257 0.7931034 0.7749194 0.7340188 0.7526646 0.5053292 0.01922275 0.010486133 0.022648767 0.01063061 0.008681254 0.014664805 0.010711191 0.02142238 ## 416 0.5150537 0.8263765 0.7197492 0.7943574 0.7779659 0.7394767 0.7570533 0.5141066 0.01881406 0.014157078 0.022886160 0.01696849 0.012896280 0.014149142 0.010802546 0.02160509 ## 526 0.5307963 0.8294258 0.7894162 0.74392633 0.5285266 0.02605025 0.0152340240.036287890 0.01480341 0.009062094 0.022029556 0.013957766 0.02791553

#查看每个CV折叠头前n个特征的变量重要性(rfe_res[[2]])

##全局变量Resample ## 1 26.76923 Gm12878-Smc3-Sydh 26 Fold1 ## 2 23.70939 Gm12878-Rad21-Haib 26 Fold1 ## 3 20.08605 Gm12878-Ctcf-Broad 26 Fold1 ## 4 17.65042 Gm12878-Znf143-Sydh 26 Fold1 ## 5 14.17874 Gm12878-Pol2-Haib 26 Fold1 ## 6 13.25481 Gm12878-Gabp-Haib 26 Fold1 ## 3

递归特征剔除结果表明，当仅考虑构建随机森林预测模型时的前四个转录因子时，模型精度开始趋于稳定(图1A)。将每个交叉褶皱中前四个TFBS的可变重要度值进行聚合后，我们可以看到，五个褶皱中最重要的前四个TFBS分别是SMC3、RAD21、CTCF和ZNF143(图1B)。这些都是已知的环挤压模型的组成部分，该模型被提出作为人类基因组三维结构的一种机制(Sanborn et al. 2015;Fudenberg et al. 2016;汉森等人。2018年)．

图1．(A)递归特征消除(RFE)表明，当只考虑前4个最具预测性的转录因子结合位点(TFBS)时，性能稳定。(B) 5个交叉折叠中前4个TFBS的总平均变量重要性(使用平均精度递减)。

3.1.3实现一个用于边界预测的随机森林

既然我们已经适当地减少了特征空间，我们就可以在上面提到的TFBS (SMC3、RAD21、CTCF和ZNF143)上简单地实现一个随机森林算法。我们可以利用preciseTAD: TADrandomForest ()函数的包装器randomForest包(Breiman 2001;利奥，维纳和其他人2002)．我们指定了训练和测试数据、超参数值、交叉验证折叠数、要考虑的性能度量(这里是准确性)、要初始化的可再现性种子、保留模型对象的指标、保留变量重要性的指标、要考虑的变量重要性度量(这里是精度的平均下降(MDA))，以及基于测试数据保留模型性能的指标。函数返回一个列表，其中包含:tadModel [[1]]), 2)一个data.frame模型中包含的每个特征的重要性是可变的(tadModel [[2]])和3)adata.frame各种模型的性能指标。

#限制数据矩阵只包括SMC3, RAD21, CTCF和ZNF143 tfbsList_filt <- tfbsList[names(tfbsList) %in% c("Gm12878-Ctcf-Broad"， "Gm12878-Rad21-Haib"， "Gm12878-Smc3-Sydh"， "Gm12878-Znf143-Sydh")] set.seed(123) tadData <- createTADdata(bounds. txt)GR = bounds，分辨率= 5000，基因组元素。GR = tfbsList_filt, featureType = "distance"， resampling = "rus"， trainCHR = "CHR1"， predictCHR = "CHR22") #运行RF set.seed(123) tadModel <- TADrandomForest(trainData = tadData[[1]]， testData = tadData[[2]]， tuneParams = list(mtry = 2, ntree = 500, nodesize = 1)， cvFolds = 3, cvMetric = "Accuracy"， verbose = FALSE, model = TRUE, importances = TRUE, impMeasure = "MDA"， performance = TRUE)

#查看模型性能性能<- tadModel[[3]]性能$Performance <- round(性能$Performance, digits = 2) rownames(性能)<- Performance $Metric kable(t(性能)，caption = "在CHR22测试数据上验证基于CHR1构建的射频时的模型性能列表。")

您可能知道，还有其他机器学习二进制分类器可以用于这种设置。例如，假设我们选择实现支持向量机(SVM)。鉴于此，这很容易实现preciseTAD: createTADdata ()方便地设置训练和测试数据集。我们使用e1071包来使用下面的示例命令运行具有径向内核、cost = 1和gamma = 0.5的SVM。我们可以看到，支持向量机模型的准确率只有0.67，而我们的随机森林的准确率是0.69(表1)。

svmModel <- svm(y ~ .， data = tadData[[1]]， kernel = "radial"， cost = 1, gamma = 0.5) svmPreds <- predict(svmModel, tadData[[2]][， -1]， positive = "Yes") #查看混淆矩阵表(svmPreds, tadData[[2]][， 1])

## ## svmPreds否是##否6278 49 ##是3076 257

3.2.1之上运行preciseTAD

假设我们想精确预测CHR22:25,500,000-27,500,000的2Mb部分的域边界坐标。为此，我们指定chromCoords = list(25500000, 27500000)．此外，我们设置了一个用于构建ptbr的概率阈值范围，包括0.975、0.99和1.0。对于DBSCAN，我们分配的范围为\ε(\ \)-邻域值包括5000、10000、15000、20000和30000，并保持3为MinPts价值。我们指定verbose = TRUE这样函数就会输出一些小的结果\ ((t) \ε)\)以及整体最佳组合

#运行preciseTAD set.seed(123) pt <- preciseTAD(genome elements . txt)GR = tfbsList_filt, featureType = "distance"， CHR = "CHR22"， chromCoords = list(17000000, 18000000)， tadModel = tadModel[[1]]， threshold = 1, verbose = TRUE, parallel = NULL, DBSCAN_params = list(eps = c(30000)， # c(10000, 30000, 50000) MinPts = c(100))， # c(3,10,100) slope = 5000, genome = "hg19"， BaseProbs = TRUE)

##[1]“使用距离类型特征空间建立bp分辨率测试数据”##[1]“建立概率向量”##[1]“确定epsilon邻域(eps)和最小点数(MinPts)的最佳组合”##[1] "初始化DBSCAN for MinPts=100和eps=30000" ## [1] " preciseTAD识别2个ptbr " ##[1] "建立ptbr " ## [1] " Cluster 1 of 2" ## [1] " Cluster 2 of 2" ## [1] " MinPts和eps的最佳组合是= (100,30000)"## [1] "preciseTAD识别出总共1678个碱基对，其预测概率等于或超过1的阈值" ##[1]"初始化DBSCAN为MinPts = 100和eps = 30000" ## [1] "preciseTAD识别出2个ptbr " ##[1] "建立ptbp " ##[1] " 2个集群中1个" ##[1]" 2个集群中2个"

#返回的是什么?ptbr .精确TAD边界区域和ptbp .边界点名称(pt)

## [1] "preciseTADparams" "PTBR" "PTBP" "摘要"

#查看结果# pt

当指定verbose = TRUE该函数为所有人提供简短的输出\ε\ \ (t \倍)组合包括每个组合产生多少ptbr / ptbp。注:建议提供verbose = FALSE在实现preciseTAD使用多个\ \ (t)而且\ε(\ \)在整个染色体上。此外，我们从详细的表述中看到\((t， \epsilon) = (1,30000)\)．

的全部输出preciseTAD: preciseTAD ()函数是一个包含4个元素的列表:1)t和eps的每个组合的平均(和标准差)归一化充实(NE)值的数据帧(仅当至少为参数提供多个值时;所有后续总结应用于(t, eps))、2)跨越每个preciseTAD预测区域(PTBR)的基因组坐标、3)preciseTAD预测边界点(PTBP)的基因组坐标的最优组合。4)包括下列统计数字的名单:PTBRWidth - PTBR宽度，PTBRCoverage -概率超过阈值的基础水平坐标与PTBRWidth的比值，DistanceBetweenPTBR -前一个PTBR结束与后续PTBR开始之间的基因组距离，NumSubRegions -每个PTBR集群中的元素数量，SubRegionWidth -跨越与每个PTBR相关的子区域的基因组坐标，DistBetweenSubRegions -前一个ptbr特异性区域结束和后续ptbr特异性区域开始之间的基因组距离，normmilizedenrichment -在模型中使用的围绕侧面ptbp的基因组注释的规范化富集，以及BaseProbs -每个对应的基础坐标的概率的数字向量。归一化富集计算方法为:与两侧预测边界重叠的峰值区域总数除以预测边界数量。图2给出了如何计算每个摘要的示意图。

图2．图解如何计算每一个诊断摘要preciseTAD ()输出。

作为边界的基本特定概率向量(BaseProbs)可能很大，取决于chromCoords设置。在我们的示例中，它有200万个条目。设置时要小心BaseProbs = TRUE以避免内存问题。为了了解前500个碱基的概率是怎样的，我们这样做:

情节(pt总结BaseProbs美元[1∶])

我们可以看到，成为边界的概率似乎在200和400个碱基附近达到峰值。它远低于1，但可能值得进一步研究。

3.2.2使用preciseTAD与Juicebox

Juicebox是艾登实验室提供的一个接口，允许在Hi-C接触图上叠加边界坐标。要可视化被预测边界包围的域，必须首先选择Hi-C映射。例如，您可以通过选择导入Rao等人2014推导的GM12878的联系矩阵文件- >打开……-> Rao和Huntley等人。-> GM12878 -> in situ Mbol -> primary．格式preciseTAD结果在Juicebox中使用，用户可以利用preciseTAD: juicer_func ()中可用的函数preciseTADR包，它将GRanges对象转换为如下所示的数据帧。

#转换pt_juice <- juicer_func(pt$PTBP)

然后，您需要将ptbp保存为BED文件write.table如下所示。保存后，将它们导入到Juicebox使用显示注释面板-> 2D注释->添加Local -> pt_juice.bed．

Filepath = "~/path/to/store/ptbps"写入。表(pt_juice、文件。path(filepath， "pt_juice.bed")， quote = FALSE, col.names = FALSE, row.names = FALSE, sep = "\t")

3.2.3Cross-cell-type预测

preciseTAD允许我们使用预先训练的模型来预测细胞类型中区域边界的精确位置，没有Hi-C数据，但有基因组注释数据。具体来说，只需要CTCF、RAD21、SMC3和ZNF143转录因子结合位点的细胞类型特异性ChIP-seq数据(BED格式)。这些数据需要使用前面提到的方法组装到grange对象中bedToGRangesList ()函数(例如，创建另一个特定于单元格类型的tfbsList_filt对象)。

pre-trained模型tadModel可以通过使用bedToGRangesList ()，extractBoundaries ()，createTADdata (),TADrandomForest ()功能描述。或者，染色体特定的预先训练模型可以从OneDrive存储 ．我们推荐使用GM12878基因组注释数据和Peakachu峰预训练的模型。注意，每个模型都注释了“CHRi_holdout”，表示在训练中保留的染色体，因此可以用来预测边界。需要注意的是，训练有素的preciseTAD以肿瘤细胞系GM12878和K562训练模型。虽然我们预计区域边界形成的基本规则(CTCF等结合)将保留在癌细胞中，但在应用预训练的模型预测正常细胞系的边界时要谨慎。

我们以Hela-S3细胞系为例。我们将下载并使用CHR22_GM12878_5kb_Peakachu模型预测Hela-S3细胞系中chr22的边界。

#读取基于chr1-chr21 tadModel构建的RF模型<- readRDS("~/Downloads/ chr22_gm12878_5kb_peakchu .rds")

hela - s3特定的bed格式TFBS数据可以在以下网站下载编码．所选的TFBS数据可从辅助存储库．它们可以转换成GRangesList。

#从存储库中读取CTCF、RAD21、SMC3和ZNF143中的hela基因组注释，并存储为列表CTCF <- as.data.frame(read.delim("http://raw.githubusercontent.com/stilianoudakis/preciseTAD_supplementary/master/data/bed/hela/ctcf.bed"， sep = "\t"，报头= FALSE)) RAD21 <- as.data.frame(read.delim("http://raw.githubusercontent.com/stilianoudakis/preciseTAD_supplementary/master/data/bed/hela/rad21.bed"， sep = "\t"，smc3 <- as.data.frame(read.delim("http://raw.githubusercontent.com/stilianoudakis/preciseTAD_supplementary/master/data/bed/hela/smc3.bed"， sep = "\t"，报头= FALSE)) znf143 <- as.data.frame(read.delim("http://raw.githubusercontent.com/stilianoudakis/preciseTAD_supplementary/master/data/bed/hela/znf143.bed"， sep = "\t"，报头= FALSE)) helaList <- list(ctcf, rad21, smc3, znf143)GR <- bedToGRangesList(bedList=helaList, bedNames=c("Ctcf"， "Rad21"， "Smc3"， "Znf143")， signal=5)

现在我们准备预测Hela-S3在chr22上的边界。

#运行preciseTAD set.seed(123) pt <- preciseTAD(genome elements . txt)GR =海拉。GR, featureType = "distance"， CHR = "CHR22"， chromCoords = list(25500000, 27500000)， tadModel = tadModel[[1]]， threshold = 1.0, verbose = FALSE, parallel = NULL, DBSCAN_params = list(30000, 100)， slope = 5000, BaseProbs = FALSE, savetobed = FALSE)