Bioperl操作指南

Bioperl 操作指南
camelbbs@
Bioperl 为许多经典的生物信息学程序提供了软件模块，这些包括： 从本地或远程数据库获取数据； 转换数据库或文件记录的格式； 操作单个序列； 搜索相似序列；
创建和进行序列比对；
搜索基因组上的基因及其它结构； 发展机器可读的序列注释；
下面的章节将描述bioperl 怎样执行这些任务；
III.1从本地和远程数据库中获取数据
bioperl 主要集中于序列操作，但是在用bioperl 操作序列之前，需要获取序列数据。

现在你可以直接将序列数据输入到bioperl 的Seq 对象，例如：
$seq = Bio::Seq->new(-seq => 'actgtggcgtcaact',
-desc => 'Sample Bio::Seq object', -display_id => 'something', -accession_number => 'accnum', -alphabet => 'dna' );
然而，在大多数时候，从在线文档及数据库中获取序列更优越。

注意在生物信息学的传统叫法中有时候被称作“数据库”的很可能是一个“索引平台文件”。

Bioperl 支持远程数据获取，也可为访问本地数据库创建索引。

有两个普通的方法完成这个。

如果你知道序列储存在什么样的数据库中（例如文本文件、本地关系型数据库或一个internet 上可访问的远程数据库），你可以写一个脚本特定地从这些数据库中获得数据。

这种方法将在III.1.1 节和III.1.2节中描述，这两节分别讲如何从远程数据库和本地的索引平台文件中获取数据。

明确地从本地关系型数据库中获取序列数据需要安装和设置bioperl-db 库和BioSQL 计划中的模块，更多介绍可见IV .3节。

另一个方法是使用最近发展起来的OBDA(Open Bioinformatics Data Access)注册系统。

使用OBDA 可以从一个数据库中输出序列而不需要知道可访问的数据库是平台文件还是关系型，甚至不管它是本地的还是仅能从网上获得的。

关于怎样安装必须的注册配置文件和获取序列数据已在doc/howto 中的BIODATABASE_ACCESS 中叙述，这里不再重复。

III.1.1 获取远程数据（Bio::DB::GenBank 等）
Bioperl 可以直接从主要的分子生物学数据库中获取序列数据。

数据可以通过序列的获取号或id 来获得。

生物秀-专心做生物
w w w .b b i o o .c o
m Killer
Digitally signed by Killer DN: cn=Killer, o, ou,
email=Killer@, c=US Date: 2011.02.25 17:39:18 +08'00'
还可以进行批量获取以方便地获取多重序列。

例如，从genbank 中获取数据的代码如下：
$gb = new Bio::DB::GenBank(); # this returns a Seq object :
$seq1 = $gb->get_Seq_by_id('MUSIGHBA1'); # this returns a Seq object :
$seq2 = $gb->get_Seq_by_acc('AF303112'); # this returns a SeqIO object :
$seqio = $gb->get_Stream_by_id(["J00522","AF303112","2981014"]);
更多信息见III.2.1关于使用SeqIO 对象。

Bioperl 当前支持从genbank,genpept,Refseq,swissprot 和EMBL 数据库中获取序列数据。

更多信息见 Bio::DB::GenBank manpage, the Bio::DB::GenPept manpage, the Bio::DB::SwissProt manpage, the Bio::DB::RefSeq manpage 和the Bio::DB::EMBL manpage 。

用户可以为一个数据库指定不同的数据库镜像－－特别相关的是SwissProt 资源有许多ExPaSy 镜像。

还有一些为那些后台防火墙指定本地的代理服务器的配置项。

获取NCBI RefSeqs 序列可以通过一个叫Bio::DB:RefSeq 的特定模块来查询EBI 服务器。

使用之前可以参考Bio::DB::RefSeq manpage ，有一些关于获取RefSeq 的告诫。

RefSeq 在Genbank 中的id 一般是以"NT_", "NC_", "NG_", "NM_", "NP_", "XM_", "XR_", or"XP_"开头
（更多信息见/LocusLink/refseq.html ）。

严格来说Bio::DB:GenBank 可被用来获取与这些id 一致的entries ，但是记住它们不是Genbank entries 。

关于获取以"NT_"开头的entries 的特定细节见Bio::DB::GenBank manpage ，这些是经过特别格式化的"CONTIG" entries 。

Bioperl 也支持从一个远程Ace 数据库中获取序列。

这个功能需要另外的AcePerl 模块。

你需要下载并安装aceperl 模块，见/AcePerl/。

另一个模块BioFetch ，可用于访问远程数据库，它会查询EBI 中的dbfetch 脚本。

适用数据库如EMBL,GenBank 或SWALL ，可以通过不同格式如objects 或streams(SeqIO objects)，或"tempfiles"来获取entries 。

详细内容见Bio::DB::BioFetch manpage 。

