这一篇哪,其实主要是复习性质的。你可以松口气,边欣赏、边回味用世界上最美的程序语言(Perl)是如何解析复杂的、神秘的生物大分子序列的。
BioPerl比普通的Perl更简单吗?
如果你只是解析简单的fasta序列的话,BioPerl未必很简单,甚至刚开始还比普通的Perl更加难懂呢~但是如果你准备解析一些很复杂的序列格式(或程序输出结果),比如Genbank,SwissProt,EMBL,blast……这时BioPerl威力就显现了。
什么时候使用BioPerl?
看一下两个例子,你就明白啦~
例子一:小明、小王、小红的各科成绩。放在文件score.txt里
name math English computer
Xiaoming 90 84 91
Xiaowang 83 88 82
Xiaohong 89 94 85
例子二:三条fasta序列。放在文件sequence.fasta里
>xiaoming
CGCTATTCAACAAGGCATTCCCCCAAGTTTGAATTCTTTGAAATAGATTGCTATTAGCTA
AGATAAGAACGAAAAGGAAGGATATGGCTAAAGAAACATCTGCAATGCTTTTATTAAAAA
>xiaowang
ACCTAGCGACTCTCTCCACCGTTTGACGAGGCCATTTACAAAAACATAACGAACGACAAT
>xiaohong
GCAGACTCGGTGGAAGAGGTGATATATGGTGCAGTGTAGCTAGCAGACTCGGTGGAAGAC
AGATAAGAACGAAAAGGAAGGATATGGCTAAAGAAACATCTGCAATGCTTTTATTAAAAA
看出区别来了吗?第一个文件score.txt,每行的信息都是独立的(第1行:标题;第2行:小明成绩;第3行:小王成绩;……),相信不用我说大家也该知道了:这样的文件很适合用Perl来处理(如果你还想懒一些,甚至连Perl程序都不用写,直接使用awk之类的Unix命令行工具都行)。第二个文件sequence.fasta,不同行的信息却是有关联的(小明的在1~3行,小王的在4~5行,……)由于Perl只能一行一行地读入信息,就比较麻烦。
Genbank序列与fasta序列有什么不同?
从前面两篇中可以看出,fasta序列“非常简单”,只有序列的名称、描述以及序列字符串。当然你也可以在名称或描述中写入其它对你来说有用的信息(比如,基因的边界在哪儿?内含子的边界在哪儿?),但是BioPerl只能摘取整个描述,还是不认得你在里面写了什么。
相比而言,Genbank序列能描述的内容就非常丰富(类似的还有SwissProt,EMBL等格式的序列),除了名称、描述和序列字符串以外还有序列号、形状(线状还是环状?)、发布日期、所属物种以及序列内包含的基因、RNA、CDS、内含子(如果有的话)……What’s more,这些信息并不是你想怎么写就怎么写的,而都是有一定的规范格式,故BioPerl可以准确地识别并提取这些信息。所以,我强烈建议你不要自行修改下载得到的Genbank序列,否则BioPerl会迷路。
如何从本地文件中提取Genbank序列?
差不多该问这个问题了。如果你已经熟记了前面两篇的内容,应该会依葫芦画瓢了(先到NCBI网站上下载一条Genbnak序列来联系吧!但要短一点的那种)。假设我们下载了一条Genbank序列保存在一个叫做ecorho.gbk 的文本文件里(扩展名不是必需的,我只是为了便于辩认),就能很快地写出如下代码:
use Bio::SeqIO; $catch_seq = Bio::SeqIO -> new(-file => 'ecorho.gbk', -format => 'genbank'); $seq_obj = $catch_seq -> next_seq;
我在这里偷懒了一点:第二句代码中创建的SeqIO对象我直接写成\$catch_seq,把后缀名seqio_obj省掉了。如果你一直不放心\$catch_seq是什么的话,还是把后缀加上去比较好。也可以使用ref函数来检查一下$catch_seq究竟是什么。
ref($catch_seq);
在这个例子中该函数会返回Bio::SeqIO::genbank,所以说明\$catch_seq是一个Bio::SeqIO类型的对象!如果$catch_seq只是一个普通的标量的话,不会有返回值。
可以接着往下说了,因为GenBank文件的信息很丰富,所以在绝大多数情况下从数据库中下载得到一个Genbank的文件一般都只有一条序列,所以我们还可以偷懒一些,把文件的读取和序列的读取合并起来:
$seq_obj = Bio::SeqIO -> new(-file => 'ecorho.gbk', -format => 'genbank') ->next_seq;
一行语句中出现两个瘦箭头,不习惯吧?但这确实是Perl高手们很习惯的写法。注意:如果一个文件中有多条序列,那是不能这么写的。请老老实实按照上一篇所给的格式写。
最后的问题是,-format参数该怎么写?读取fasta文件时我们写-format => ‘fasta’,读取Genbank文件时我们写-format => ‘genbank’,那么其它的格式呢?比如SwissProt格式,是不是写-format => ‘swissprot’呢?不是的,其实应该写成-format => ‘swiss’,关于Bio::SeqIO模块究竟支持哪些序列格式,以及格式的写法,请到http://www.bioperl.org/wiki/HOWTO:SeqIO看一看。
如何获取序列的信息?
怎样从$seq_obj对象里把genbank序列的信息提取出来呢?我们还是依葫芦画瓢,照着上一篇提到的几项内容,代进去试试看:
$display_name = $seq_obj->display_name; # 序列的名称 $desc = $seq_obj->desc; # 序列的描述 $seq = $seq_obj->seq; # 序列字符串 $seq_type = $seq_obj->alphabet; # 序列的类型(dna还是蛋白质?) $seq_length = $seq_obj->length; # 序列的长度
你就会愉快地发现它们都能返回一些有趣的信息,$seq_obj能调用的方法还有一些,在http://www.bioperl.org/wiki/HOWTO:Beginners的Table3中有详细的描述。下面这张图列举了调用$seq_obj的某些方法之后能返回什么。
