基于深度学习的染色质交互作用预测

1、基于深度学习的染色质交互作用预测摘要：染色质交互作用一般发生在生物基因表达过程中，可以对染色质高维空间如染色质环、拓扑结构域等产生影响，通过 特定实验的高通量测序技术如5C、Hi-C等可获得染色质交互作用信息,预测染色质交互作用信息具有重要价值,可 减少测序技术高昂的代价。提出基于组蛋白修饰数据来预测Hi-C数据的深度学习方法，从不同角度对结果进行分 析，表明该方法可以较好预测原始染色质交互作用情况，为组蛋白修饰预测染色质交互作用数据提供一种可行路线。 关键词：组蛋白修饰；Hi-C测序技术；深度学习；染色质交互作用0引言染色质是由DNA、组蛋白、非组蛋白等多种物质组 成的遗传物质，其结构复杂，

2、难以直接观察，但对细胞 遗传过程的基因表达有重要影响。自3C技术问世以 来，众多方法被陆续报道用于捕获染色质构象，其中 Hi-C技术是捕获染色质相互作用频次的最新最常用 方法o Hi-C原始互作数据可以通过交互频次的读取 序列映射到对称矩阵中,并且利用这种矩阵热图可以 表示并构造为染色质的高级结构TAD121、隔间和染色质 环等。染色质的高级结构与其功能密切相关,对基因 表达和生物遗传有重要影响131，如三维结构变化可能诱 导肿瘤发育产生1 41 o目前对染色质结构的研究非常丰富，由于Hi-C数 据测量非常昂贵且耗时，但对与基因表达、转录和疾病 状态相关的染色质结构有重要意义151，所以有很多

3、研究 开始关注减少实验进行预测，对染色质结构研究存在 三维结构构建，二维结构预测以及基因表达和功能性 研究。对染色质三维构建目前存在多种方法,分别使 用了多种距离模型算法来构建三维空间结构问。最新 的三维结构构建方法还可以分析其基因功能，7|o而二 维结构分析同样基于Hi-C数据区域分割域DNA序 列预测，91来构建TAD、染色质环等结构，识别其区域 性。染色质测序技术的发展,还对表观遗传学方面的 实验应用非常重要1 101,表观遗传学包括组蛋白修饰等方 面,对基因表达调控和染色质重塑有重要影响111。分析 组蛋白修饰的功能作用121，研究染色质结构功能和基因 表达131的影响，例如对染色质开

4、放性141和染色质状态的 识别151，具有重要研究价值。利用深度学习预测染色质交互作用能够有效获取 输入数据的前后关联信息和局部特征，目前对人类细 胞类型GM128781 161和果蝇细胞类型171都有相关报道，使 用的数据一般围绕DNA序列和表观遗传学数据。为 了深入研究组蛋白修饰和染色质交互作用两种数据相 关性，本文提出了一种从组蛋白修饰数据中预测人类 基因组中Hi-C数据的方法，基于卷积神经网络，针对 常见人类细胞类型IMR90建立了深度学习模型比较预 测，通过线性相关系数皮尔逊系数以及图相似性系数 等评估，并最终在预测结果与原始结果之间表现出高 相关性。1数据和方法1.1数据预处理Hi

5、-C数据:在GEO公开数据集上，可以通过访问 代码GSE63525下载IMR90细胞类型的数据。我们从 实验原始观测的Hi-C序列交互数据生成原始矩阵数 据,并根据分辨率确定每段基因的长度，计算对应段位 置(例如本文采用10kb分辨率,k为一千数量,b指代 碱基段即一段碱基序列，就会将每个交互位置数据除 以104，获得其bin序号,每条染色体按细胞类型和染 色体不同有上百万或上亿碱基长度)。本文关注研究 染色质内交互作用，即同一序号染色体间的交互作用 (同染色质间交互)，最终获得22条染色体的Hi-C原 始交互作用热图(去除性染色体影响)。组蛋白修饰数据:从Roadmap上可以下载每种细 胞类

6、型的所有表观遗传修饰数据， HYPERLINK https:/egg2.wust https:/egg2.wust!. edu/roadmap/data/byFileType/signal/consolidated/ macs2signal/fO!dChange/。不同的细胞格类型对应有不同的组蛋白修饰。对于 IMR90细胞类型，可以下载以下多种修饰因素：H3K23ac、H3K79Me1、H3K27Ac、H3K79me2、H3K27me2、 H3K79me3、H3K27me3、H3K9acH3K3K36me1、H3K9me1、 HP4、RPD3、H1、H3K36me2、H3K36me3、H3K