III.1.2索引和访问本地数据库(Bio::Index::*, bp_index.pl, bp_fetch.pl, Bio::DB::*)
通过Bio:Index 或Bio::DB::Fasta 对象，Bioperl 允许索引本地序列文件。

下面的序列数据格式支持Bio::Index: genbank,swissprot,pfam,embl 和fasta 。

一旦使用Bio::Index 索引一组序列，就可用与访问远程数据库相似的句法来获取单条序列。

例如，如果想设置一个fasta 格式的索引平台文件数据库，并获取一个文件，可以这样写代码：
# script 1: create the index
use Bio::Index::Fasta; # using fasta file format
use strict; # some users have reported that this is necessary
my $Index_File_Name = shift; my $inx = Bio::Index::Fasta->new(
生物秀-专心做生物
w w w .b b i o o .c o m
-filename => $Index_File_Name, -write_flag => 1);
$inx->make_index(@ARGV);
# script 2: retrieve some files use Bio::Index::Fasta;
use strict; # some users have reported that this is necessary
my $Index_File_Name = shift;
my $inx = Bio::Index::Fasta->new($Index_File_Name); foreach my $id (@ARGV) {
my $seq = $inx->fetch($id); # Returns Bio::Seq object # do something with the sequence }
为了方便创建和使用更复杂或更灵活的索引系统，bioperl 在scripts/index 文件夹下分配了两个样本脚本，bp_index.PLS 和bp_fetch.PLS 。

这些脚本能够被用作模板来发展用户化的本地数据文件索引系统。

Bioperl 还提供了Bio:DB:Fasta 作为索引和查询Fasta 格式文件的手段。

它与Bio::Index::Fasta 的精髓是一致的但是提供了更多的方法，例如：
use Bio::DB::Fasta; use strict;
my $db = Bio::DB::Fasta->new($file); # one file or many files my $seqstring = $db->seq($id); # get a sequence as string my $seqobj = $db->get_Seq_by_id($id); # get a PrimarySeq obj
my $desc = $db->header($id); # get the header, or description line
关于这个模块的全部特征的更多信息见Bio::DB::Fasta manpage 。

这两个模块还可以让用户在fasta 头文件中指明一个特定字符串作为识别id ，如字符串"gi|4556644|gb|X45555"中的gi number 。

看下面的fasta 格式序列，它储存在文件"test.fa"中：
>gi|523232|emb|AAC12345|sp|D12567 titin fragment
MHRHHRTGYSAAYGPLKJHGYVHFIMCVVVSWWASDVVTYIPLLLNNSSAGWKRWWWIIFGGE GHGHHRTYSALWWPPLKJHGSKHFILCVKVSWLAKKERTYIPKKILLMMGGWWAAWWWI
默认情况下，Bio::Index::Fasta 和Bio::DB::Fasta 会使用他们在fasta 头文件中遇到的第一个词作为获取关键词，如本例中的"gi|523232|emb|AAC12345|sp|D12567"。

一个单独的id 作为一个关键词会更有用。

下面的代码将索引"test.fa"并创建一个索引文件"test.fa.idx"，其关键词是Swissprot ，或"sp"标识符。

$ENV{BIOPERL_INDEX_TYPE} = "SDBM_File"; # look for the index in the current directory $ENV{BIOPERL_INDEX} = "."; my $file_name = "test.fa";
my $inx = Bio::Index::Fasta->new( -filename => $file_name . ".idx",
生物秀-专心做生物
w w w .b b i o o .c o m
-write_flag => 1 ); # pass a reference to the critical function to the Bio::Index object $inx->id_parser(\&get_id); # make the index
$inx->make_index($file_name);
# here is where the retrieval key is specified sub get_id {
my $header = shift;
$header =~ /^>.*\bsp\|([A-Z]\d{5}\b)/; $1; }
这里是关于怎样获取序列的方法，如用一个Bio::Seq 对象：
my $seq = $inx->fetch("D12567"); print $seq->seq;
如果你想用一个Swissprot id 或一个gi number 获取一条序列，而fasta 头文件实际上是多个gi 和Swissprot 串联的头文件，该如何操作？
>gi|523232|emb|AAC12345|sp|D12567|gi|7744242|sp|V11223 titin fragment
修改那些传递到id_parser 方法的功能：
sub get_id {
my $header = shift;
my (@sps) = $header =~ /^>.*\bsp\|([A-Z]\d{5})\b/g; my (@gis) = $header =~ /gi\|(\d+)\b/g; return (@sps,@gis); }
Bio::DB::Fasta 模块使用相同的原理，但是句法稍微有所不同，例如：
my $db = Bio::DB::Fasta->new('test.fa', -makeid=>\&make_my_id); my $seqobj = $db->get_Seq_by_id($id); sub make_my_id {
my $description_line = shift;
$description_line =~ /gi\|(\d+)\|emb\|(\w+)/; ($1,$2); }
bioperl 的核心安装文件并不支持访问储存在关系数据库中的序列和数据。

