基因组医学中的机器学习-计算问题与数据集综述

Submitted by chenrouyu on Fri, 05/25/2018 - 16:13
13

一、简介

首先,我们介绍机器学习是如何用来解决基因组医学的关键问题。基因组学研究活细胞的DNA序列所编码的功能与信息结构,而精确医学的出发点就是期望基于患者的所有相关信息(包括基因组信息)来进行对患者进行个性化治疗。这两个领域都在快速发展,尤其是在数据方面[1] - [4]。我们认为上述领域的实际问题给机器学习提供了一个非常好的展示其重要性的机会[5]。在未来,机器学习有可能延长数百万人的寿命并提高他们的生活质量,使得他们避免遭受遗传性疾病或失调症之苦。

一个基因组可以看做是建立一个有机体的说明书。自 1953 年以来,人们已经弄清楚 DNA 分子是遗传信息存储的物理媒介[6],但直到 2001 年人类基因组计划才初步构造了一个典型的人类基因组的原始信息草图[7],[8]。然而,更大的挑战是如何解释遗传信息自身的结构,功能和意义。生物学家埃里克·兰德对当前的研究现状做了一个简单的总结:「基因组像本天书,难读」。不过,人们对于遗传信息是如何被组织成不同的基因已经有了很多认识。每个基因就像说明书中的一章,其描述了如何建立一个特殊的分子系。所谓的蛋白质编码基因[5]描述了如何从氨基酸链(蛋白质)构建大分子,而非编码(noncoding)基因描述了如何从核糖核酸(RNA)链构建小分子。文献[9]和[10]介绍了分子遗传学和细胞生物学。大体来说,人类基因组包含20000个蛋白质编码基因[11]和 25000 非编码基因[12]。一些基因对于生命非常重要、一些则对健康非常重要、而有的即使被剔除掉也无明显的危害性。

一个典型基因的最重要的一个信息结构存在这交替区域,也被称为内含子和外显子。这些区域的边界由核苷酸序列的模式来确定,许多致病突变就是通过破坏这些模式来施加其作用。脊髓性肌萎缩(SMA)是在北美地区导致婴儿死亡的主要遗传原因[13]。如果婴儿的基因组缺失 SMN1 基因,或者该基因被破坏了,那么容易导致运动神经元存活(SMN)蛋白产生不足。该基因的另一种形式,被称为SMN2,可以补偿SMN蛋白的产生。图1显示了来自于蛋白质编码基因 SMN2 的外显子7的核苷酸序列。由于在所显示的四个位置中的核苷酸存在差异,细胞无法识别外显子,从而产生一个不具备正常功能的蛋白,因此也就无法补偿SMN蛋白的产生。研究人员评估了修复外显子7中的SMN2功能的疗法[14],[15]。对脊髓性肌萎缩的研究目前已经非常深入,已经可以通过外在的症状来诊断,但基因检测是确认和研发治疗方案的关键。很多其他的遗传性疾病,致病机制更加复杂。癌症是异质性疾病的一个最具代表性的例子,也即对于同一个疾病,多种病因都会导致相似的症状,但却需要不同的治疗方法[16]。对于癌症,基因组数据正成为研制更多的针对性的诊断和靶向治疗的关键[17]。

8

图1. 利用机器学习算法确定的外显子和调控指令。如果一个婴儿在某个版本的运动神经元存活基因是纯合的,那么会导致脊髓性肌萎缩,这也是婴儿死亡的主要原因。图中位于基因组指令内的三个核苷酸,是被机器学习技术确认是对构建蛋白质时包含这种外显子非常重要。

精确医疗不是一个全新的概念,早在一个多世纪前医生们就已经开始用血型来定制输血[18]。现在不同的是,可以较为廉价地从患者和以及更广泛的人群中收集基因组数据,因此基因组数据在非常迅速的增长,并且各个数据资源之间的共享已经开始在推进。基因组的复杂性和体量是传统实验室中使用数据的20-50倍[17]。本文重点关注基因组医学中的机器学习应用。在基因组医学里,通过评估一个基因组特性来找到靶向治疗,或者从现成的治疗方法中找到比较合适的方法,也可以用来确定疾病风险从而制定相应的预防措施。
我们认为,为了实现基因组医学的目标,必须开发出高效的计算系统,使其能够像细胞里的基因解读机制一样去准确地解读基因组文本。这虽然是一个巨大的挑战,但一旦实现,我们就能够快速地、廉价地去探索遗传变异并制定潜在的有效疗法。这很可能会比使用实验和模型生物等实践方法更加准确。

