Sequential regulatory activity prediction across chromosomes with convolutional neural networks
Sequential regulatory activity prediction across chromosomes with convolutional neural networks
Sequential regulatory activity prediction across chromosomes with convolutional neural networks 基于卷积神经网络的染色体序列调控活动预测
摘要
基因预测表型
从DNA序列预测大型哺乳动物基因组中细胞类型特异性（见CSDN）的表观遗传和转录谱。
利用CNN对启动子和远端调控元件进行识别，综合其内容，进行基因表达的预测
我们发现，对基因组变异对基因表达影响的模型预测与人类群体中eQTLs基础的因果变异很好地吻合，并可用于生成机制假设，以实现疾病位点的精确定位。
背景
尽管许多研究表明，在一系列人类疾病和特征中，基因型和表现型的变异之间存在很强的关系，但这种关系运作的机制仍不完全清楚(Boyle等，2017)。非编码变异尤其抑制了进展;通过全基因组关联研究(GWAS)，大多数与表型统计相关的基因组位点不会改变编码序列，但只有极少数的机制已被彻底描述。 大量证据表明，许多非编码变异通过改变基因表达影响性状。反过来，基因表达决定了多细胞生物中细胞类型和状态的多样性(表观基因组学研究路线图等，2015)。因此，基因表达提供了一个易于处理的中间表现型，而改进的模型将有很大的价值。大型财团和许多个体实验室绘制了多种细胞的表观遗传和表达谱。 此外，最近人们认识到，许多数据可以使用机器学习作为潜在DNA序列的功能精确建模。 成功的转录因子(TF)结合预测模型、易访问的染色质和组蛋白修饰为基因组变异提供了机制洞察力和有用的解释。
尽管取得了这一进展，但在复杂的生物体中，从DNA序列预测细胞类型特异性基因表达的模型仍然难以找到。 现有的模型使用实验性的注释作为输入(例如，各种已知的监管属性的峰值调用)，允许它们阐明这些注释之间的关系(Cheng et al. 2012;González等人。2015)，但没有分析潜在序列在确定这些注释时的因果作用。 即使有实验内的训练数据来推断相关的序列基序，远端调控的复杂性(增强子可以与启动子在几十万个核苷酸上相互作用)挑战了当前的方法(Levine 2010;Long等人2016年)。 然而，从增强子生物学和3D染色体接触域的研究中得到的成熟的基因调控原则尚未完全纳入表达性机器学习模型中(Mifsud等人2015;Dekker和Mirny 2016年)。 对更大的序列和不同的实验数据进行建模，为提高预测精度提供了一条前进的道路。 更有效的模型将使研究人员能够描述一种组织或细胞类型的一个实例，并将其投射到具有不同基因组序列的个体。 在这里，我们使用新颖的机器学习算法，仅使用DNA序列作为输入，来学习预测数百种人类细胞类型的数千个表观遗传和转录谱。 通过使用该模型，我们预测了这数千个数据集的基因组变异的两个等位基因之间的差异，特别关注基因表达的预测变化。 我们证明了这一观察在鉴别GWAS基因座内可能的因果变异和机制方面的相当大的潜在价值。
结果
Baenji
Basenji是Basset的改进版，Basset建模基于“峰”的染色质轮廓，尤其关注DNase I超敏感，预测某位点是否开放。
改进： (1)建模远端调控相互作用模型。
(2)预测更精细的分辨率，定量(相对于二进制)基因组谱，更适合于基因表达的动态范围(图1)。
不同：
对count data建模而不是peak data，所以预处理需要做的更~~，
预处理：
需要的数据：
- BigWig coverage tracks
- Genome FASTA file
下载hg19 FASTA文件。hg19.ml.fa、hg19.ml.fa.fai
从FANTOM5获取一些与心脏生物学相关的CAGE数据集（.bw文件）CNhs11760.bw、CNhs12843.bw、CNhs12856.bw。
接下来，我们希望选择基因组序列形成随机梯度下降batch，将它们分成训练/验证/测试集，并构建TFRecords方便后续使用。
最相关的选项是:
在膨胀卷积层之后，每个128 bp的区域对齐到一个向量，该向量考虑了大范围序列中的相关调控元素。最后，我们应用最后的width-one卷积层，对提供的每个数据集的该区域对齐读的归一化计数参数化多任务泊松回归(Hashimoto et al. 2016)。（是不是最后的全连接层预测一个值）也就是说，模型的最终目标是预测长染色体序列128 bp bins的read覆盖率。覆盖率高意思就是对上了？？？
下载使用的数据：
593ENCODE DNase-seq
1704ENCODE histone modification CHIP-seq
356Roadmap DNase-seq
603Roadmap histone modification CHIP-seq
973FANTOM5 CAGE
这些数据是由
529个不同的细胞/组织---->DNase-seq
1136---->CHIP-seq
595---->CAGE
将这些数据进行预处理。PS：这个预处理过程中包含了对multimapping reads（在参考组上有多个匹配点，以前一些组织的测序图谱啥的扔掉了多匹配的reads，使得一般的基因组没有注释）额外的计算，以及对GC bias的归一化（GC bias 是指DNA序列中GC含量过高从而会对某一个研究有影响。比如，测序，使得测序难度增加。因为GC偏好可能会对特定的分析结果造成影响，放大变异影响真实信息，所以需要校正。有一种简单的校正，就是先统计每个GC含量（0, 1, 2, 3,…, 100%）下的特定bin的平均覆盖度，再计算所有bin的平均覆盖度，用来校正测序得到的覆盖度。）。We processed these data with a pipeline that includes additional computation to make use of multimapping reads and to normalize for GC bias(Methods)。 模型结构：
原文中说的是用四个卷积层，每两个卷积层之间有pooling，分别是2、4、4、4，7个膨胀卷积，以及最后一个卷积层预测4229个覆盖的数据，贝叶斯优化参数。
但是在实际实验过程中，作者使用的模型略微有些出入。
预测精度
与一组128-bp bins的held-out test sequences（实验结果）进行比较，结果见A，最左侧分别代表不同方式获得的数据，右边的深色的数据是实验结果，浅色的是预测结果，可以看出大致能预测出peak，但是不同的数据预测的芦苇有些差异，punctate peak data往往更依赖于底层序列（没有太理解底层序列是啥意思？？？），Basenji预测的结果解释了在DNase-seq和CHIP-seq预测活跃调控区域的最大差异。
H3K79me2和H3K4me3预测的结果不是很好，这是因为他们依赖的是较远的序列信号和还没有被完全解释的传播机制。
另外因为噪声（在测序中主要有哪些噪声？？？）的存在，basenji倾向于在低覆盖高估，高覆盖度低估（见C）。
另外和Basset进行了比较：以949DNase-seq，将长的序列分成1024-bp长的序列，利用peak calling 算法将每个短序列进行分类，Basenji获得了更高的AUPRC（在所有的试验中），平均值从0.435--》0.577，中值从0.449--》0.591（这些值不能和以前的Basset论文里面的比较，因为用的数据不一样）。
Basenji预测复制内匹配复制一致性。（Basenji predictions within-replicate match replicate concordance.没有特别理解这个复制的意思）。绘制log-log Pearson相关曲线，然后比较重复实验和实验与预测之间的相关性，平均和中位数交叉重复预测(即重复1对重复2的实验数据的预测)也超过了重复之间的相关性(配对t检验P-value<7 × 10−7)（图D）。P741再好好看一下。
Cell-type–specific gene expression
数据：注释的TSSs（GENCODE v25）
比较方法：在每个TSS旁边128-bp bin的实验值和预测值。
对于训练集之外的基因，对不同的TSS的值算总和来计算accuracy statistics。
做完log2转换后，在不同细胞类型的平均的Pearson值预测的达到>0.85，对于nonzero基因也达到了0.77（中位数是0.86和0.78）（见图3AB）。
CAGE数据集，跟TSSs对齐的reads更多，相关性更强，说明CAGE低测序深度会影响许多模型的精确性。（见A）
在不同类型的细胞预测的也不相同，表明模型学习到了细胞类型特异性转录。为验证实验和预测的一致性，用bootstrap高斯混合模型聚类（见下图），调整后的Rand index分布表明集群存在一致性（见C）。
量化预测特定细胞类型特异性的难度：
在所有CAGE样本中，根据变异系数四等分，计算皮尔逊相关系数。发现（见下图）变异程度一般的相关性还可以，但是四分位数就相关性很差，虽然这些基因突变很多，但是basenji的预测解释了其中相当大的差异，说明学习到了一定的细胞类型特异性。
Distal regulatory elements
量化远端序列如何影响basenji的预测，并且生成基因区域的saliency map（特征图，找出哪些是相对模型重要的特征）.
找到增强子、启动子等~
gene gain：给缺乏某种基因的细胞“添加”上去 gene loss： 将细胞中某种基因进行缺失
验证：在全基因组范围内检测这些调控元件的能力：
数据：ENCODE for GM12878的几个注释:启动子、增强子(不是重叠启动子)和CTCF结合位点(不是重叠启动子或增强子)(ENCODE Project Consortium 2012)。 方法：计算了这些注释的重叠实例的最大显著性分数，并对背景集进行了洗刷。 启动子、增强子和CTCF位点的得分与背景组显著不同(Kolmogorov-Smirnov检验p值分别为9 × 10−64、2 × 10−183和4 × 10−61;方法)(图4 b)。启动子在高端和低端都有更极端的得分，这表明它们通常包含可能调节基因表达率的抑制片段。这种特性也出现在增强剂中，但在压制性终点的幅度要小得多。
Expression QTLs
QTL全称是Quantitative Trait Loci。顾名思义就是基因组上的一些位点，对一些特定性状具有某种量化的影响。这里的重点是，我们关心的性状是一个有多个不同水平的可定量（Quantitative）的性状，而不是定性（Qualitative）的性状。最典型的例子就是疾病（得病/没得病）。一般任何复杂性状（由多个基因决定的性状）都可以认为是数量性状（我其实觉得这里Quantitative翻译成"数量"怪怪的，不过约定俗成了，大家明白就好），例如人的身高/体重/IQ，农作物的株高/产量等。
QTL定位的基本原理，就是测定一群个体的某个数量性状（表型），以及它们的基因型（就是基因组上的一些遗传标记，例如SNP/RFLP等等，但不一定是全基因组），然后寻找基因型和表型的对应关系。 在没有观察GTEx等数据的情况下，经过训练的Basenji模型可以通过比较不同SNP等位基因的模型输出来预测哪些单核苷酸多态性(SNPs)是eQTLs。 基准测试：下载了GTEx V6p版本，并将重点放在与FANTOM5 CAGE匹配的19个组织上。
给定一个SNP-基因对，我们将其SNP表达差异(SED)评分定义为该基因的两个等位基因的TSS(或多个备选TSS之和)预测CAGE覆盖率之间的差异(（注意是预测值之间的差异）见5A)。
LD：可以检测遍布基因组中的大量遗传标记位点，或者候选基因附近的遗传标记来寻找到因为与致病位点距离足够近而表现出与疾病相关的位点，这就是等位基因关联分析或连锁不平衡定位基因的基本思想。
因为连锁不平衡的存在，（可能有另外的位点相关），使得与eQTLs statistics比较比较复杂。
disease-associated loci
deepSEA模型以前对疾病变异有应用，它是预测变异对TF的结合和染色质的影响。
利用CNN，预测ENCODE和Roadmap DNase-seq和CHIp-seq peak calls。
数据集：NIH GWAS Catalog dataset的12296个biallelic（双等位）SNPs
以及一个采样相同大小的带有匹配的小等位基因频率的negative set。（a negative set with matched minor allele frequencies that we sampled down to the same size.）
因为这个实验跟细胞/组织无关，对每个SNP等位基因周围128-bp范围利用模型倒数第二层的预测值。将每个倒数第二层单元的等位基因值之间的最大log2比值分配给SNP。
10折交叉验证
与basenji比较，basenji效果更好
预测出了在哪一段上的基因会影响变异，各种例子，这里不做解释，因为我也看不懂那都是些子啥基因。
discussion
背景：转录调控是跨细胞类型和状态的基因表达特异性的主要驱动因素。基因组研究界需要更有效的模型，了解大型哺乳动物基因组的序列如何决定转录，以便了解基因组改变如何影响这些基因组下游的物理输出。 模型：从DNA序列预测表观遗传和转录谱。 结果：模型解释了这些数据中相当大的差异，包括细胞类型特异性活动。对包含不同版本变异等位基因的序列的预测在统计上与在人类群体中所做的测量结果一致，并接受eQTL分析。
数据集：因为信噪比和测序的技术影响，数据集质量不一。观察到测序深度越深，在replicates之间更大的consistency，模型的预测的精度就越高。
P747里面第一句话讲特定的细胞类型和状态将增强我们彻底检测所有调控元件的能力，并提供调控活动何时何地发生的精确预测，那什么样子的是特定的呢？？？
模型改进的可能：
1、多任务训练的方案可能不是最优的，各个数据集共享信息可能会更好。（我理解的是迁移学习，迁移不同数据集的特征，可能对模型的效果会更好）
2、所有工作的一个大假设是参考基因组指定了所研究的功能注释的DNA（THE reference genome specifies the DNA underlying the functional annotations studied.）Advances that enabled such genome-specific training would likely improve the model 3、虽然用到了膨胀卷积，但是肯定不能把所有的调控元件都捕获到，所以怎么扩大模型的视野对提高模型准确性和解释性会好很多。
尽管只关注转录，而没有考虑转录后调控和调控层之间的相互作用(Skalska等，2017)，我们发现我们的预测信息非常丰富，其中基因组变异将被强调为测量RNA丰度水平的人口研究中的eQTLs。
我们预计在进一步整合功能基因组学训练的调节活性模型与群体基因分型和表型方面有相当大的潜力。
这些正交方法都提供了非编码基因组如何工作的观点，他们的共同考虑应该会使这些观点更加尖锐。
我们认为巴桑吉是朝着这个方向迈出的重要一步
methods
data processing
下载FASTQ格式的文件（FANTOM5得到的973个CAGE实验、593 DNase-seq and 1704 histone modification ChIP-seq performed by ENCODE , and 356 DNase-seq and 603 histone modification ChIP-seq performed by the Epigenomics Roadmap。
用Bowtie 2比对reads，返回最大10个multimapping alignments。
用EM算法将multireads在这10个地方按比例分配。
为了模拟跨GC%谱的趋势，而不排除GC%富集的活性区域，对GC%归一化。
用高斯filter为每个位置分配一个估计的相关的GC%值（为附近的核苷酸分配更大的权重，因为越近越相关）
然后将三次多项式回归拟合到log2的覆盖率的估计值。
然后为了突出显示GC%模型无法解释的剩余覆盖率，重新配置了覆盖率估计。
利用python脚本可以将对齐的BAM文件转化为推测覆盖率的BIGWIG文件。
为了避免组装的间隙以及1kb的不可测区域，所以用131kb的染色体上的不重叠序列。
丢掉了>35%无法映射的序列，剩下14533个序列。
训练：测试：验证=90:5：5
在每个序列中，对128bin的覆盖估计值求和，作为模型要预测的信号。
model architecture and training
利用CNN预测覆盖度
卷积、池化、膨胀卷积、卷积，所有层都batchnormalization、relu、dropout.
标准的卷积层用最大池化到128-bp bins。
一共用了7个膨胀卷积，视野达到32kb。
Glorot 初始化权重，Adam，随机梯度下降，贝叶斯优化 。
1、每隔一段时间反向补充已经完成的DNA序列，并且把序列也反向，2、校序列左右移动0\1\2\3长度来增加数据和减少过拟合。
Peaks binary classification comparison
为了和basset比较，将训练数据和测试数据转换为在更短序列的binary peak calls（131kb---->1024bp）方法：For each dataset, we called peaks on the smoothed, normalized counts in the center 256 bp of the subsequences using a Poisson model parameterized by the maximum of a global and local null lambda similar to the MACS2 approach, and applied a 0.01 FDR cutoff。
然后用已经训练过的最好的参数训练basset。
gene expression cluster comparison
2000个最容易变异的基因，选择了1000个，做10聚类的高斯混合模型聚类，用调整后的Rand index statistic来量化聚类后的簇内的相似性，分布呈现正态分布，这样估计分布的均值和方差，来计算一个分布大于零的P-value。
regulatory element saliency maps
测量远端序列对基因表达预测的影响
（之前有人做过的工作）我们计算experiment-specific saliency maps，作为128 bp的bin表示的点积和模型预测的梯度对序列中这些bin表示的求和。 数据：GM128768启动子(ENCFF 492VIP)、增强子(ENCFF811LAE)和CTCF (ENCFF002COQ)位点的编码注释作为BED文件。从增强子中减去启动子，同时从CTCF位点中减去启动子和增强子。在131-kb的区域内对这些注释进行打乱，从而形成了一个background set，我们为这些区域计算分数并删除与true element重叠的实例。 我们计算了每个saliency score的P-value，一个正态分布调整到background score将有一个更极端值的概率，并使用Benjamini Hochberg false-discovery rate将多个假设修正为0.05。
GTEX eQTL analysis
从GTEX v6版本下载eQTL分析，研究x^2 statistic和significance calls。
因为LD，所以临近的变异有高相关的统计量，basenji将单独的变异分离出来。
计算标记了的LD文件