这种功能由附属的bioperl-db 文库提供。

更多信息见IV .3节。

III.2转换数据库或文件记录格式
生物秀-专心做生物
w w w .b b i o o .c o m
III.2.1转换序列文件(SeqIO)
生物信息学的一个普通的任务就是在广泛使用的序列数据格式之间进行转换。

Bioperl 的SeqIO 对象使这件事变得轻而易举。

SeqIO 能够阅读大量的序列格式（在一个或多个文件中）：Fasta,EMBL,GenBank,Swissprot,PIR,GCG,SCF, phd/phred, Ace, fastq, exp, chado, or raw (plain sequence)。

SeqIO 也能够分析alf, ztr, abi, ctf, 和ctr 格式的文件，一旦序列数据被SeqIO 读取，就变成bioperl 的Seq 、PrimarySeq 或RichSeq 对象所能利用的形式，选用这几个对象要看序列的来源。

此外，序列对象可以任何支持的数据格式被写入另一个文件(也用SeqIO)，使数据转换简单并可执行，例如：
use Bio::SeqIO;
$in = Bio::SeqIO->new(-file => "inputfilename", -format => 'Fasta');
$out = Bio::SeqIO->new(-file => ">outputfilename", -format => 'EMBL');
while ( my $seq = $in->next_seq() ) {$out->write_seq($seq); }
另外，perl 的"tied 句柄"句法是SeqIO 可用的，你可以用standard<> 和打印操作来读写序列对象，例如：
$in = Bio::SeqIO->newFh(-file => "inputfilename" , -format => 'fasta'); $out = Bio::SeqIO->newFh(-format => 'embl'); print $out $_ while <$in>;
如果不使用"-format"语句，Bioperl 会依赖于文件后缀来决定格式，这是一个较慢的办法。

如果没有后缀可用，SeqIO 会基于实际内容推测序列的格式。

这里是当前通用的一组后缀：
Format Suffixes Comment fasta fasta|fast|seq|fa|fsa|nt|aa Fasta
genbank gb|gbank|genbank|gbs|gbk Genbank
scf scf SCF tracefile pir pir PIR embl embl|ebl|emb|dat EMBL raw txt plain gcg gcg GCG ace ace ACeDB bsml bsm|bsml BSML XML game GAME XML swiss swiss|sp SwissProt phd phd|phred Phred fastq fastq Fastq
Locuslink LL_tmpl format qual Phred quality file chado Chado XML
tinyseq NCBI TinySeq XML exp exp Staden experiment file abi* abi ABI tracefile
生物秀-专心做生物
w w w .b b i o o .c o m
alf* alf ALF tracefile ctf* ctf CTF tracefile ztr* ztr ZTR tracefile
pln* pln Staden plain tracefile
*这些格式需要bioperl-ext 包和来自于Staden 包的io_lib 文库
更多信息见Bio::SeqIO manpage 或SeqIO HOWTO (/HOWTOs/html/SeqIO.html)
III.2.2转换比对文件(AlignIO)
数据文件储存的多序列比对也具有不同的格式。

AlignIO 是bioperl 中转换比对文件格式的对象。

AlignIO 基本用法类似于SeqIO 对象，它的许多命令的名字与SeqIO 相同。

如在SeqIO 和AlignIO 对象中都可以用"-file"和"-format"来创建文件：
use Bio::AlignIO;
my $io = Bio::AlignIO->new(-file => "receptors.aln", -format => "clustalw" );
如果"-format"语句不起作用， Bioperl 就通过文件后缀来决定格式，下面是目前常用的一组后缀：
格式 后缀 内容 bl2seq
clustalw aln
emboss* water|needle
fasta fasta|fast|seq|fa|fsa|nt|aa maf maf
mase Seaview mega meg|mega meme meme metafasta
msf msf|pileup|gcg GCG nexus nexus|nex pfam pfam|pfm
phylip phylip|phlp|phyl|phy|phy|ph interleaved prodom
psi psi PSI-BLAST selex selex|slx|selx|slex|sx HMMER stockholm
*water, needle, matcher, stretcher, merger, 和supermatcher 见IV .2.1的EMBOSS 。

与SeqIO 不同的是，AlignIO 不能创建每种格式的输出文件。

AlignIO 当前支持6种输出格式： fasta, mase, selex, clustalw, msf/gcg 和phylip (交叉存取)。

AlignIO 与SeqIO 的另一个重要区别是AlignIO 一次只能操作一个比对IO ，而SeqIO.pm 可以在单串中操作
生
物秀-专心做生物
w w w .b b i o o .c o m
多个序列。

