一次一个染色体：使用Perl，第1部分

2003年9月11日 by

James D. Tisdall是即将发布的 《精通生物信息学Perl》一书的作者。

长期以来，Perl在生物学领域的应用已经成为标准。Perl仍然是生物学家在众多编程任务中首选的语言。Perl之所以能够在系统管理员中成为成功的故事，以及在Web和CGI编程的早期成为一大成功故事，同样的原因也使其成为生物信息学编程的通用语言，被称为生物信息学。

Perl能够同样很好地处理DNA和蛋白质序列数据的原因之一是它非常容易声明和使用字符串。您只需使用它，无需担心内存分配，或者随着字符串的增减管理内存。DNA、蛋白质和其他标准生物数据几乎总是以字符串的形式在Perl中表示，因此这种对字符串的便利直接转化为对DNA和蛋白质的便利。

例如，假设您有一个子程序get_chromosome，它返回一个包含人类染色体中所有DNA的字符串。在人类中，这可能是一个大约100Mb长度的字符串。以下代码片段调用get_chromosome来初始化标量变量$chromosome1，其中包含总结人类染色体1的DNA序列数据的字符串

$chromosome1 = get_chromosome( 1 );

这种编程简单得就像做蛋糕一样。我的意思是，简单得就像馅饼一样。嗯，你知道我的意思。

但在这种奇妙编程的简便性之下，隐藏着一个问题。这是一个可能导致您在处理染色体和基因组时的出色、直观的代码（在您的打印输出中看起来如此优雅，并且看起来被巧妙地分成直观令人满意的、相互作用的子程序）变得低效混乱的问题，当它不是完全堵塞您的计算机时，几乎无法运行。

因此，在这篇简短的文章中，我将向您展示一些技巧，这些技巧可以帮助您编写处理大量生物序列数据的代码——在这种情况下，非常长的字符串——同时仍然从程序中获得令人满意的速度。

内存是瓶颈

问题究竟是什么？通常归结为这一点：通过处理非常大的字符串，每个字符串都占用您计算机用来存储运行程序的主要内存的很大一部分，您可以很容易地超负荷使用可用的主要内存。

当您的计算机上的程序（Linux、Unix或Mac OS X计算机上的一个进程）耗尽主要内存时，其性能开始严重下降。它可能试图通过征用部分磁盘空间来存储它再也无法适应的运行程序的部分来克服快速和高效的主要内存的缺乏。

但是，当一个程序开始写入和读取到硬盘内存时，它可以变得非常慢，而且会很快。根据计算的性质，程序可能会开始“震荡”，即反复在主内存和硬盘之间写入和读取大量数据。您优雅的程序变成了一个贪婪、懒惰的悲观主义者，它抓取了所有可用的资源，然后似乎只是坐在这里。你创造了一个怪物！

上述代码段调用了get_chromosome函数。在不知道该子程序具体功能的情况下，我们可以合理推测它可能从某处获取了100MB的数据，例如磁盘文件、关系型数据库或网站。为了完成这一任务，它至少需要100MB的内存。然后，当它将数据返回并存储在$chromosome1中时，程序还需要再使用100MB的内存。现在，你可能想对染色体进行正则表达式搜索，并将所需的表达式保存到带有括号的特殊变量中，如$1、$&等。这些特殊变量可能非常大，这意味着你的程序将需要更多的内存。

由于这是你编写的优雅、简单的代码，你可能在寻找治愈致命疾病的特效药的艰难过程中多次复制染色体数据或其部分。结果代码可能清晰、易于理解且正确——所有这些都是优秀的代码特性——但字符串复制的数量可能会让你陷入困境。不仅复制大型字符串会占用内存，而且实际的复制过程本身也可能很慢，特别是当复制量大时。

空间效率

更多内容请参阅作者的其他文章

• 使用Perl解析蛋白质结构域
• 生物信息学入门
• 为什么生物学家想编程

当运行程序中有大量数据时，你可能需要在程序设计中添加一个新的约束条件。这个约束是“使用最少的内存”。通常，一个几乎不运行并且占用计算机数小时的数据量，可以通过重新编写算法，使其仅使用一小部分内存来运行几分钟。

这是通过首先减少空间来减少时间的情况。（天体物理学家请注意。）

参考文献

有一种简单的方法可以减少程序中大型字符串的数量。

如果你需要一个子程序返回大型字符串，就像我在示例中使用的get_chromosome子程序一样，你可以使用引用来消除一些内存使用。

