基因组学研究的数据分析之六：从fasta文件中提取信息

    在生物学研究中会经常遇见fasta格式文件，下面是一个fasta格式文件的例子，一条序列信息分两行，第一行以>开头包括序列相关的注释信息，第二行就是DNA或蛋白质序列，一个文件可能有一个或成千上万条这样的序列信息：

>gi|224514618|ref|NT_077402.2| Homo sapiens chromosome 1 genomic contig, GRCh37.p9 Primary Assembly
TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAAACCCTAAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCAACCCCAACCCCAACCCCAACCCCAACCCCAACCCTAACCCCTAACCCTAACCCTAACCCTACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTCGCGGTACCCTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTACCACCGAAATCTGTGCAGAGGACAACGCAGCTCCGCCCTCGCGGTGCTCTCCGGGTCTGTGCTGAGGAGAACGCAACTCCGCCGTTGCAAAGGCGCGCCGCGCCGGCGCAGGCGCAGAGAGGCGCGCCGCGCCGGCGCAGGCGCAGAGAGGCGCGCCTGGCGCAGGCGCAGAGACGCAAGCCTACGGGCGGGGGTTGGGGGGGCGTGTGTTGCAGGA

    那么我们该如何提取我们需要的信息呢？下面一段代码应用Bioperl从fasta文件提取序列信息：

#!/usr/bin/perl

use strict;

use warnings;

use Bio::SeqIO;

my $in_obj = Bio::SeqIO->new(-file => "/data/test.fasta",  -format => "fasta");

while ( my $seq_obj=$in_obj->next_seq() ) {

  my $display_id=$seq_obj->display_id();

  my $description=$seq_obj->desc();

  my $sequence=$seq_obj->seq();

  print "$display_idt", "$descriptionn";

}

    一个进行Perl编程的好习惯是进行数据分析时尽量使用已有的Perl模块，因为许多时候使用他人提供的模块往往比自己设计模块更划算，特别是许多模块已经“久经考验”。使用Bioperl就是一个例子，在Bioperl可以找到许多有用的模块，学习Bioperl最好的办法是学习已有的代码和软件帮助文档，可以在www.bioperl.org找到学习Bioperl的教程，在search.cpan.org找到相关的帮助文档。

   上面这段代码如果从生物学角度分析，就是这样：打开文件test.fasta，逐次提取序列信息，每一条序列信息包括display_id（就是“gi|224514618|ref|NT_077402.2”，包括GI号和登录号）, description（就是“Homo sapiens chromosome 1 genomic contig, GRCh37.p9 Primary Assembly”）和DNA序列，每提取一条序列信息同时输出该序列信息的display_id 和description。这段代码比较简单，最好对这些代码的布局烂熟于心，因为以后应用Bioperl的模块大部分都是类似写法，到时依样画葫芦就行，不过呢如果能够再深入了解一点编程思想，则对以后应用Bioperl会有很大帮助。

   如果从编程角度分析，上面这段代码展示典型的面向对象模块的Perl编程布局。关于面向对象编程(object-oriented programming, OOP)的学习，可以方便找到相关书籍，起码应该了解类class、对象object、属性attribute、方法method、封装envelopment、继承inheritance、多态性polymorphism、构造函数constructor和析构函数destructor等概念，这里面类class的概念相当与模块module或包package， 方法method的概念相当于函数function或子例程subroutine, 将执行一项相对完整的任务所需要的对象、数据和方法封装到一个数据结构即称为类，数据用来描述对象（即对象的属性），方法用来创建、访问、操纵或销毁对象的，类之间或者方法之间因为在一个类内或方法内对其它类或方法的调用而产生继承，因为这种继承使得扩展类或方法，由此产生多态性。

   如果从编程角度分析，上面这段代码作如下解读：代码use Bio::SeqIO表示调用类Bio::SeqIO，就是名为SeqIO.pm的文件，类SeqIO包括多个方法，如new、next_seq和write_seq等。一般的类中的第一个方法叫new（当然也可以起其它的名字），任何时候第一调用类当然包括类中的方法，必须首先调用new，作为构造函数方法负责定义（创建或者构造）指定对象，明确对象所属的类同时返回指向该对象的引用（指针），就是$in_obj = Bio::SeqIO->new(-file => "/data/test.fasta",  -format => "fasta")，注意符号"->"，左边使用类名，->表示把类名作为第一个参数传递给方法new，-file => "/data/test.fasta"和-format => "fasta"分别是方法new的第二和三个参数，同时返回指向一个匿名散列的引用$in_obj，代码$seq_obj=$in_obj->next_seq()、$display_id=$seq_obj->display_id()、$description=$seq_obj->desc()和$sequence=$seq_obj->seq()等都是在调用不同的方法，注意方法display_id、desc、seq和next_seq之间是有继承性的，它们是属于类Bio::Seq的，只要明确父类的名字Bio::SeqIO实际运行时会自动调用子类，但是在写代码时必须逐层调用，否则编译时会返回错误信息。