AlignIO 的句法基本上与SeqIO 一致：
use Bio::AlignIO;
$in = Bio::AlignIO->new(-file => "inputfilename" , -format => 'fasta');
$out = Bio::AlignIO->new(-file => ">outputfilename", -format => 'pfam');
while ( my $aln = $in->next_aln() ) { $out->write_aln($aln); }
仅有的区别是返回的对象的条目，$aln ，是针对SimpleAlign 对象而不是Seq 对象。

AlignIO 也支持上面描述的SeqIO 的tied 句柄的句法。

更多关于SimpleAlign 的信息见Bio::AlignIO manpage, Bio::SimpleAlign manpage 和III.5。

III.3 操作序列
Bioperl 包括许多序列分析模块。

如果在bioperl 中没有找到你要找的功能，就在EMBOSS 或PISE 中找，它们可在bioperl-run 附加文库中得到（见IV .2.1）。

III.3.1使用Seq 方法操作序列
现在我们知道了如何获取序列并将它们作为序列对象访问。

再来看看怎样使用序列对象去操作我们的序列数据并返回信息。

Seq 提供了多种方法来执行许多普通（有时候不普通）的操作序列和返回数据的任务。

这里列了一些最有用的：
这些方法返回字符串或者被用来设定值：
$seqobj->display_id(); # the human read-able id of the sequence $seqobj->seq(); # string of sequence
$seqobj->subseq(5,10); # part of the sequence as a string $seqobj->accession_number(); # when there, the accession number $seqobj->alphabet(); # one of 'dna','rna','protein'
$seqobj->primary_id(); # a unique id for this sequence irregardless # of its display_id or accession number $seqobj->desc(); # a description of the sequence
值得一提的是有些值对应于给定格式的特定区域。

例如，display_id 方法返回一个Genbank entry 的LOCUS 名字，在一个Fasta 文件中，(\S+)跟随>，ID 来自于一个SwissProt 文件，等等。

desc()方法将返回一个 Genbank 文件中的DEFINITION 行，在Fasta 文件中，这一行跟着display_id ，SwissProt 文件中在DE 区。