将引用传递给子程序的做法对经验丰富的Perl程序员来说很熟悉。在我们的例子中，我们可以重写子程序，使其返回值被放置在作为参数传递的字符串中。但我们不是传递字符串的副本——我们传递字符串的引用，这几乎不占用额外空间，并且仍然允许子程序将整个1号染色体的DNA提供给调用程序。以下是一个示例

load_chromosome( 1, \$chromosome1 );

这个新的子程序有两个参数。假设1将告诉子程序我们想要哪个人类染色体（我们想要最大的染色体，即1号染色体）。

第二个参数是一个标量变量的引用。在子程序内部，引用最可能被用来初始化一个类似这样的参数

my($chromnumber, $chromref) = @_;

然后，DNA数据通过调用它$$chromref被放入字符串中，例如像这样

$$chromref = 'ACGTGTGAACGGA';

不需要返回值。在子程序调用后，主程序会发现$chromosome1的内容已更改，现在由“ACGTGTGAACGGA”组成。（当然，染色体的长度要长得多。）

使用引用也是将大量数据传入子程序而不复制它的好方法。然而，子程序可以更改引用数据的内部内容，这是需要注意的地方。有时你只是想让子程序能够访问数据，但你在子程序查看数据后，仍希望包含数据的变量具有相同的数据。因此，当你将数据引用传递给子程序时，你必须注意你的操作，并确保你知道你的需求。

使用缓冲区管理内存

处理非常长的字符串的最有效方法之一是分批处理。

以下是一个示例程序，用于在整个染色体中搜索特定的12碱基模式，同时使用非常少的内存。（一个碱基是构成DNA的主要构建块之一的四种分子之一。在Perl字符串中，四种碱基用字符A、C、G和T表示。你经常听到生物学家谈论“兆碱基”而不是“兆字节”的字符串。如果你听到那个，你很可能是在和生物信息学家说话。）

当编写一个将在染色体中搜索任何正则表达式的程序时，很难想象你如何避免将整个染色体放入字符串中。但通常，你搜索的内容大小是有限的。在这个程序中，我在寻找12碱基模式“ACGTACGTACGT”。我将从磁盘文件中获取染色体数据。

我的技巧就是一次只读取一或两行染色体数据，搜索模式，然后重复使用内存来读取下一行或两行的数据。

我在编程中需要额外做的额外工作首先是，我需要自己跟踪已经读取了多少数据，以便我可以报告染色体中成功搜索的位置。其次，我需要意识到我的模式可能从一行结束并在下一行的开始处完成，所以我需要确保我在跨越行断点的同时也能在输入文件的数据行内进行搜索。

以下是一个小的程序，它读取命令行上作为参数给出的FASTA文件，并在任何数量的DNA中搜索我的模式——一整个染色体，一整个基因组，甚至所有已知的遗传数据，只是假设数据在命令行中命名的FASTA文件中。我将我的程序命名为find_fragment，假设DNA在名为human.dna的FASTA文件中，我将像这样调用它

[tisdall@coltrane]$ perl find_fragment human.dna

为了测试目的，我创建了一个非常短的FASTA DNA文件，human.dna，它包含

> human dna--Not!  The fragment ACGTACGTACGT appears at positions 10, 40, and 98
AAAAAAAAAACGTACGTACGTCCGCGCGCGCGCGCGCGCACGTACGTACG
TGGGGGGGGGGGGGGGCCCCCCCCCCGGGGGGGGGGGGAAAAAAAAAACG
TACGTACGTTTTTTTTTTTTTTTTTTTTTTTTTTT

以下是程序find_fragment的代码

#!/usr/bin/perl

#
# find_fragment : find 'ACGTACGTACGT' in a very large DNA FASTA file
# using minimal memory
#
#  N.B. This example program does no checking of the input to ensure
#       that it is DNA data in FASTA format; it just assumes that
#       it is. This program also assumes there is just one FASTA
#       record in the input file.
#
#  Copyright (c) 2003 James Tisdall
#

use warnings;
use strict;
use Carp;

# Make sure the program is called with one argument, presumably a
# FASTA file of DNA
my $USAGE = "perl find_fragment file.FASTA";
unless(@ARGV == 1) { croak "$USAGE:$!\n" }