当前的基因组医学的最新进展是什么?目前,人类基因组中理解得最透彻的区域是蛋白质编码外显子[6]。蛋白质的通用遗传代码在50年前就已经被实验所证实[19],并且探测一个编码突变如何改变相应的氨基酸序列是一个基因组诊断流程的标准特性。例如,如果一个突变引入了一个「终止密码子[7]」进入蛋白质序列(也称为「无意义」突变),那么一般情况下该蛋白质将被截断。然而,突变的预测是否会破坏最终的蛋白质分子的稳定性或结构是一个长期悬而未解的问题[20]。此外,即使有证据表明基因组约 5.5% 的位置经受了纯化选择[21],但编码区只占了人类基因组的约 1.5%。研究发现越来越多的致病突变是在蛋白质编码区之外[22],这说明了只对编码区进分析的工具存在着明显缺陷。很多功能性非编码位置是调控序列,这意味着他们指示着细胞如何来调控重要的过程(例如对基因表达和外显子的可靠识别)。这就凸显出计算模型的重要性,因为计算模型可以自动识别和理解基因组中的调控指令(如图1中的指令)。调控序列显著地增加了细胞生物学的复杂性,这不仅仅归因于基因的绝对数量(实际上香脂杨树的基因总数是人类基因的两倍多[23])或编码区本身(小鼠和狗的基因与人类基因的不同已编码区域只有不到 1%[24])。

我们怎样才能学会自动「解读基因组」?与常见的认知任务如视觉对象检测和语音识别不同,人类先天不具备感知和解读基因组序列的能力,并且先天也不清楚一个活细胞内的所有机制、途径和互动。因此必须开发具备超越人类先天能力的计算系统。下面将会介绍一些知名的机器学习与基因组生物学研究组所开发一些基因组解读的技术。现在逐步达到共识的是,在基因组学与生物学里,越来越多的资源将用于计算新技术的开发而不是纯粹的数据采集----其实一些生物学家对此一直在争论了多年(生物学的「文化包袱」,即授予了超越其他学科的数据特权,在制约着技术的发展[2])。例如,研究人员在癌症基因组图谱(TCGA)项目上花费了近十亿美元后,开始考虑该项目是继续把重点放在测序上还是应该转变到功能分析上[25]。

能够解读基因组文本的计算系统在基因组医学应用上有很多种方式。例如近期「基因编辑」技术的突破允许科学家改变活细胞的基因组,其具体疗效甚至几年前都没有人能预料到。基因治疗现在包括了针对性的基因修改,如删除有害的变异甚至是在基因组的预定位置插入新的序列。基因组编辑技术[26],[27]为基因组医学打开了一扇前所未有的机遇之门。因此在计算机上预测基因组编辑的具体效果也变得前所未有的重要。换言之,知道如何编辑和知道编辑什么是不一样的。

二、机器学习在基因组解读中的应用

从生物标记(如基因组)中预测表型(例如性状和疾病风险)本质上是一个有监督的机器学习问题。其输入是一个与底层生物学相关的DNA序列(基因型)段,输出是表型。但这个过程(如图2(a))对复杂的表型和疾病是不太适用的。原因有两个:首先,一个完整的基因型和表型之间的关联是异常复杂的。即使在单个细胞,基因组一般通过包含多层的复杂和内部链接的生物物理过程和控制机理来调控细胞的状态,这些控制机理已经被进化所重塑。试图仅仅通过对基因组和表型的观察来推断这些复杂的调控过程很像努力让计算机仅仅从二进制代码以及输赢标签而非对弈过程来学习下棋;其次,即使计算模型可以用于推断(也即用以预测疾病的风险),模型的隐变量很可能不符合对应的生物学机理。搞清楚致病机理不光对于研发新的治疗方法特别重要,而且能够为表型筛选提供补充信息。传统上,表型筛选只利用期望得到的生物效果来确认化合物,而没有利用精确的靶位信息[28]。

10

图2. (a)基因组医学的一个主要目标是从基因型(也即某一个生物个体全部基因组合的全体)预测表型,例如疾病风险。(b)利用训练的模型来预测由一系列DNA序列定义的基因型是如何影响细胞变量,如蛋白质浓度等。(b)过程简化并模块化了(a)中的机器学习问题,并有助于探索针对一些关键细胞变量的治疗方法。