下面的方法返回一组Bio::SeqFeature 对象：
$seqobj->get_SeqFeatures; # The 'top level' sequence features $seqobj->get_all_SeqFeatures; # All sequence features, including sub- # seq features
对于一条comment 注释，你可以使用：
生物秀-专心做生物
w w w .b b i o o .c o m
use Bio::Annotation::Comment;
$seq->annotation->add_Annotation('comment',
Bio::Annotation::Comment->new(-text => 'some description');
对于一条reference 注释，可以使用：
use Bio::Annotation::Reference;
$seq->annotation->add_Annotation('reference',
Bio::Annotation::Reference->new(-authors => 'author1,author2', -title => 'title line', -location => 'location line', -medline => 998122 );
序列特征将在III.7的机器可读序列注释中讨论。

对这个对象的总体描述在Bio::SeqFeature::Generic manpage 中可以找到，一个相关的高水平的注释描述在Bio::Annotation::Collection manpage 中。

另外的关于获取序列特征的样本代码可在gb2features.pl 中找到，它在子文件夹examples/DB 下。

最后，还有关于特征和注释的HOWTO (/HOWTOs/html/Feature-Annotation.html) 和一节在FAQ 中关于特征的介绍(/Core/Latest/faq.html#5)。

下面的方法返回新的序列对象，但是并不把起始对象的特征传递到结果特征中：
$seqobj->trunc(5,10); # truncation from 5 to 10 as new object $seqobj->revcom; # reverse complements sequence $seqobj->translate; # translation of the sequence
注意到一些方法返回字符串，一些返回队列，一些返回对象。

更多信息见Bio::Seq manpage 。

许多方法自己都带有注释。

然而，bioperl 灵活的方法准许有更多的注释。

生物信息学中的翻译意味着两种稍有不同的事件：
1.从头到尾翻译一段核苷酸序列。

2.重视mRNAs 中真实编码区的局限性。

bioperl 关于执行序列的翻译可以很简单的完成这些任务的第一步。

任何不是蛋白质字母表的序列对象，都可以通过返回一个蛋白序列对象的方法被翻译：
$translation1 = $my_seq_object->translate;
然而，翻译方法也可以通过传递许多可选参数来改变作用。

例如，translate()的前两句可以用来修饰用于代表终止(default'*')和未知氨基酸('X')的参数。

(正常情况下最好是左边未被改变)。

第三句决定翻译框。

默认的框是"0"。

用其他两种前移框来翻译,可以这样写：
$translation2 = $my_seq_object->translate(undef,undef,1); $translation3 = $my_seq_object->translate(undef,undef,2);
生物秀-专心做生物
w w w .b b i o o .c o m
translate()的第四个语句使它可以选择使用遗传密码。

目前有16个编码平台，包括
'V ertebrate Mitochondrial', 'Bacterial', 'Alternative Yeast Nuclear' 和'Ciliate, Dasycladacean and Hexamita Nuclear' 翻译平台。

这些平台在Bio::Tools::CodonTable 对象中，被翻译方法所使用。

例如，线粒体翻译：
$human_mitochondrial_translation = $seq_obj->translate(undef,undef,undef,2);
如果我们想翻译全部编码区(CDS)，主要是核酸数据库EMBL 、GenBank 和DDBJ 这样做，翻译方法不得不执行更多的技巧。

特别是，'translate'需要确定序列的开头和末尾有适当的起始和终止密码子，而序列内部没有终止密码子。

另外，如果使用的遗传密码中有一个非典型(非A TG)起始密码子，翻译方法需要将起始氨基酸转换为亮氨酸。

这些检查和转换可以通过设定翻译方法的第五句求真值而触发。

如果第五句设为真而没有一个标准来找合适的CDS ，默认情况下就会出现一种警告。

通过设定第六句来求真值，如果发现不合适的CDS ，就通知程序进行消除。

例如：
$protein_object = $cds->translate(undef,undef,undef,undef,1,'die_if_errors');
相关细节见Bio::Tools::CodonTable manpage 。

III.3.2获得序列的基本统计数据(SeqStats,SeqWord)
为了增加Seq 对象中的直接可用的方法，bioperl 提供了多个帮助对象来测定一条序列更多的信息。

例如，SeqStats 对象提供了来获取序列的分子质量的方法，以及每种残基（核酸的碱基或蛋白质的氨基酸）的出现数量。

对核酸来说，SeqStats 也可以返回使用的密码子的数量统计。

例如：
use SeqStats;
$seq_stats = Bio::Tools::SeqStats->new($seqobj); $weight = $seq_stats->get_mol_wt();
$monomer_ref = $seq_stats->count_monomers();
$codon_ref = $seq_stats->count_codons(); # for nucleic acid sequence
注意：有时候序列会包含模糊密码。

因为这个原因，get_mol_wt()返回一个条目给一个包含分子质量的最大下限和最小上限的两个元件的队列。

SeqWords 对象与SeqStats 相似并提供一种方法计算序列中"单词"（例如tetramers 或hexamers ）的出现频率。

更多信息请见Bio::Tools::SeqStats manpage 和Bio::Tools::SeqWords manpage 。

III.3.3识别限制性酶切位点(Bio::Restriction)
操作核酸序列的另一个常见任务是定位限制性酶切位点。

Bioperl 为这个目的提供了Bio::Restriction::Enzyme ，Bio::Restriction::EnzymeCollection 和Bio::Restriction::Analysis 对象。

这些模块取代旧的模块Bio::Tools::RestrictionEnzyme. 一种新的收集酶的对象定义如下：
use Bio::Restriction::EnzymeCollection;
my $all_collection = Bio::Restriction::EnzymeCollection;
生物秀-专心做生物
w w w .b b i o o .c o m
Bioperl 默认的Restrcition::EnzymeCollection 对象返回500多个不同的II 型限制性酶的数据。

用available_list()方法可以得到可用酶的列表，但这些仅是名称，不是功能对象。

你也可以访问酶的子集。

例如选择所有的具有6碱基长识别位点的可用酶对象，代码如下：
my $six_cutter_collection = $all_collection->cutters(6); foreach my $enz ($six_cutter_collection){
print $enz->name,"\t",$enz->site,"\t",$enz->overhang_seq,"\n"; # prints name, recognition site, overhang }
还有一些方法也可用来选择酶对象组，如unique_cutters()和blunt_enzymes()。

你也可以通过名称来选择一个酶对象，如：
my $ecori_enzyme = $all_collection->get_enzyme('EcoRI');
一旦选择了一个合适的酶对象，就可用fragments()方法获得给定核酸序列上的酶位点。

执行这个任务的句法：
use Bio::Restriction::Analysis;
my $analysis = Bio::Restriction::Analysis->new(-seq => $seq); # where $seq is the Bio::Seq object for the DNA to be cut @fragments = $analysis->fragments($enzyme); # and @fragments will be an array of strings
为了得到同裂酶，甲基化位点和微生物资源的信息，你需要从一个REBASE 文件里直接创建你的EnzymeCollection ，如：
use Bio::Restriction::IO;
my $re_io = Bio::Restriction::IO->new(-file=>$file,-format=>'withrefm'); my $rebase_collection = $re_io->read;
正确格式的REBASE 文件可以在ftp:///pub/rebase 找到，一般它有一个名字如"withrefm.308"。

如果需要的话你也可以创建新的酶，如：
my $re = new Bio::Restriction::Enzyme(-enzyme=>'BioRI',-seq=>'GG^AA TTCC');
更多信息请看 Bio::Restriction::Enzyme manpage, Bio::Restriction::EnzymeCollection manpage, Bio::Restriction::Analysis manpage, 和 Bio::Restriction::IO manpage.
III.3.4识别氨基酸裂解位点（Sigcleave ）
对于氨基酸序列我们有兴趣知道氨基酸序列是否包含一个可分裂的信号序列来指导细胞内蛋白质的运输。

SigCleave 是一个预测信号序列的程序（最初是EGCG 分子生物学包的一部分），同时基于von Heijne 算法识别裂解位点。

初始设置可以控制分数报告。

如果用户没有传递初始值，密码子默认为一个报告值3.5。

SigCleave 将只返
生物秀-专心做生物
w w w .b b i o o .c o m
回满足初始限制的分数/位置信息。

使用Sigcleave 句法如下：
# create a Seq object, for example:
$seqobj = Bio::Seq->new(-seq => "AALLHHHHHHGGGGPPRTTTTTVVVVVVVVVVVVVVV"); use Bio::Tools::Sigcleave;
$sigcleave_object = new Bio::Tools::Sigcleave ( -seq => $seqobj, -threshold => 3.5,
-desc => 'test sigcleave protein seq', -type => 'AMINO' );
%raw_results = $sigcleave_object->signals; $formatted_output = $sigcleave_object->pretty_print;
注意SigCleave 对象中"Type"是指"amino"，而在Seq 对象中它是"protein"。

更多细节请看Bio::Tools::Sigcleave manpage 。

III.3.5显示各种复杂的序列信息：OddCodes,SeqPattern
OddCodes:
有时候有必要显示一条氨基酸序列的各种信息，如疏水性氨基酸的位置或者是带正电的地方。

Bioperl 可以用Bio::Tools::OddCodes 实现这个目的。

例如，快速查看带电氨基酸在序列中的位置可以这样实现：
use Bio::Tools::OddCodes;
$oddcode_obj = Bio::Tools::OddCodes->new($amino_obj); $output = $oddcode_obj->charge();
每条序列都可被转换为成一种三字母序列（A ，C ，N ）即负电（酸性），正电（碱性）和中性氨基酸。

例如ACDEFGH 会成为NNAANNC 。

为了得到更完全的序列的化学特征，可以用chemical()方法将序列转化成一种8字母的化学字母表 { A (acidic), L (aliphatic), M (amide), R (aromatic), C (basic), H (hydroxyl), I (imino), S (sulfur) }:
$output = $oddcode_obj->chemical();
在这种情况下例子序列ACDEFGH 就成为LSAARAC 。

OddCodes 也提供到字母表的翻译以显示氨基酸序列的一些特征如疏水性、“功能性”或使用Dayhoff 定义来分组。

详细内容见文档Bio::Tools::OddCodes manpage 。

生物秀-专心做生物
w w w .b b i o o .c o m
SeqPattern 对象使用perl 正则表达式来操作序列。

设计SeqPattern 的目的是为了产生一条包括模糊碱基和（或）正则表达式的核酸序列模式的反转互补序列。

当提交的一条序列的正义和反义链都需要模式匹配的时候，就可以用这种方法来实现。

典型的SeqPattern 句法如下。

一些有趣的例子在examples/tools 文件夹的Seq_Pattern.pl 代码中可以提供更多的信息。

use Bio::Tools::SeqPattern;
$pattern = '(CCCCT)N{1,200}(agggg)N{1,200}(agggg)'; $pattern_obj = new Bio::Tools::SeqPattern(-SEQ => $pattern,
-TYPE => 'dna'); $pattern_obj2 = $pattern_obj->revcom();
$pattern_obj->revcom(1); # returns expanded rev complement pattern.
更多的细节可见Bio::Tools::SeqPattern manpage 。

III.3.6 同等系统转换（Coordinate::Pair, RelSegment ）
同等系统转换是一个普通的需求，例如，当你想查看序列特征的相关位置并转换这些相关位置到分离同等通过一条染色体或一个contig 。

尽管同等的转换听起来很漂亮，当包括转换到同等的可能性在负链的时候或者有一个同等的系统末端因为你到达一个克隆或contig 的末端。

Bioperl 有两个不同的方法来进行同等系统的转换（基于模块Bio::Coordinate::Pair 和 Bio::DB::GFF::RelSegment ）。

Coordinate::Pair 方法是有时是更“低水平”。

用它，你可以定义一个输入的同等的系统和一个输出的同等系统，每种情况下一个同等系统就是一个起始位置、终止位置和链的三重体。

终止位置是特别重要的当处理一个未结束的类（当达到一个克隆或contig 末端时同等系统所终止的地方）。

一旦两种终止系统被定义，就需要定义一个Coordinate::Pair 来在他们之间绘图。

然后一个人可以在终止系统之间进行图谱定位，代码如下：
$input_coordinates = Bio::Location::Simple->new
(-seq_id => 'propeptide', -start => 1000, -end => 2000, -strand=>1 ); $output_coordinates = Bio::Location::Simple->new
(-seq_id => 'peptide', -start => 1100, -end => 2100, -strand=>1 ); $pair = Bio::Coordinate::Pair->new
(-in => $input_coordinates , -out => $output_coordinates ); $pos = Bio::Location::Simple->new (-start => 500, -end => 500 ); $res = $pair->map($pos);
$converted_start = $res->start;
在这个例子中$res 也是一个Bio::Location 对象，如你期望的那样。

更多细节见Bio::Coordinate::Pair 和 Bio::Coordinate::GeneMapper 文档。

Bio::DB::GFF::RelSegment 方法被设计为操作序列特征的同等转换而不是武断的转换同等系统。

使用Bio::DB::GFF::RelSegment 你可以定义一个相关的同等系统给一个特别的特征（称作“refseq ”）。

你也可以访问单独的同等系统（典型的如整条染色体）。

通过重新定义相应的参考特征（如“refseq ”）你能简单决定一个与其他特征相关的特征的位置。

代码如下： 生物秀-专心做生物
w w w .b b i o o .c o m
$db = Bio::DB::GFF->new(-dsn => 'dbi:mysql:elegans', -adaptor => 'dbi:mysqlopt'); $segment = $db->segment('ZK909');
$relative_start = $segment->start; # $relative_start = 1;
# Now retrieve the start position of ZK909 relative to feature ZK337 $segment->refseq('ZK337');
$relative_start = $segment->start;
# Now retrieve the start position of ZK909 relative to the entire chromosome $absolute_start = $segment->abs_start;
这个方法很方便因为你不必直接跟踪相同的踪迹，你只需跟踪一个特征的名字（在相配系统起始时即被标记）。

还有，这个方法不需要你已经储存所有的序列特征为GFF 格式。

更多的，Bio::DB::GFF::RelSegment 已经主要发展并在一个储存了所有序列特征的Bioperl-db 关系型数据库中测试了其应用。

然而，如果一个人想使用Bio:DB::GFF 机制（包括它的同等转换能力）而不建立一个本地关系型数据库，就可以通过定义'databse'有一个别号叫'memory'，例如：
$db = Bio::DB::GFF->new( '-adaptor' => 'memory' );
更多关于bioperl 中同等转换和其它GFF 相关能力的细节见Bio::DB::GFF::RelSegment manpage ，Bio::DB::GFF manpage ，和测试文件t/BioDBGFF.t 。

III.4 搜索相似序列
分子生物学的一个基本任务就是识别那些与所感兴趣的序列相似的序列。

由NCBI 最先发展的Blast 程序，被广泛用于识别这样的序列。

bioperl 和bioperl-run 包提供了许多模块来使Blast 运行更加方便，同时可以分析Blast 产生的长篇报告。

III.4.1 运行BLAST （使用RemoteBlast.pm ）
Bioperl 通过RemoteBlast 对象来执行NCBI 的远程blast 。

运行远程blast 的代码骨架如下：
$remote_blast = Bio::Tools::Run::RemoteBlast->new (
-prog => 'blastp',-data => 'ecoli',-expect => '1e-10' ); $r = $remote_blast->submit_blast("t/data/ecolitst.fa"); while (@rids = $remote_blast->each_rid ) {
foreach $rid ( @rids ) {$rc = $remote_blast->retrieve_blast($rid);} }
你可能想改变这种远程工作的一些参数，这个例子显示了如何改变其来源：
$Bio::Tools::Run::RemoteBlast::HEADER{'MA TRIX_NAME'} = 'BLOSUM25';
关于各种CGI 参数的描述可以见：/BLAST/Doc/urlapi.html 生物秀-专心做生物
w w w .b b i o o .c o m
注意这个脚本分两个部分。

Blast 提交和随后的结果反馈。

当NCBI Blast 负荷过大时，提交和获取结果之间有一段时间间隔。

$rc 对象包括blast 的报告，可以用Bio::Tools::BPlite 或 Bio::SearchIO 来处理。

1.0版之后的默认的返回对象是SearchIO 。

对象的类型可以通过-readmethod 参数来改变，但是最好的Blast 解析子是Bio::SearchIO ，其它的在后来的版本中都不支持。

注意，要使这个代码实际可用，应该加上一些细节如检查从Blast 返回的代码来看它是否成功并且是否有一个"sleep"循环等待对NCBI 服务器的连续的请求。

见附录中的代码示例22可以找到你可以使用的一些工作代码，更多的细节见Bio::Tools::Run::RemoteBlast manpage 。

还要注意创建一个远程blast 工作的句法与创建StandAloneBlast,Clustalw ，和T-Coffee 工作的句法稍微有所不同。

RemoteBlast 特别需要一个导向的连字符来传递参数，如 '-prog' => 'blastp'，而别的程序不使用连字符来传递参数。

III.4.2用Search 和SearchIO 解析BLAST 和FASTA 报告
不管Blast 如何运行（本地或者远程，使用或不使用perl 界面），它们都返回大量待筛选的数据。

Bioperl 提供了许多不同的对象－Search.pm/SearchIO.pm 和BPlite.pm （连同其辅助修饰，BPpsilite 和BPbl2seq ）来解析Blast 报告。

这一节将要描述bioperl 主要的Blast 和FASTA 报告解析器Search 和SearchIO ，III.4.3描述旧版的BPlite 。

我们推荐使用SearchIO ，它在以后的所有包中都可用。

Search 和SearchIO 模块提供了一个统一的界面，用来解析BLAST （以标准和BLAST XML 格式），PSI-BLAST, RPS-BLAST, bl2seq 和 FASTA 产生的序列相似性搜索报告。

SearchIO 模块也为HMMER 报告提供解析，可以设想Search/SearchIO 将来会发展成为一个统一的界面来进行更大范围的解析包括解析Genscan 。

用Search 和SearchIO 可直接解析序列相似性报告。

一开始，一个SearchIO 对象指定一个文件包含报告。

next_result 方法将下一条报告读到Search 对象，就像SeqIO 的next_seq 方法将文件中的下一条序列读到一个Seq 对象。

一旦某个报告（如一个SearchIO 对象）被读入并可用于脚本，报告的总体特征（如query ）可被决定并且它的每个hit 可用next_hit 方法访问。

每个hit 的每个高分片段都可用next_hsp 方法来访问。

除了读取一个文件中的多重报告需要额外的句法外，其它的Search/SearchIO 解析句法与被代替的BPlite 对象中的用法很相似。

读取一个BLAST 报告的代码如下：
# Get the report
$searchio = new Bio::SearchIO (-format => 'blast',
-file => $blast_report); $result = $searchio->next_result; # Get info about the entire report $result->database_name;
$algorithm_type = $result->algorithm; # get info about the first hit $hit = $result->next_hit; 生物秀-专心做生物
w w w .b b i o o .c o m
# get info about the first hsp of the first hit $hsp = $hit->next_hsp;
$hsp_start = $hsp->query->start;
关于怎样使用SearchIO 的更多的细节可见/HOWTOs/html/SearchIO.html 或docs/howto 的子目录。

更多文档可见Bio::SearchIO::blast manpage, Bio::SearchIO::psiblast manpage, Bio::SearchIO::blastxml manpage, Bio::SearchIO::fasta manpage, 和Bio::SearchIO manpage 。

在examples/searchio 文件夹中也有一个样本代码阐明如何使用SearchIO 。

最后在FAQ(/Core/Latest/faq.html#3)中也有一节关于SearchIO 的问答。

III.4.3用BPlite ，BPpsilite 和BPbl2seq 解析BLAST 报告
Bioperl 的旧版BLAST 报告解析器－BPlite, BPpsilite, BPbl2seq 和Blast.pm 不再被支持，但是作为Bioperl 代码初创的产物，它们还会在Bioperl 中保留一段时间。

BPlite 的大多数用户界面非常类似于Search ，如访问下一条hit 或HSP 用的next_Sbjct 和next_HSP 方法，分别对应于Search 的next_hit 和next_hsp 。

BPlite
BPlite 的使用句法如下，nextSbjct()用来取得hits ，获取高分片段的方法叫做nextHSP()：
use Bio::Tools::BPlite;
$report = new Bio::Tools::BPlite(-fh=>\*STDIN); $report->query;
while(my $sbjct = $report->nextSbjct) { $sbjct->name;
while (my $hsp = $sbjct->nextHSP) { print $hsp->score,"\n"; } }
对这个模块的完整描述可见Bio::Tools::BPlite manpage 。

BPpsilite
BPpsilite 和BPbl2seq 是分别用来解析（多重迭代）PSIBLAST 报告和Blast bl2seq 报告的方法。

它们都在BPlite 对象的基础上有较小的变化。

见Bio::Tools::BPbl2seq manpage 和Bio::Tools::BPpsilite manpage 。

解析多重迭代PSIBLAST 报告的句法如下。

唯一重要的是在BPlite 中增加了决定blast 迭代次数和访问每一次迭代的结果的方法。

每次迭代的结果都被作为一个（完整的）BPlite 对象以相同的方式来解析。

use Bio::Tools::BPpsilite;
$report = new Bio::Tools::BPpsilite(-fh=>\*STDIN); $total_iterations = $report->number_of_iterations; $last_iteration = $report->round($total_iterations) while(my $sbjct = $last_iteration ->nextSbjct) { 生物秀-专心做生物
w w w .b b i o o .c o m。