# $fragment: the pattern to search for
# $fraglen:  the length of $fragment
# $buffer:   a buffer to hold the DNA from the input file
# $position: the position of the buffer in the total DNA
my($fragment, $fraglen, $buffer, $position) = ('ACGTACGTACGT', 12, '', 0);

# The first line of a FASTA file is a header and begins with '>'
my $header = <>;

# Get the first line of DNA data, to start the ball rolling
$buffer = <>;
chomp $buffer;

# The remaining lines are DNA data ending with newlines
while(my $newline = <>) {

    # Add the new line to the buffer
    chomp $newline;
    $buffer .= $newline;

    # Search for the DNA fragment, which has a length of 12
    # (Report the character at string position 0 as being at position 1,
    # as usual in biology)
    while($buffer =~ /$fragment/gi) {
        print "Found $fragment at position ", $position + $-[0] + 1, "\n";
    }

    # Reset the position counter (will be true after you reset the buffer, next)
    $position = $position + length($buffer) - $fraglen + 1;

    # Discard the data in the buffer, except for a portion at the end
    # so patterns that appear across line breaks are not missed
    $buffer = substr($buffer, length($buffer) - $fraglen + 1, $fraglen - 1);
}

以下是运行命令perl find_fragment human.dna的输出

Found ACGTACGTACGT at position 10
Found ACGTACGTACGT at position 40
Found ACGTACGTACGT at position 98

代码如何工作

在强烈推荐启用use strict和use warnings，并加载Carp模块以便程序在需要时“崩溃”之后，程序变量被声明和初始化。

FASTA文件的第一行是标题，这里不需要，因此它被读取但未使用。然后读取DNA数据的第一行到缓冲区中，并删除其newline字符。我从这里开始，因为我想即使片段被新行打断，也要搜索片段，所以我至少需要查看前两行；这里我得到第一行，在下面的while循环中，我首先将第二行添加到缓冲区中。

然后是while循环，这是程序的主要工作，它开始读取命令行上命名的FASTA文件中名为的下一行，在这种情况下是FASTA文件human.dna。使用“chomp”删除newline，并将新行添加到缓冲区中。

然后是短while循环，它对缓冲区中的$buffer中的$fragment进行正则表达式模式匹配。它有“g”修饰符，表示*全局搜索（片段可能在缓冲区中出现多次）；和“i”修饰符，表示*不区分大小写搜索，即大写或小写的DNA数据（例如ACGT或acgt）。

当找到片段时，程序简单地打印出位置。 $position 保存缓冲区在总 DNA 中的起始位置，这是我要跟踪的内容。 $-[0] 是一个特殊变量，它给出字符串中最后一个成功匹配模式的偏移量。我还加 1，因为生物学家总是说 DNA 序列的第一个碱基在位置 1，而 Perl 则说字符串的第一个字符在位置 0。所以我将 Perl 的位置加 1 以获得生物学家所使用的位置。

代码的最后两行通过消除其开头部分来重置缓冲区，并相应地调整位置计数器。缓冲区被缩短，只保留输入文件行之间的可能匹配模式的结尾部分。这将是从缓冲区尾部开始的部分，其长度比片段长度短一个碱基。

这样，程序最多在 $buffer 中保留两行 DNA，但仍能搜索整个基因组（或染色体或 FASTA 文件中的任何内容）以查找片段。与读取整个基因组并在这个过程中耗尽内存的程序相比，它运行得非常快。

何时需要麻烦

一个空间效率低下的程序可能在您的较好计算机上运行良好，但当您需要在安装较少主内存的计算机上运行它时，它将完全无法运行。或者，它可能在快速基因组上运行良好，但在人类基因组上运行得非常慢。

一般来说，如果您知道您将处理大量数据集，请在设计和编码时将程序使用的空间量视为一个重要的约束条件。这样，当大量 DNA 被扔向您时，您就不需要回过头来重写整个程序。

编辑注：敬请期待本系列的第二部分，将在本月底发布。在第二部分中，James 将更深入地探讨空间效率，并包括一个使用缓冲区的更通用版本的程序。特别是，第二部分将涵盖使用最小空间运行子例程、完全消除子例程以及具有有限长度的序列基序。

O'Reilly & Associates 将很快发布（2003 年 9 月）《精通 Perl 生物信息学》。

标签

科学
生物信息学
生物序列数据
生物学

浏览他们的文章

反馈

这篇文章有什么问题？请在 GitHub 上打开一个问题或拉取请求来帮助我们。