我们采用了一个自认为性能更好的方法。这个方法首先训练一个模型来预测可度量的中间层细胞变量(也被称为分子表型),这些变量可以进一步链接到表型(如图2(b))。例如,在前面提到的脊髓性肌萎缩例子里,细胞变量可以是频率。当基因被复制而产生蛋白的时候,外显子由频率来决定其是否包含在内。细胞变量还可以是一个蛋白质与一个DNA的结合位置,或者是一个细胞里的一个基因(转录物)的拷贝数,或者是蛋白质沿转录的分布,以及蛋白的浓度等。细胞变量的更多例子将在下一节进行介绍。

上面提到的方法解决了前面所述的两个问题。由于与表型变量相比,细胞变量与基因型的关系更加紧密,并且更容易根据基因序列来确定,因此学习从基因型到细胞变量的映射模型相对是更简单的。高通量测量技术已经可以产生不同条件下的可用于分析感兴趣细胞变量的大数据。这些数据集能够用来训练更大的、更为精确的模型。此外,由于细胞变量对应于中间层的生化活性数量(如一个基因转录物的浓度),因此它们是非常好的靶向治疗目标。如果与健康人作对比就可以断定出疾病的高风险与一个细胞变量的变化关联关系,那么一个有效的治疗手段就是将该细胞变量恢复至正常状态。另外,在上述的脊髓性肌萎缩例子中,通过修复基因组指令来增加外显子在蛋白中包含的频率这一方法目前正在临床试验中[14]。

本文的余下部分将对机器学习在细胞变量计算模型的构建以及疾病遗传因素的理解上做一个全面的综述。我们提出将细胞变量计算模型的学习过程作为一个中间步骤,并将解释这么做会更有利于处理日益增加的不同类型的数据。我们会详细描述我们已经建模过的两种类型的细胞变量,并简要总结我们是如何理解脊髓性肌萎缩、癌症和孤独症谱系障碍的。为了更好地介绍我们的方法,我们首先介绍现有的疾病风险的计算模型。此外,我们将描述数据集和相关的机器学习问题,以便于数据科学家更好地参与到这个非常重要的领域。

11

图3. 基因组医学中的生物学家、数据科学家、医学家协同工作方式的简单示意图。机器学习的作用主要利用高通量测量数据构建面向细胞变量(对细胞相关功能的一个量化)的专门或者通用的预测模型。通过搞清楚变异是如何通过细胞变量来影响疾病,诊断学家和遗传药理学家能够更加容易发现疾病的直接关联,从而研发新的治疗措施并对患者制定个性化的靶向疗法。

三、细胞生物学、机器学习与基因组医学

这一章将介绍基因组医学涉及到的不同群体的整体工作流程。如图3所示。为了建立一个特定细胞变量的计算模型,首先必须有相应的生物量测量方法,并且还要收集不同条件下的训练数据。20世纪90年代,生物测量技术通常需要多个手动操作步骤并只能获得少量的数据。这样的技术可以用来进行假设检验,但是无法提供足够多的数据来构建准确的具有复杂性输出的预测模型。随着高通量测量技术的商业化,现在能够以较低的成本来获得几十万规模的细胞变量的测量数据。例如,早在数十年前微阵列技术已经被用于窥视活细胞[29],直到今天相关的分析新方法和化学新方法仍在不断的推出,如通用蛋白结合微阵列(PBMs)[30],[31],ChI-芯片[32],[33],和RNAcompete[34],[35]等。高通量测序技术的应用也很广泛[36],包括蛋白结合位点的定位、进化研究中不同生物体的基因组测序、医学研究中的个体基因组分析、以及感兴趣区域或整个基因组变异的发现。

