Bioops

Bioinformatics=(ACGAAG->AK)+(#!/bin/sh)+(P(A|B)=P(B|A)*P(A)/P(B))

Bioperl从文件中提取序列(Bioperl HOWTO翻译4)

| Comments

从文件中提取序列

英文原文

新手会不太习惯使用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()是提取BLASTHMMER结果中的下一个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的判断。

Comments