乘风原创程序

  • Perl 如何拯救了人类基因计划
  • 2008/5/24 14:54:05
  • 原始出处: http://www.tpj.com/
    本文作者: Lincoln Steinx
    本文译者: tcwu

    Reprinted courtesy of the Perl Journal, http://www.tpj.com/
    Lincoln Stein's website is http://stein.cshl.org/

    位置:英国剑桥,欧洲最大的 DNA 定序研究中心的会议室

    场合:本中心与美国最大的 DNA 定序研究中心的电脑科学家的高阶层会议

    问题:虽然这两个中心使用几乎相同的实验技术,几乎相同的数据库,还有几乎相同的资料分析工具,却仍然无法交换或者比较出有意义的结果。

    解决方案:Perl

    人类基因计划在大约 1990 开始雄心勃勃的开始发展,期望能以国际合作的力量找出人类以及部分实验动物的完整 DNA 序列。这项工作的性质兼具科学与医药。藉由了解生物起源的构造,以及种种令人苦恼的细节,我们希望能够了解有机体如何由一颗卵发展到复杂的多细胞生物,食物如何代谢以及被转换成身体的一部分,以及神经系统如何流畅的工作。从医学的角度来看,从全盘了解完整的 DNA 序列得来的知道,能够有效的加速找出人类疾病(以及可能的疾病)的原因。

    经过六年的努力,基因计划已经超越原先的计划日程了。人类以及实验动物的详细基因地图已经完成了(用一连串的标记来安排 DNA 是决定完整序列的第一步)。最小的有机体,酵母,已经近乎完成,接著下一个,细小的蠕虫也不远了。几个月以前大规模的人类 DAN 定序工作已经在各中心展开,并且在整年中都会全力进行。

    从资讯处理的角度来看,DNA 是一个非常长的字串且由四个字母构成,分别是 A,T,G,C(这四个字母分别为四种化学物质的缩写,是构成双股螺旋的基础。字串的长度令人印象深刻但也不是让人特别惊讶: 3 x 10^9 的字母长,或者如果你用一个位元来储存每个字母的话,需要 3GB 的储存空间。

    3GB 不小但是以今天的标准无疑的可以被管理。不幸的是,这只是储存结果所需要的空间。要得到这些结果所需要的储存空间远远大的多。根本的问题在於,目前的 DNA 定序技术一次最多只能读到 500 个连续的字母。为了决定更长的序列,必须将 DNA 定序为一小块一小块部分重叠的片段,称为”reads”,这些拼图般的区域用演算法来检查并找出相配的地方。因为 DNA 序列并非随机的(相似但并非完全一样的基本图案在基因中出现很多次),而且因为 DNA 定序技术有杂讯且倾向於错误,必须针对某个区域作 5 到 10 次,才能可靠的把”reads”片段重组成真正的序列,这增加了要管理的资料的数量。在这之上的是与实验室工作相关的资讯: 谁进行了实验,什麽时候完成的,被定序的基因区段,用哪个软件以及哪个版本重组序列,对实验加上的注解,等等。除此之外,一般人通常想要将定序用的机器产生的原始资料储存下来,每 500 个字母会产生 20-30kilobytes 长度的资料档案!

    这还并非全部。只决定 DNA 序列是不够的。在这些序列里面,有功能的区域散布在很长的没有功\能的区域中。人类 DNA 中有基因,控制区域,结构区域,甚至有些病毒卷入人体中并且变成了化石遗迹。因为基因跟控制区域负责健康与疾病,研究者会想要标示以其注解它们当它们被重组的时候。这类的注解产生了更多的资料需要被追踪与储存。

    某些人评估居会有 1 至 10TB 的资讯需要被储存才能看出人类基因计划的结论。

    所以 Perl 能做什麽呢?从一开始研究者就了解到资讯学将会在基因计划中扮演很重要的角色。一个整合每个基因中心的资讯学核心被建立了。这些核心的任务有两个: 提供电脑支援和数据库福对给他们相关的实验室,还有发展资料分析和管理软件给整个基因研究社群。

    公平的来说,资讯科学小组一开始的结果是好坏参杂的。有些小组试图在复杂的关联数据库上面建立大型的系统,它们一次又一次的受到生物研究的高度动态的阻挠。当一套系统里里外外运作都与复杂的实验室协定正常配合,被实作出来且经过除错,新的科技已经取代了旧的协定,於是软件工程师又得回到设计版前。

    然而,大多数的小组学会了建构模组化的,松散结合的系统,这些系统的部分可以被拿出来或放进去而不需要重新更换整套系统。举例来说,在我的小组里面,我们发现许多资料分析工作牵涉到一连串的半独立的步骤,假想某人想要操作一个刚刚被定只好的位元(图一)。首先需要一个基本的品质测试:序列够不够长 ? 含糊不清的字元是不是在最大限制以下? 接著有”向量检查”。为了技术上的理由,人类 DNA 在被定序前必须经过细菌处理(这就是”复制”的程序)。并不常发生,但是人类 DNA 有时候会在处理的过程中遗失,而且整个序列包含了细菌的向量。向量检查确保只有人类的基因被记载到数据库里面。接下来是反覆序列的检查。人类 DNA 充满了重复的元素,使得将拼图安装在一起变的富有挑战性。重复序列测试试图由一个已知重复序列的数据库找出新序列,倒数第二个测试则是试图找出新的序列靠著对照於一个很大的社群的 DNA 序列数据库。在这时候一个比对成功通常能提供线索找到新序列的功\能。完成这些测试之後,序列以及它的资讯已经被收集且载入了实验室的区域性数据库。

    将 DNA 序列通过这些独立分析步骤的过程看起来很像一根水管,因此不用多久我们就了解 Unix 系统上面的”pipe”可以操作项工作。我们发展了一个简单的以 Perl 为基础的资料交换格式,叫做”boulderio”,这格式允许松散连结的程序加入资讯到以导管(pipe)为基础的输入/输入字串流。 Boulderio 的基础是一对一对标签/值。Perl 模组让获得输入很容易,取出他们有兴趣的标签,在上面作一些事情,然後丢回输出字串流。所有该程序不感兴趣的标签被传到标准输出,所以其他的程序可以使用它们。

    在这种形式的设计下,分析一个新 DNA 序列看起来像是这样(这并非我们真正在用的程序源代码,但已经很接近了)


    name_sequence.pl < new.dna |
    quality_check.pl |
    vector_check.pl |         
    find_repeats.pl |
    search_big_database.pl |
    load_lab_database.pl  

    1 一个包含新的 DNA 序列的档案会一个叫做 "name_sequence.pl" 的 perl script 处理,该程序的任务只是给予这个序列一个新的且唯一的名字并且将他输出成 boulder 格式,它的输出看起来像:


    NAME=L26P93.2    
    SEQUENCE=GATTTCAGAGTCCCAGATTTCCCCCAGGGGGTTTCCAGAGAGCCC......  

    来自 name_sequence.pl 的输入接著被传入品质检查程序,它会找寻 SEQUENCE 标签,跑品质检查演算法,然後把它的结论写到资料字串流里面。资料字串流现在看起来像:


    NAME=L26P93.2    
    SEQUENCE=GATTTCAGAGTCCCAGATTTCCCCCAGGGGGTTTCCAGAGAGCCC......    
    QUALITY_CHECK=OK  

    现在资料字串流进入了向量检查。它会从字串流中取出 SEQUENCE 标签并且执行向量检查演算法。资料字串流现在看起来:


    NAME=L26P93.2
    SEQUENCE=GATTTCAGAGTCCCAGATTTCCCCCAGGGGGTTTCCAGAGAGCCC......
    QUALITY_CHECK=OK
    VECTOR_CHECK=OK
    VECTOR_START=10
    VECTOR_LENGTH=300

    这些东西继续通过导管(pipeline),直到最後 "load_lab_database.pl" 程序核对所有收集的资料,将一些最後的结论是否适合未来适合标记起来,以及把所有的结果送入实验室的数据库。Boulderio 格式的一个很好的特性就是许多个序列的纪录可以被连续的在 Unix 导管里面被处理。用一个”=”表示表示一笔记录的结尾,以及下笔资料的开始:


    NAME=L26P93.2
    SEQUENCE=GATTTCAGAGTCCCAGATTTCCCCCAGGGGGTTTCCAGAGAGCCC......
    =
    NAME=L26P93.3
    SEQUENCE=CCCCTAGAGAGAGAGAGCCGAGTTCAAAGTCAAAACCCATTCTCTCTCCTC...
    =  

    也可以在纪录中再建立子纪录,这使得我们可以使用结构化的资料型态。

    以下为范例源代码,这个源代码处理 boulderio 格式。它有物件导向的风格,纪录被从输入字串流中拉出来,修改,然後丢回字串流中。


    use Boulder::Stream;
    $stream = new Boulder::Stream;    
    while ($record = $stream->read_record('NAME',