除了大规模基因测序之外,高通量测量技术还可用于细胞变量的测量,例如不同转录物的丰度[8][37]。尽管对于癌症和某些神经系统疾病,体细胞突变在受孕后可能会在DNA里面发生改变[38],[39],一个个体的基因组是相对稳定的。另外一方面,不同细胞的转录组[9]是有变化的,并且还受细胞周围环境(如它代表的组织)的影响。以前,微阵列被用来测量大规模转录物,但现在通常选择高通量测序方法。另一个高通量测序的应用是分析蛋白质如何与DNA [40]的特定区域进行交互。蛋白质的结合可以影响基因组的指令是如何执行的,这展示了细胞生物学调控的复杂性。这些所关注的特定细胞变量的测量数据让我们可以在一个生命体所定义的指令的最基本层面来理解细胞的基本工作原理。高通量测量技术能够涵盖大部分基因组的各种不同细胞状态下的细胞变量,包括疾病状况。很多这方面的数据现在是公开的。这给数据科学家带来了一个难得的机会,可以使用机器学习技术来构造细胞变量的预测模型。

预测模型的输入包括从一段DNA得到的特征,例如特定核苷酸的频率或某些模体[10]的存在,其中一些特征可以从DNA序列本身得到[41]。为了解释DNA里面的编码指令是如何通过生化过程和结构来对细胞变量产生影响,需要用到一些其他的特征,例如,蛋白质结合的DNA和RNA,核小体定位和占有情况[42],以及RNA的二级结构[43]。通常,良好的模型输入应该是可以直接从DNA序列中提取到的。而对于一个计算模型是否能够在基因组医学上发挥作用,很多程度上就是取决于模型输入是否容易得到。鉴于全基因组测序的成本在持续下降,越来越多的基因组数据将可用于模型训练,并在基因组医学的应用背景中,它可能会成为获取病人基因组的标准工具。

图3所示的过程中的一个重要方面就是利用机器学习得到的模型能够泛化(或推广)到新的遗传背景下。例如,我们可以利用公开的基因组和健康组织提供的数据记录来构建计算模型,然后将其应用到患病细胞的基因组上去判定患病细胞的变化分布。「泛化」是模型构建时的一个非常关键的性能指标。从建模的角度来看,所构建的计算模型必须具备强大的泛化能力,以便能够在新的遗传背景下处理细胞状态数据。因此,模型构建的一个重要方面就是在验证阶段使用模型从未见过的DNA序列以及那些在训练的时候没有碰到过的细胞状态数据。很难构建一个对与训练数据完全不同的DNA序列与细胞状态数据都能够比较准确处理的模型,所以在验证阶段,我们必须要挑选有代表性的测试样例来更好地评价模型的可靠性。

如果一个模型具有很好的泛化性能,那么它能够处理导致细胞变量变化的变异DNA序列。这或许能够在不需要对患病细胞进行检测的条件下预测疾病状态。在实践中,基于一个在参照的基因组和正常组织数据上训练好的模型,零样本学习已经成功应用于确定导致不同疾病的基因变异。

当用于预测细胞变量的模型没有直接考虑到疾病的相关信息时,如果模型准确地反映了指令在基因组中是如何执行的,那么它应该能够找到改变细胞变量的基因变异所导致的疾病。这种方法已被证明能够较好地处理大量的突变与疾病[44],不过它也会出现错误。如果利用细胞变量改变的程度来对基因变异进行打分,那么当一个非致病性变异极大地改变了细胞变量,就会出现假阳性。例如,一个改变细胞变量从而导致头发颜色变化的基因突变就是假阳性。而有些变异所改变的细胞变量没有在模型中考虑的,假阴性就会出现。上述两种错误都是由于模型的不完全精确造成的。当研究特定疾病的时候,对变异分值的计算可以将疾病相关的数据(如人群数据)纳入进来。利用这种方式就可以从候选变异集里面挑选出对细胞变量最可能有影响的部分。更一般地,这些分值可以作为某些疾病专门模型的输入特征,很多基因组区域的分值都可以用到模型里面。

在总结完基因表达的过程后,我们提供几个具体的例子。重点是剪接和蛋白质核酸结合,但与疾病相关的细胞变量千变万化,如转录率[45],DNA甲基化[46],[47],聚腺苷酸化[48],染色质结构[49],[50],RNA折叠[51],和蛋白质折叠[52]。

12

图4.基因表达有三个主要步骤构成。首先,基因转录(transcription)制造了一个RNA分子(本质上是DNA的拷贝)。在这个步骤里,RNA分子称为前体信使RNA(前mRNA)。RNA过程然后修改前mRNA。这个过程包括剪接掉长的序列片段(称为内显子)并连接侧面区域(称为外显子)。在这个步骤里,RNA分子称为信使RNA(mRNA)。基因转录通过读取mRNA序列的三个字母编码制造了蛋白质分子(一个氨基酸链)。其他的过程包括多聚腺苷酸化(将腺嘌呤基扩展到mRNA的尾部)、mRNA稳定化(对mRNA分子进行处理使得它不容易被降级)、mRNA定位(mRNA被移到一个容易转录的位置)、以及蛋白质定位(将蛋白质迁移到细胞里面的一个特殊类型的位置)。

