利用bioperl读取复杂序列

Genbank序列描述的内容就非常丰富(类似的还有SwissProt,EMBL等格式的序列),除了名称、描述和序列字符串以外还有序列号、形状(线状还是环状?)、发布日期、所属物种以及序列内包含的基因、RNA、CDS、内含子(如果有的话)……What's more,这些信息并不是你想怎么写就怎么写的,而都是有一定的规范格式,故BioPerl可以准确地识别并提取这些信息。所以,不要自行修改下载得到的Genbank序列,否则BioPerl可能会无法识别。

如何从本地文件中提取Genbank序列?

[perl]
#!usr/bin/perl -w

use Bio::SeqIO;
$catch_seqio_obj->new(-file => 'wrky.gbk', -format => 'genbank');
$seq_obj = $catch_seqio_obj->next_seq; #后缀加上去比较好,也可以使用ref函数来检查一下$catch_seqio_obj究竟是什么。ref($catch_seq); 该函数会返回Bio::SeqIO::genpept,所以说明$catch_seqio_obj是一个Bio::SeqIO类型的对象!
#$seq_obj = Bio::SeqIO -> new(-file => 'wrky.gbk', -format => 'genbank') -> next_seq; 上边两句合并,注意:如果一个文件中有多条序列,那是不能这么写的。
[/perl]
[perl]
$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; # 序列的长度
[/perl]
$seq_obj能调用的方法还有一些,在http://www.bioperl.org/wiki/HOWTO:Beginners的Table3中有详细的描述。

我们更想要物种名以及序列里面包含的基因、mRNA、CDS等元件的位置、名称、序列等信息。BioPerl有另外一套机制来获取它们。

处理多个Genbank文件:只是模块要换一下,把Bio::SeqIO改成Bio::SeqIO::MultiFile,
[perl]
use Bio::SeqIO::MultiFile;
$catch_seqio_obj = Bio::SeqIO::MultiFile -> new(-files=>['a.gbk','b.gbk'], -format=>'genbank'); # 读取a.gbk和b.gbk两个文件
while($seq_obj = $catch_seqio_obj -> next_seq) {
do sth...
}
[/perl]
或者:
[perl]
use Bio::SeqIO::MultiFile;
$catch_seqio_obj = Bio::SeqIO::MultiFile -> new(-files=>[glob "*.gbk"], -format=>'genbank'); # 读取当前目录所有后缀名为gbk的文件。
while($seq_obj = $catch_seqio_obj-> next_seq) {
do sth...
}
[/perl]
从命令行中传入参数

如果我们想用前面的代码处理另一个Genbank文件(比如seq.gbk),我们必须打开源代码文件,把-file => 'wrky.gbk' 改成 -file => 'seq.gbk',这样子太麻烦了。事实上这些要处理的文件的名字应该由用户来决定,而不是由你写在源代码里。有没有办法像“钻石操作符”(<>)那样通过命令行参数读入文件名呢?当然可以请看:
[perl]
use Bio::SeqIO;
$filename = shift @ARGV; # 从命令行中读取文件名,而不是写在程序里
$catch_seqio_obj = Bio::SeqIO -> new(-file => $filename, -format => 'genbank');
[/perl]
呵呵,这样子这几行代码就可以处理世界上所有的Genbank文件了,再也不需要修改源代码了。

读取多个文件
[perl]
use Bio::SeqIO::MultiFile;
@filenames = @ARGV;
$catch_seq = Bio::SeqIO::MultiFile -> new(-files => [@filenames], -format => 'genbank');
[/perl]

 

Table 3: Values from the Sequence object (Genbank)
MethodReturns
display_idECORHO
descE.coli rho gene coding for transcription termination factor.
display_nameECORHO
accessionJ01673
primary_id147605
seq_version1
keywordsattenuator; leader peptide; rho gene; transcription terminator
is_circular
namespace
authority
length1880
seqAACCCT...ACAGGAC
divisionBCT
moleculeDNA
get_dates26-APR-1993
get_secondary_accessionsJ01674

 

关于Bio::SeqIO模块究竟支持哪些序列格式,以及格式的写法如下:(from http://www.bioperl.org/wiki/HOWTO:SeqIO

 

Table 1: Bio::SeqIO modules and formats supported
NameDescriptionFile extensionModule
abiABI tracefileab[i1]Bio::SeqIO::abi
aceAce databaseaceBio::SeqIO::ace
agaveAGAVE XMLBio::SeqIO::agave
alfALF tracefilealfBio::SeqIO::alf
asciitreewrite-only, to visualize featuresBio::SeqIO::asciitree
bsmlBSML, using XML::DOMbsmlBio::SeqIO::bsml
bsml_saxBSML, using XML::SAXBio::SeqIO::bsml_sax
chadoxmlCHADO sequence formatBio::SeqIO::chadoxml
chaosCHAOS sequence formatBio::SeqIO::chaos
chaosxmlChaos XMLBio::SeqIO::chaosxml
ctfCTF tracefilectfBio::SeqIO::ctf
emblEMBL databaseebl|emb|datBio::SeqIO::embl
entrezgeneEntrez Gene ASN1Bio::SeqIO::entrezgene
excelExcelBio::SeqIO::excel
expStaden EXP formatexpBio::SeqIO::exp
fastaFASTAfast|seq|fa|fsa|nt|aaBio::SeqIO::fasta
fastqquality score data in FASTA-like formatfastqBio::SeqIO::fastq
flybase_chadoxmlvariant of Chado XMLBio::SeqIO::flybase_chadoxml
gameGAME XMLBio::SeqIO::game
gcgGCGgcgBio::SeqIO::gcg
genbankGenBankgbank|genbankBio::SeqIO::genbank
interproInterProScan XMLBio::SeqIO::interpro
keggKEGGBio::SeqIO::kegg
largefastaLarge files, fasta formatBio::SeqIO::largefasta
lasergeneLasergene formatBio::SeqIO::lasergene
locuslinkLocusLink LL_tmplBio::SeqIO::locuslink
metafastaBio::SeqIO::metafasta
phdPhredphredBio::SeqIO::phd
pirPIR databasepirBio::SeqIO::pir
plnPLN tracefileplnBio::SeqIO::pln
qualPhredBio::SeqIO::qual
rawplain texttxtBio::SeqIO::raw
scfStandard Chromatogram FormatscfBio::SeqIO::scf
seqxmlSeqXML sequence format using XML::LibXML and XML::WriterxmlBio::SeqIO::seqxml
striderDNA Strider formatBio::SeqIO::strider
swissSwissProtspBio::SeqIO::swiss
tabtab-delimitedBio::SeqIO::tab
tableTableBio::SeqIO::table
tigrTIGR XMLBio::SeqIO::tigr
tigrxmlTIGR Coordset XMLBio::SeqIO::tigrxml
tinyseqNCBI TinySeq XMLBio::SeqIO::tinyseq
ztrZTR tracefileztrBio::SeqIO::ztr

 
本文来自:http://blog.163.com/guojinyan2008@126/blog/static/3074706720111019101034907/

  • 本文由 整理发表
  • 网站部分文章源自互联网,若未正确标注来源,请联系管理员更新。文章转载,请务必保留本文链接

您必须才能发表评论!