7、9me3、H3K- 4me 1、H4、H3、H3、3、3me 1、H4K3、3k4me 1、H4K16ac以上数据可在ENCODE项目中公开访问，首先对 组蛋白修饰的数据进行预处理,保持与Hi-C数据同样 分辨率大小(例如大小为10kb)。1.2实验环境GPU: NVIDIA TU102 GeForce RTX 2080 Ti Rev. A (rev a1)CPU: 48 英特尔至强 CPU E5-2650 v4 2.20GHz 内存:128GBPython版本和依赖包环境：Python 3.6,基于Ten- sorFlow 的 Keras。1.3模型与评估本文通过多种的深度学习方法评估训练结

8、果。训 练神经网络并使其拟合的过程重，使用MSE(Mean Square Error)作为损失函数，使用 MAE (Mean Absolute Error)作为目标函数。线性相关性上使用皮尔逊相关系数(PCC)和斯皮 尔曼相关系数(SPCC)来评估预测结果,结果数值将显 示预测结果与原始数据之间的线性相关性，相关系数 的绝对值越接近1,其相关性就越强。N匕=1皮尔逊相关系数公式: ( YrealredYpred)SPCC是基于PCC的一种相关系数计算方法，给 定变量基础上给出一个排序差异集合d,由两个变量集 合中的每个元素计算，最后使用PCC公式来获取排序 变量的结果。除序列评估外，还可使用计

