省得自己写代码，循环、遍历，纠结于各种算法的效率上。

直接拿来主义，用已有的模块，List::Compare。

既然发在CPAN上，估计是优化过的，而且反正perl的速度本来就不快，又只是自己用，简单易用省心，最好。

1
2
3
4
5
6
7
8
9

use List::Compare; #需要安装
@Llist = qw(abel abel baker camera delta edward fargo golfer);
@Rlist = qw(baker camera delta delta edward fargo golfer hilton);
$lc = List::Compare->new(\@Llist, \@Rlist); #使用数组的reference

@intersection = $lc->get_intersection; #子集
@union = $lc->get_union; #并集

print "@intersection\n@union\n";

List::Compare还有很多其他用法。

比如，get_unique()，返回第一个数组特有的那些元素；

而get_complement()，返回第二个数组特有的那些元素；

以及get_symmetric_difference()，返回属于某一数组但不同时属于两个数组的那些元素。

详见List::Compare。

Steve Jobs (1955-2011), “stay hungry, stay foolish.”

too busy in these days.

用惯了apt-get和yum，现在转到mac上，总想着sudo apt-get install一下安装某些软件。其实mac也有类似的软件，MacPorts和fink。

今天先介绍MacPorts。

先安装xcode和X11，然后到这里下载.dmg。安装以后就可以在终端里类似apt-get那样安装软件了（包括GNU软件、开源软件等），如：

sudo port install wget

sudo port selfupdate       #升级port
sudo port sync             #同apt-get的update。
port list                  #列出所有软件
port search XXX            #查找XXX软件
port deps XXX              #查看XXX软件的依赖
sudo port install XXX      #安装XXX软件
sudo port uninstall XXX    #卸载

更多的细节可以在MacPorts官网获得。

输出序列文件

下面将展示如何利用两个模块创建一个序列文件。上例中，已经有了一个序列对象$seq_obj，然后需要创建另一个用于读写文件的对象，SeqIO对象。IO表示Input-Output（输入输出）。Bio::SeqIO可用于读取和输出各种Bioperl支持的序列格式文件（支持的序列格式列表详见SeqIO HOWTO）。创建Bio::SeqIO对象和前述使用new()创建序列对象的方法类似，如下所示：

1
2

use Bio::SeqIO;
$seqio_obj = Bio::SeqIO->new(-file => '>sequence.fasta', -format => 'fasta' );

注：在-file参数中，“>”符号表示要创建一个名字为“sequence.fasta”的文件用于内容输出。这和一般Perl脚本中，使用函数open()来写文件也是用“>”。（译者注：< 输入、> 输出）。在“-format”参数中声明序列格式是“fasta”，所以创建的会是一个fasta格式的序列文件。

现在来把刚才的两个例子放在一起：

1
2
3
4
5
6
7
8
9

#!/bin/perl -w
use Bio::Seq;
use Bio::SeqIO;
$seq_obj = Bio::Seq->new(-seq => "aaaatgggggggggggccccgtt",
                                           -display_id => "#12345",
                                           -desc => "example 1",
                                           -alphabet => "dna" );
$seqio_obj = Bio::SeqIO->new(-file => '>sequence.fasta', -format => 'fasta' );
$seqio_obj->write_seq($seq_obj);

最后一行write_seq()是个新东西，是不？在这一行中，序列对象$seq_obj作为write_seq()的参数传递给了SeqIO对象。从另一个角度看，SeqIO对象能够识别并处理序列对象，并将此序列对象以fasta格式输出到一个文件中。来试着运行一下这个脚本：

perl seqio.pl

>#12345 example 1
aaaatgggggggggggccccgtt

LOCUS       #12345                    23 bp    dna     linear   UNK
DEFINITION  example 1
ACCESSION   unknown
FEATURES             Location/Qualifiers
BASE COUNT        4 a      4 c     12 g      3 t
ORIGIN       1 aaaatggggg ggggggcccc gtt
//

从文件中提取序列

新手会不太习惯使用Bio::SeqIO来处理序列文件。从某些方面可以理解，因为学过了Perl一般的处理文件方式，会对Bioperl的方式不习惯，或者认为比较复杂。但使用open()函数会让提取序列变得非常复杂。请相信SeqIO！它可以处理几乎所有的序列格式、读写序列文件并且和其他很多Bioperl模块相辅相成。

前面我们已经创建了一个“sequence.fasta”的序列文件，现在我们来用SeqIO读取它，如下：

1
2
3

#!/bin/perl -w
use Bio::SeqIO;
$seqio_obj = Bio::SeqIO->new(-file => "sequence.fasta", -format => "fasta" );

