博文
Perl代码实例分析之从下载序列·一
|||
从网络上下载DNA，mRNA或protein序列，是分子生物学和生物信息学的基本功。最简单的下载序列的方式莫过于直接从NCBI (http://www.ncbi.nlm.nih.gov/) 检索到序列，然后通过网页上的Send to选项下载之。本文的目的是通过对一例序列下载Perl语言代码来简要介绍Perl及BioPerl的一些基础知识。
Perl编程语言是生物信息学最流行的动态解释计算机语言(Dynamic Interpreted Computer Languages)。另外常用的动态解释语言有Ruby和Python。但是，Perl相较于其它动态解释语言最大的优势在于拥有一个包含大量标准模块的仓库(CPAN, http://www.cpan.org/)，可供开发者使用。另外，Perl语言本身的特点是，不同程度学习者都可以用来解决问题，比较适合研究生在生物信息学工作中通过写较少行数的Perl代码解决一些小问题。另外，Perl语言也被称为编程语言中的“瑞士军刀”，功能强大。
这里介绍的Perl代码实例使用了Bio::Perl模块来下载序列，并将序列写入文件，Bioperl(http://search.cpan.org/~cjfields/BioPerl-1.6.901/Bio/Perl.pm, http://www.bioperl.org/wiki/Main_Page )是一个专门为生物信息学应用而开发的一个Perl模块套件。
写程序必须要考虑数据的输入和输出，因为常用的生物信息学平台为GNU/Linux操作系统，命令行状态，所以，这里没有考虑Windows操作系统和图形用户界面(GUI)。这里使用三个参数跟在程序名后作为数据输入的方式
第一个参数是，来源数据库，支持的有swiss, embl, genbank, genpept, refseq。
第二个参数是，输出文件的格式，可以是fasta和genbank。
第三个参数是，序列编号文件，为一个编号一行的纯文本文件。
确定了以上三个参数后，还要进行环境的判断，假如参数没有三个怎么办呢？这里要说明Perl的一个默认数组@ARGV，用来保存程序的参数，因此，运行程序时，程序名(默认变量$0)后的跟的第一个参数就是$ARGV[0]，相应的第二个是$ARGV[1]，第三个是$ARGV[2]。
判断运行环境时首先定义了一个子过程，功能是打印使用方法（&usage，如果子程序定义在前，&符号可以省略），然后退出。接着进行了以下判断。
如果参数不是三个，调用usage子过程，退出程序。
如果第一个参数不是swiss, embl, genbank, genpept, refseq之一，调用usage。
如果第二个参数不是fasta, genbank之一，调用usage。
如果第三个参数指定文件打不开，调用usage。
通过以上对参数的判断后，就可以按照预想的方式下载序列了。从参数中，程序了解到从哪个数据库下载哪些序列，保存为何种格式的文件。这样，程序就通过循环分别下载序列名称列表中的每一条序列，并写入到相应序列名称的文件中。
例子：
echo NG_032174.1 JQ180353.1 | perl program_name genbank fasta -
perl program_name genbank fasta List.txt
程序中涉及的关键内容
Bio::Perl
@ARGV
子程序
STDERR
get_sequence(‘数据库’,‘序列编号’)
write_squence(“>文件名”，“文件格式”，序列对象)
对象
源代码：
#!/usr/bin/env perl
use Bio::Perl;
sub usage{
print STDERR "Usage: $0 <swiss|embl|genbank|genpept|refseq> <fasta|genbank> <List_of_Acc>n";
exit;
}
unless(@ARGV==3){&usage;}
$db=$ARGV[0];
unless($db == 'swiss' || \$db == 'embl' || \$db == 'genbank' || \$db == 'genpept' || $db == 'refseq'){
print STDERR "Database name error: $dbn";
usage;
}
$file_type=$ARGV[1];
unless($file_type == 'fasta' || $file_type == 'genbank'){
print STDERR "File type error: $file_typen";
usage;
}
$file=$ARGV[2];
unless(open FILE,$file){
print STDERR "Open file error: $filen";
usage;
}
while(my $acc = <FILE>){
$acc =~ s/[rn]//g; #去掉换行符号
print STDERR "Downloading $acc ... ";
$seq=get_sequence('genbank',"$acc");
if(write_sequence(">./$acc.\$ARGV[1]", "$ARGV[1]",$seq)){
print STDERR "Writing to file $acc.$ARGV[1]";
}
print STDERR "n";
}
https://blog.sciencenet.cn/blog-656335-563847.html
下一篇：Perl代码实例之词频统计软件word_freq