9、算峰值信噪比(Peak Signal to Noise Ratio)和结构相似度指数(Structural SIMilarity)来进行评估，这两者都常用于图像处理和 去噪。PSNR表示图像信噪比，此值越大代表失真越少， MAX=max Ypred min Ypred 。PSNR = 20*log10( MAX )4MSESSIM值的范围从0到1,衡量两幅图的相似度，判 断预测图像是否接近于原始图像：SSIM -QreidNpred + 1)(realpred + C2)(4 + 说red + 回(以 + 优red + C2)2方法2.1数据预处理对于输入数据和输出数据,本文分别使用不同方 法

10、进行预处理，因为不同数据实验测序方法不同，其表 达值彼此差异较大,难以直接定性分析。Hi-C数据可使用标准化函数将其归一化到0,1 的 区间范围，表示每个位置交互的可能性,其原始序列交 互数据可从Rao等人(GSE63525)1181所做的公开可用的 Hi-C实验中获得。原始数据为每个染色质每段碱基 部分间的交互作用强度，在10k碱基分辨率下，每个作 用强度数据就表示染色质上按顺序排列的两段10k长 度碱基的交互频率。由于每个染色质长度彼此不同， 使用字母标识i和j表示两个段碱基序号，它们之间的 计数nij表示染色质段上相互作用次数,通过碱基位置 对应关系可构成对称矩阵，矩阵大小为N表示染色质

11、 长度L/分辨率R(如图1所示)。对于组蛋白数据，首先通过bwtools和指定bin长度（分辨率）生成组蛋白修饰序列信号数据矩阵H,此 分辨率与Hi-C数据分辨率一致,且进行截取使实验数 据长度相同，矩阵H中每列为组蛋白修饰类型，共M 列，使用最大最小值归一化使数据在0,1区间，每列数 据表示组蛋白修饰与染色质产生作用的可能性。图1 Hi-C交互作用热图图1 Hi-C交互作用热图神经网络介绍：卷积神经网络（CNN）一般用于图 像处理或自然语言处理等高维特征自动提取,可学习 到平移不变性等特征，在生物信息研究上也广泛应 用。利用CNN可以快速提取相应染色质交互的相关 因子特征序列，构建交互作用概

12、率分布输入。在全连 接层部分,使用循环神经网络考虑染色质上下游数据 带来的影响，并使所有神经元之间的参数共享，这些参 数针对染色质基相互作用进行了优化,可生成用于不 同组蛋白修饰的最佳滤波器，均方误差（MSE）作为损 失函数，平均绝对误差（MAE）作为目标函数。对输入 输出数据采用取对数值获取-0,1区间值域范围,并作 为神经网络模型的标签和样本，可理解为对每段染色 质交互作用的可能性预测。本文提出的模型基于每对序列对应的方式,对矩 阵中数据一一预测，如图3所示，使用一个w大小的窗 口来获取每个交互基因距离下的所有数据，由于Hi-C 矩阵中的对称特性,交互数据可以按列或行获取。因 为染色质的交

13、互和高维结构受到碱基段的上下游影 响，有明显的区域性，对于输入的组蛋白修饰序列数 据,每个交互位点的上下游各一个碱基段作为输入，共 三个碱基段长度。因此，每个样本关注w大小的染色 质交互作用，使用x-1到x+1段（x为交互作用发生的 位置）的3个长度的组蛋白修饰作为输入来预测相应 的Hi-C交互作用情况。本文对模型输入部分进行划 分,构建一个输入模块获取两对不同位置的输入数据, 此模型将在Hi-C矩阵中的对角线区域附近生成每个 bin的交互作用数据。最后结果用热图重建方法来还 原预测矩阵。窗口大小w的选取,可使用的Hi-C矩阵 为10kb（ 104 ）分辨率下的实验数据，设定w=50,因此 基

14、因组距离为500kb,即每个碱基段包含500k （500X 103 ）的碱基。这样可以观察交互作用密集区 域,排除交互作用发生不明显的区域和较远距离的稀 疏数据,使得预测结果更有价值。Bini Input | Bin2 Input| ConvlD(k：rnel=32) Bini Input | Bin2 InputIBatchNormalizationI |BatchNormalization| : :| ConvlD (k：rne1=64厂| | ConvlD(kernel=64门 : : ； |ConvlD(kerne1=128) | |ConvlD(kernel=128)|ConvlD(

15、kernel=l)ConvlD(kernel=32)ConvlD(kernel=l)BatchNormalizationConvlD(kernel=l)ConvlD(kernel=32)ConvlD(kernel=l)ConvlD(kernel=64) BatchNormalization| ConvlD (kernel=128)|BatchNormalization| ConvlD (kernel=128)ConvlD (kernel=l)Dense(w*8)Dense(w*4)Dense(w*2)RNN(w*2)Dense(=w)Dense(=1)图2染色质交互作用预测神经网络模型结构图图

16、3两段bins输入模型预测过程包括三个阶段，第一阶段为输入数 据的卷积和特征聚合阶段，对每两对输入的组蛋白修 饰信号矩阵，通过多层一维卷积神经网络获取其多种 类型修饰的一维聚合特征序列，代表不同修饰类型共 同作用的综合权重分布。第二阶段为拼接层,将两对 输入构成二维矩阵，再次利用卷积神经网络进行不同 位置间的交互作用影响特征提取。最后一层为全连接 层和循环神经网络层，利用碱基上下游序列的影响，转 化为时序序列的前后特征，来预测相应染色质碱基段 的交互作用结果。由于神经网络预测过程中是针对各 个位点进行的，需要根据位置和对称性重建矩阵结果 进行对比评估，预测结果中可得到多段w大小的预测 数值排列

17、，根据提取过程的顺序可以依次对应到交互 作用发生的位置上，并使用重构算法将其恢复成数值 矩阵。因此，最终结果也是对角线区域在指示基因组 距离上的完整矩阵，并存在一定的拓扑结构区域。重构矩阵热图伪代码：M为最后结果矩阵For i in N:if i N - w + 1: Mi, i: i+w = Input1:w; Mi: i+w, i=M i, i: i+welse: x = N - I; Mi, i: i+w = Input1:x; Mi: i+w, i=Mi, i: i+w3结果我们对模型进行了多次训练，设定超参数为训练 轮次设定为30轮,批次大小每次100个样本，优化器 是RMSProp

18、。训练完毕后可以获得一维滤波器的各项 权重，表示每对固定组蛋白修饰序列的综合作用特 征。最终结果分布表示序列相关性和矩阵相似性。最终结果分别使用线性相关和图像相似性进行分 析。数据集使用GEO数据库中的GSE63525访问代码 获取的10kb分辨率的IMR90原始交互作用数据，重构 成交互作用矩阵。其中训练集使用1-17号染色体的 Hi-C样本，测试集使用18-22号染色体的Hi-C样 本。最后结果显示，在500kb的基因组距离下，对Hi- C样本进行预测结果分析。线性相关性分析中测试集 PCC最好达到0.85,SPCC最好达到0.8，SSIM的测试 样本中最好可达到0.98。表1细胞类型IMR90染色质编号的测试数据集为染色体18号-22号/基因组距离（500kb）染色体序号整体PCC整体SPCC平均SS