13

图5. 基因在选择性剪接过程中能产生不同的蛋白质,这里基因指令根据细胞的情况(如细胞类型)来决定一个外显子是否包含还是排除。利用RNA-seq,特定细胞类型包含的成千上万的每一个外显子的频率都能够被测量到。这些数据能够用来训练一个可以发现控制剪接指令并且预测剪接的计算模型。

A 基因表达

在基因表达里,基因首先被复制(转录)来得到一个信使RNA(mRNA),然后mRNA被翻译得到一个蛋白质。含有外显子和内含子的DNA序列首先转录成RNA,这被称为前体的mRNA(前mRNA),如图4所示。「前体」是指前mRNA需要在细胞核内进行进一步的处理以得到一个成熟的mRNA。在RNA处理的过程中,会发生各种修改,其中之一是剪接[53]。剪接移除来自前mRNA的内含子并与外显子连接在一起。另一步骤是聚腺苷酸化来将一个腺嘌呤碱基的序列扩展到mRNA的尾部[48]。在标准模型中,剪接移除内含子并保留所有外显子,如图4所示,但大多数基因可能以不同的方式被剪接,以便外显子有时去除和/或内含子被保留,这会增加蛋白质的种类。剪接是一个关键的细胞过程,已被集成在基因调控网络里[54]。剪接后可能会出现在转录完成后,但是它经常与转录同时发生,从而使转录和剪接过程相互作用[55]。最后,mRNA转运出细胞核至核糖体,可以将mRNA转换成成蛋白质。这些是基因表达的主要过程,其余也包括mRNA的稳定和蛋白的定位。

B 剪接的计算模型

尽管能够产生高度复杂的组织和器官,除了功能分子如微RNA,人类基因组仅包含2万多个基因[11]。人类基因组实现复杂性功能的一种方法是通过指令来使用不同的方式去剪接同样的基因。在选择性剪接里,一些外显子可在剪接过程中移除,这取决于细胞的上下文(如细胞所在的组织)以及附近的基因组指令。这个过程的控制称为「剪接调控」,它取决于在DNA和前mRNA中的许多基因元素的相互作用,这是调控点附近的特征,以及蛋白质或者其它与基因元素交互的分子等反式因子[56]。鉴于一个基因内有多个外显子,选择性剪接可以以不同方式来连接外显子,并产生不同的蛋白质亚型(如图5上部),从而「扩大」细胞的蛋白质指令[57]。最新统计数据显示,平均每个蛋白编码基因的转录物数量大约是4 [11]。在显著情况下,白细胞介素7受体(IL7R)的外显子的选择性剪接能够得到转化的水溶性或水不溶性蛋白质和局部膜[58]。选择性剪接的重要性已经得到了证据的支持,人类基因组中至少95%多外显子基因是通过选择性剪接得到的[59]。许多已发现的人类疾病是受到规则剪接指令的变异影响的[60]。

选择性剪接在更复杂的细胞里出现得更加频繁。在包括不同的脊椎动物物种发现的所有细胞里,人类大脑细胞具有最复杂的选择性剪接模式[61]。自然地,异常剪接被认为是一些精神疾病主要诱因。例如转录组分析发现在泛自闭症障碍(ASD)的皮层区域里的选择性剪接模式有一致性偏差[62],[63],这说明了剪接失调是一种泛自闭症障碍机制。

Lo'pez-BIGAS等人估计,多达60%的遗传疾病是由与剪接过程存在的缺陷导致的突变引起的[64]。特别是,近三分之一的致病突变会改变剪接位点[65],很多会导致异常外显子跳读[11][66],[67]。此外,近45%的疾病/性状相关的变种驻留在内含子,其中大部分被认为是调节剪接模式[22]。反常的剪接主要在两个方面导致异常:1)它导致失活或更低效的蛋白亚型;或2)它破坏蛋白亚型的平衡。和剪接相关的主要疾病包括神经和精神疾病,囊性纤维化,帕金森病,脊髓性肌萎缩,肌强直性营养不良,肌萎缩性侧索硬化症,过早老化,以及数十种癌症[56],[66],[68]。