与前面写文件脚本明显不同的一点是-file参数值没有了“>”符号，意思是读取“sequence.fasta”文件。现在我们加上最关键的一行，通过next_seq方法把序列文件中的序列生成序列对象。

1
2
3
4

#!/bin/perl -w
use Bio::SeqIO;
$seqio_obj = Bio::SeqIO->new(-file => "sequence.fasta", -format => "fasta" );
$seq_obj = $seqio_obj->next_seq;

next_seq()方法的作用是将下一个存在的序列提取到SeqIO对象中。在上例中，提取的是第一个序列。同时，将SeqIO对象变量命名为$seq_obj，此SeqIO对象和最开始的例子中，利用Bio::Seq自己手动创建的SeqIO对象是相似的。next_<something>方法在Bioperl中经常出现，都是类似的功能。比如，Bio::AlignIO中的next_aln()表示提取序列比对结果中的下一个比对，Bio::SearchIO中的next_hit()是提取BLAST或HMMER结果中的下一个hit。

如果文件中有多个序列，就应该用一个循环来一个一个地读取，示例：

1
2
3
4

while ($seq_obj = $seqio_obj->next_seq){
    # print the sequence   
    print $seq_obj->seq,"\n";
}

在读取序列时候，并不需要规定“-format”参数来指明文件格式。如果为了安全起见，可以声明。如果没有指定，SeqIO可以通过文件后缀名决定序列格式（Bioperl支持的序列后缀名详见SeqIO HOWTO）。在上例中，后缀“fasta”已经表明了这是一个fasta文件，SeqIO自动将其作为fasta文件读取。如果是未知的后缀名或没有后缀名，SeqIO会根据序列文件内容猜测可能的格式，但这种猜测无法保证100%正确。

指明“-alphabet”参数会比较有帮助，这样就避免Bioperl去猜测序列类型（”dna”, “rna”, “protein”），因为Bioperl有时候可能猜不对，例如MGGGGTCAATT这样的蛋白质序列会被当作DNA序列，或者类似“-”符号出现在序列中也会影响Bioperl的判断。

从网络数据库中提取一个序列

英文原文

Bioperl强大的功能之一就是可以从各种类型的资源中提取序列而不用考虑其格式，序列文件、网络数据库、本地数据库等。在此举例说明如何利用Bio::DB::Genbank模块从Genbank中提取一个序列条目。

先声明要使用的模块：

1

use Bio::DB::GenBank;

也可以从下列数据库中提取序列：SwissProt (Bio::DB::SwissProt)、GenPept (Bio::DB::GenPept)、 EMBL (Bio::DB::EMBL)、SeqHound (Bio::DB::SeqHound)、Entrez Gene (Bio::DB::EntrezGene)、 RefSeq (Bio::DB::RefSeq)等。

然后创建一个对象：

1
2

use Bio::DB::GenBank;
$db_obj = Bio::DB::GenBank->new;

这里创建了一个“数据库对象”，但并没有任何参数。再来看数据库对象的一个有用的功能：

1
2
3

use Bio::DB::GenBank;
$db_obj = Bio::DB::GenBank->new;
$seq_obj = $db_obj->get_Seq_by_id(2);

get_Seq_by_id方法 识别Genbank的GI号。另外，get_Seq_by_acc 可识别accession号 (e.g. “A12345”)，get_Seq_by_version可识别带版本号的accession号(e.g. “A12345.2”)。相应的方法只能识别相应的条目标识号。

Bioperl可以从数据库中提取一个或多个序列。但如果要提取大批量的序列时，要避免使用循环，否则NCBI会认为是滥用其服务资源而封掉使用者的IP。有其他很多更好更快的方法来大批量提取，例如，下载GenBank的某一部分数据（译者注：比如可以从NCBI ftp只下载人类的mRNA序列），然后从中直接提取需要的序列；或者使用格式化过的数据库（使用本地BLAST中formatdb格式化），然后用fastacmd（也在本地BLAST程序中）提取所需要的序列。

从网络数据库中提取多个序列

虽然有很多更复杂的方法从Genbank中提取多个序列，这里只展示一个比较“生物学”的办法，即利用Bio::DB::Query::GenBank模块。查询Genbank中拟南芥拓扑异构酶（Arabidopsis topoisomerases）的核苷酸序列，代码示例：

1
2
3

use Bio::DB::Query::GenBank;
$query = "Arabidopsis[ORGN] AND topoisomerase[TITL] and 0:3000[SLEN]";
$query_obj = Bio::DB::Query::GenBank->new(-db    => 'nucleotide', -query => $query );

注：类似这种通过字符串获取序列的办法只适用于Bioper 1.5版之后，且只能查询Genbank，至于其他数据库，如Swissprot或EMBL，只能通过标识号或accession号查询。

另一个例子，获取布氏锥虫（Trypanosoma brucei）的  EST序列

1
2
3

$query_obj = Bio::DB::Query::GenBank->new(
                      -query   =>'gbdiv est[prop] AND Trypanosoma brucei [organism]',
                      -db      => 'nucleotide' );

关于Genbank的查询字符串格式，详见此处。

刚才的例子中只是创建了一个查询对象（query object），但并没有真正提取到序列。还需要创建一个数据库对象，然后就可以按照以前讲过的在数据库对象中提取序列信息。

1
2
3
4
5
6
7
8
9
10

use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
$query = "Arabidopsis[ORGN] AND topoisomerase[TITL] and 0:3000[SLEN]";
$query_obj = Bio::DB::Query::GenBank->new(-db    => 'nucleotide',  -query => $query );
$gb_obj = Bio::DB::GenBank->new;
$stream_obj = $gb_obj->get_Stream_by_query($query_obj);
while ($seq_obj = $stream_obj->next_seq) {
    # do something with the sequence object    
    print $seq_obj->display_id, "\t", $seq_obj->length, "\n";
}

$stream_obj 和get_Stream_by_query会比较陌生。需要提取多个序列信息时，可以使用这种数据流（stream）的形式。和get_Seq_by_id比较类似，只不过可以提取多个序列。

注意上述代码是一步一步分开执行的：先是查询对象，然后传递给数据库对象，最后用一个循环来输出序列信息，循环内用的是序列对象。

警告：一定要注意查询字符串的设置。现在核苷酸数据库非常庞大，如果返回的查询结果包含了非常多的序列的话，程序会长时间运行下去，甚至耗尽内存。可以通过设置[SLEN]来限制结果数量。

创建序列（序列对象）

第一个脚本先来学习如何创建一个序列，准确的说是一个序列对象。Bioperl是“面向对象”方式的。至于为什么要用“面向对象”，Why introduce these odd or intrusive notions into software that should be biological or intuitive?（不会翻译）。原因在于模块化或者面向对象化会使得处理复杂数据的时候更灵活更简单。一旦越过了这个坎儿，你会发现使用“对象”是理所当然的。

有一种简单的方式去理解“什么是对象”，对象就是一个装载数据的箱子。举一个典型的例子，一个生物学序列包含了不同类型数据（序列本身、序列标识和注释等等），这个序列就类似一个对象。

在Bioperl中，所有的对象都是有特定的模块来创建的，如果要创建某个对象，必须先告诉Perl需要加载哪个模块。例如：

1
2

#!/bin/perl -w
use Bio::Seq;

表示用Pel来加载“Bio::Seq”模块，这个模块在电脑中是一个“<bioperl module directory>/Bio/Seq.pm”文件。Bio::Seq模块可以用于创建一个Bio::Seq对象。这个模块是Bioperl的一个核心模块。在Bioperl中，“Bio::Seq object”、“Sequence object”或者“Seq object”都是指的同一个意思，这个对象包含了一个生物学序列及其名称、标识符和属性。先来看如何创建一个简单的序列对象：

1
2
3
4

#!/bin/perl -w
use Bio::Seq;
$seq_obj = Bio::Seq->new(-seq=> "aaaatgggggggggggccccgtt",
-alphabet => ‘dna’ );

变量$seq_obj就是一个简单的，包含一个序列的序列对象。另外也需要声明这个序列是DNA序列。（-alphabet可以是‘dna’、‘rna’或者‘protein’），如果不声明，Bioperl会猜测一个合理的序列类型。一般能猜测正确，但如果序列中有太多N或X表示未知序列的字符，Bioperl很有可能猜错从而导致更严重的问题。

如前所示，可以自己创建一个对象，Bio::Seq也可以通过读取特定格式的文件来自动创建对象。比如：序列比对结果、数据库记录以及BLAST结果文件。

必须用new()方法来创建一个新的对象。“模块名称->方法名称new(一些参数名称=>参数值)”，是Bioperl的标准语法格式。

注：如果你有一些编程经验，你可能会知道“函数”或“子程序”的概念。在面向对象编程中，这个概念应该叫做“方法”。

前面说过，对象就是一个装载数据的箱子，但不止如此，这个箱子还有其他功能。一个对象会有多个方法用于调用。例如，Bio::Seq模块可以调用seq()方法来输出 Bio::Seq对象中的序列。举例说明：

1
2
3
4