在过去的十年中,我们的团队研发了精确的剪接计算模型,还提出了剪接生物学机制的新观点、并确定了不同疾病和神经系统疾病的遗传原因。通过精确的剪接模型和AS计算,我们已经能够预测剪接是如何受基因变异的影响,并以此评估一个基因组中的突变是否会影响患病风险。剪接的计算模型的输入是基因组的剪接发生的邻近区域特征。计算模型预测对应的外显子被mRNA保留或排除的频率[44],[69] - [71]。为了训练模型,RNA测序(RNA-SEQ技术用于确定哪些亚型存在于不同的组织类型的细胞里。下面对RNA-SEQ做一个简要的描述,细节可以在[72]中找到。RNA-SEQ涉及了所研究的细胞群体的mRNA的第一次分裂,然后对这些分裂进行序列化。RNA-SEQ片段,或称为「reads」,随后被映射到生物体的参考基因组上,这个映射允许出现少量的误匹配。「reads」的数量被用于确定样品中的同类型测序片段的相对丰度[73]。对于剪接,被映射到外显子的边界之间的连接处的「reads」将用于分析[44],[74]。这能够计算出一个外显子包括或排除在mRNA内外的频率(图5,底部)。

最近,我们展示了如何将上述方法应用于发现变异对剪接的新的影响方式以及如何导致疾病的[44]。为了分析一个特定的遗传性变体,我们将正常的DNA序列和突变的DNA序列输入到计算模型,把输出在剪接水平上的变化作为突变是否可能是有害的指标。我们的方法能多对临床变异和一系列疾病的合成变异的效果进行精细的预测,包括脊髓性肌萎缩,遗传性非息肉性大癌症和非自闭症障碍等。实验表明,模型的预测与实验数据相当匹配[44]。通过分析泛自闭症障碍患者所有的罕见基因组变异,熊等人确定了在中枢神经系统发展、突触传递,神经元投影中的19种反常剪接的基因[44],这显著扩大了公认的与泛自闭症障碍有关的基因集。

这个例子阐述了图3所示的过程。首先在实验室收集感兴趣的细胞变量的数据,然后传给数据科学家来构建模型去分析遗传病。

C 蛋白质DNA和蛋白质RNA结合的计算模型

一个「蛋白质序列相互作用」是一种蛋白质分子和DNA或RNA链之间的化学吸引力,可能发生在活细胞(在体内),也可能发生在一个受控培养基(在体外)。如果具有很多训练数据,那么蛋白质序列交互作用的强度是细胞变量一个低层次的分类。这些交互作用对于精确建模非常重要,因为他们影响到细胞里面的许多关键的流程。DNA和RNA结合基因分别选择性地结合到DNA或RNA链上,并且对结合点的突变非常敏感。目前已经知道很大一类由改变结合点的突变导致的疾病,包括胆固醇产生过剩[75],黑素瘤[76],和前列腺癌[77][参见图 5(H)]。蛋白质序列结合的精确模型对于解释基因组和预测突变的影响是非常重要的,机器学习在这些地方可以发挥重要的作用。

人类基因组编码了至少1400种DNA结合蛋白[78]和至少 1500 种RNA结合蛋白(RBPs)[79],使得这些成为最大类别的蛋白质。许多DNA结合蛋白被称为转录因子(TF),因为当它们结合时,它们影响特定基因转录到RNA的速率。RNA结合蛋白可以影响RNA的后续处理,例如,通过折叠该RNA[80]或支配某些外显子的mRNA剪接方式[81]。实际上,RBPs的计算模型提供了我们上面提到的剪接模型的关键输入特征。图6(顶部)只是描绘了这些蛋白质DNA和蛋白质RNA相互作用从而影响细胞状态的两个重要方面。

生物学家已经开发了高通量技术来测量个体蛋白质的特异性序列。这些测量可以用于训练蛋白质序列交互的计算模型:输入是基因序列,输出则是结合值[82]。但是,在训练预测模型前,最好理解这些测量的表达含义以及他们是怎么获取的,可能会存在哪些偏差。基于序列的方法和基于微矩阵的方法是两个流行的高通量实验方法。面向这两种方法的大数据集都存在。下面首先回顾一下两个方法的过程,然后介绍流行的计算模型。其他的一些蛋白质DNA测量技术,请参加Levo和Segal的论文[83]。

14

图6. DNA-和RNA-结合蛋白调控细胞里的很多过程,包括基因转录率(左上部)和剪接率(右上部)。为了判定特定蛋白结合的模体,生物学家通过对活细胞的蛋白质边界片段进行测序(中部)或者通过暴露标注的蛋白质到一个微阵列(底部)上的合成片段来获取特异性数据。

1)测序:测序方法的工作原理是首先分离出与感兴趣蛋白质结合的DNA或RNA片段(长度约100),然后对这些片段进行序列化,最后将其映射至一个参考的基因组。其基本思想是,无论映射的「reads」在哪里重叠并堆积形成一个「尖峰」,我们就可以推断,该蛋白质偏好于在尖峰(图6中)附近进行结合。一个常见的DNA结合-蛋白质实验方案是ChIP-seq[40],对于RNA结合-蛋白则有RIP-和CLIP-SEQ方案[84]。每个测序实验会产生100-100000个所有的从参考基因组得到的短序列(每一个对应一个不同的峰值)。机器研究者会将这些「尖峰」标记为弱标签,因为外观和结合点的位置是未知的,很像ImageNet分类的挑战[85]。从序列的这个列表中,我们可以学习到一个模型来区分蛋白质结合序列和蛋白质很可能忽略的序列。计算生物学家开发了一些工具用于在一个尖峰集合里面查找重复的模体。一个流行的软件工具是MEME-ChIP[86],它使用期望最大化方法来从一个背景序列(隐形)的统计模型里面找到一个或者多个尖峰(阳性)。该工具自20世纪90年代开始一直沿用至今。

测序的一个难点是所得到的尖峰通常包含了与实验中的目标蛋白质的结合位置邻域相同的其他蛋白质的模体。对于计算机视觉研究人员,这类似一个无监督的对来自于真实文本图像的字符‘q’进行建模,但发现最后的模型对字符‘u’的响应也非常强。那么,我们不应该建立只能识别一个特定的手写字符的模型,而应该去建立一个能够识别一个特定字符趋于出现的上下文的模型。而生物学家关心的是正确的模体应与正确的蛋白质相关联,因此计算模型的可解释性是一个非常重要的方法。目前,生物学家完全依靠人工来核查MEME芯片这样工具的输出,可视化每个发现的模体,然后基于背景知识来推理出哪些模体可能属于靶蛋白质。

图7. Ray等人构建的RNA-结合蛋白RBFOX1的位置频率矩阵(PFM)[35]. 一个PFM能够用来沿着一个目标序列对每个结合位置进行打分(图右)。生物学家通常用PFMs来可视化一个序列。

2)微阵列(基因芯片):高通量蛋白结合芯片(PBM)实验将感兴趣的蛋白质暴露至40000-250000个不同的「探针」序列[87]。它芯片最初始不知道蛋白质在从哪个或哪些范围的探针序列开始结合的。该实验依赖于一个特制的滑板,其表面是含有微小凹槽的矩阵---每个凹槽对应着一种类型的蛋白质。蛋白质的拷贝通过滑动来引入,以便它们可以有选择地与探针结合。在此之后,未结合的蛋白被冲走。然后利用一个荧光团的激光来激发未被冲走的蛋白质,使其发射可见光。整个滑动过程被拍摄下来,并且测量每个凹槽的亮度来反应每个凹槽的蛋白质浓度(见图6,底部)。一个PBM实验的最终结果是一系列序列-得分对数据,这些数据可以用来训练预测蛋白质序列特异性的计算模型。

探针序列被明确设计成能够覆盖住所有可能的长度为k的子序列。k通常选择为所涉及的高通量阵列长度的约10倍。探针组从一个k等级的de Brujin序列开始产生,然后被重叠性地分割成期望的探针序列,每个序列的长度约为40[35]。因此,与测序数据不同,微阵列数据不存在无关模体频繁出现在阳性集这种情况。原始的亮度测量微阵列实验仍然可能有噪声、非自然数据和偏差。因此亮度测量值需要进行进一步后处理,包括异常值去除,空间趋势剔除和分位数标准化等(见[35],补充资料)。