#!/bin/perl -w
use Bio::Seq;
$seq_obj = Bio::Seq->new(-seq => "aaaatgggggggggggccccgtt", -alphabet => 'dna' );
print $seq_obj->seq;

这个脚本会在屏幕上输出“aaaatgggggggggggccccgtt”。->符号表示一个对象调用它的某个方法。

再来举一个更实际点的例子。一个序列对象一般都要包含标识符、描述和序列。

1
2
3
4
5
6
7

#!/bin/perl -w
use Bio::Seq;
$seq_obj = Bio::Seq->new(-seq => "aaaatgggggggggggccccgtt",
                          -display_id => "#12345",
                          -desc => "example 1",
                          -alphabet => "dna" );
print $seq_obj->seq();

其中aaaatgggggggggggccccgtt、 #12345和 example 1叫做“参数”。这个例子展示了如何将各类参数传递给new()方法。

Beginners HOWTO

原作者

briano at bioteam.net

版权

摘要

入门简介

不可避免会有一些晦涩的编程术语。如果你对编程越了解越好。不太了解也不必担心，这里只是介绍一些基本的编程知识。另外会有一些模块化和面向对象编程的基础内容，如果你需要写一些复杂的脚本，这一点是必须要掌握的。

当学习一门新技术的时候，所面临的挑战之一就是专业术语的学习。编程也不例外。耐心是最好的解决办法！你要知道，学计算机的再去学生物照样十分头疼，有木有！！！

注：bioperl中有一个模块，Bio::Perl，适用于Bioperl新手入门。（译者注：把bioperl一些基本应用转化成函数的形式调用，而不是大量的模块，更易于新手理解。）但这个HOWTO文档并不介绍这个模块，因为这只是其中一个很小的模块，不是真正的面向对象的，且不具备扩展性。这里会主要介绍Bioperl的核心，以及如何举一反三利用好Bioperl。

安装Bioperl

1. 在windows下，会出现“Error: Failed to download URL http://bioperl.org/DIST/GD.ppd”，或者“Error: Failed to download URL http://bioperl.org/DIST/GD.ppd”等错误信息。这说明有一些安装Bioperl必须的Perl模块没有被安装，你需要手动安装这些模块。详见“在windows上安装Bioperl”。
2. 在Unix下，类似“Can’t locate <some module>.pm in @INC…”的错误信息说明这个模块没有安装。详见“在Unix下安装Bioperl”。
3. “Tests Failed”错误信息表明在安装过程中有些模块安装是测试不成功，可能会影响今后的使用。Bioperl大约有1一千多个模块，安装时会进行一万多个测试。（译者注：复杂性难以想象，很难避免出错。）一般这种情况出现在GD模块上，GD只和Bio::Graphics相关，不影响其他的模块；XML parser也可能出现这种情况，但也只是影响到读取XML文件。如果要把所有模块安装成功，需要大量的工作。（译者注：新手勿入）

获取帮助

需要提供的信息包括但不限于以下几点：

1． 所使用的Bioperl版本号

2． 操作系统

3． 你的目的

4． 代码

5． 所有出现的错误信息

总是会有人给bioperl-l发邮件抱怨自己的问题没有被回复。一般原因是因为没有提供以上的所有信息，通常是没有提供代码和错误信息。

关于Perl

1. Learning Perl是最流行的perl入门书籍。（译者注：google之，很多免费下载链接，包括中文翻译版《perl语言入门》。）
2. Perl in a Nutshell也很好。提供的例子可能并不那么好，但覆盖了大量的内容。
3. Perl有自己的文档。试着在命令行中输入“perldoc perl”查看perl帮助文档。对于特定的模块，输入“perldoc <模块名字>”，可以查看关于此模块的帮助文档。举例说明：

perldoc Bio::SeqIO

写脚本

在Unix环境下，一般在命令行（shell）环境下工作。可在命令行中用这个命令查看Perl版本（译者注：windows下同样有效。以下命令除特殊说明外，默认在Unix系统下）：

perl –v

This is perl, v5.10.0 built for cygwin-thread-multi-64int

Copyright 1987-2007, Larry Wall Perl may be copied only under the terms of
either the Artistic License or the GNU General Public License, which may be
found in the Perl 5 source kit.

Complete documentation for Perl, including FAQ lists, should be found on this
system using “man perl” or “perldoc perl”.  If you have access to the Internet,
point your browser at http://www.perl.org/, the Perl Home Page. 最好使用最新版本，5.4一下的版本可能会有些问题。查看Perl程序所在的位置

which perl

/bin/perl

开始写脚本：

emacs seqio.pl

1

#!/bin/perl