3)基本计算单元:早期的对结合点建模的方式是通过「一致序列[12]」来实现,其中假定蛋白质偏好于以一种固定的模式来选择结合点位置,并同时允许有时候会出现一个或两个错配的位置。如今,对结合点位置进行建模的主流方法是位置频率矩阵(PFM)以及在此基础上的若干变种[88]。一个蛋白质的PFM的参数(i,j)代表了当PFM分配该基因结合点的时候,第i个基出现在第j个位置的频率[89]。每一列的频率归一化为概率。例如,图7(左)显示了一个用于人体RBFOX1蛋白质建模的PFM。例如,RBFOX1是一个影响神经元发育过程中的剪接的RBP (图6,上部),并且是与患有自闭症[90]有关。论文[35]中显示的PFM参数是从一个基于微阵列实验的高密度探针所获取的重复模体数据中推断得到。RBFOX1的一致序列普遍认为是gcaug,但[35]给出的PFM表明,它更偏好于结合到较长的模体ugcaug上,这称为是二次偏好,有时也写作为(U)gcaug。一旦一个蛋白质的PFM确定之后,它可以充当用于沿着基因组的模板匹配的软模式。给定一个序列s = {s1, …, sn}并且其中si的取值范围为{1, 2, 3, 4}(分别代表a,c,g,t),那么一个包含4×m个参数的矩阵M的归一化PFM的结合值zi为

15

当对每个位置进行打分之后,整个序列就可以用不同位置对应分值的均值、最大或者分值和来作为整个序列的分值。对于熟悉深度神经网络的研究者来说,PFMs扮演了一个类似于卷积神经网络中滤波器的功能[91], 这里输入可以看做是一个1维图像每个通道表示{a, c, g, t}是否出现。如图7(右)所示。

在结合复杂度上,图7显示的RBFOX1的例子是一个最容易用PFMs建模的蛋白质之一。许多蛋白质的序列特异性是通过2个PFMs的组合来进行更精确的建模,其中次级PFMs用来描述结合的次要模式[92]。使用默认设置,MEME-ChIP软件套件[86]可以输出最多三个不同的PFMs列表,其中每一个相对于MEME的建模假设是有显著统计意义的。MatrixREDUCE软件套件[93]可输出高达20个类PFM的矩阵,被称为面向特定位置的仿射矩阵(PSAMs)。利用PFMs和PSAMs进行建模有一个难点就是模体的长度(矩阵的列数)必须事先通过人为指定。而这个指定可能影响结果的质量;通常的做法是设置一个带有余量的长度值,然后在矩阵的后处理步骤中去掉高熵列(即「不关心」位置)。

4)其它的非基于PFM的模型:PFMs最明显的缺点是其基于了一个很强的假设,也即每个位置对于结合强度的贡献是独立的。这个假设对于一些蛋白质大体是成立的,但对许多来说并非如此。计算生物学家积极开发了新的模型来捕捉更具挑战性的蛋白质特异性。例如有学者开发的模型可能获得剪接位点的全局特征[94],有的模型则专门用于建模某些特定的蛋白家族如C2H2型锌指状物[95],[96],还有关于「结构选择」的RBPs模型,它对RNA自身折叠的方式很敏感[97]。蛋白质也可以结合在链的邻近位置,从而导致蛋白质-蛋白质相互作用进一步增加了剪接特异性的复杂度[83]。一些更为典型的模型通常依赖于生物学知识,如蛋白质的晶体结构或关于蛋白质表面上结合的领域知识。Stormo和Zhao[98]提供的一个关于DNA-结合蛋白实验与基于PFM计算模型的综述。论文[84]则对于近期的RBP实现和模型给出了一个综述。Slattery等[99]提供一个关于TF-DNA结合的更多深入问题的精彩回顾。

对于有兴趣研究序列特异性问题的机器学习研究者,一个很好的入手点是DREAM5 TF-DNA模体识别挑战赛[100]。在这个挑战赛里,Weirauch等人提供了PBM和ChIP-seq数据并评价了26个计算生物学领域的专门算法的性能,包括MEM-ChIP和MatrixREDUCE。 Alipanahi等[41]最近利用卷积神经网络来参与这一挑战。在PBM和Chip-seq数据相关的两个挑战赛上,他们得到的结果要优于原来的算法。DREAM5挑战赛只是一系列大型计算生物学挑战赛的一部分。详情请见[101][102]。