sub parse_genbank {
my $filename = shift;
open GBFH, $filename;
my ( $fflag, $features, $flag, $seqflag ) = ( 0, 0, 0, 0 );
my ( $seq, $head, $id, $acc, $cds );
my $data = {};
foreach ( <GBFH> ) {
if ( /^LOCUS/ && $flag == 0 ) {
/^LOCUS\s*(.*\S)\s*$/;
( $id = ( $head = $1 ) ) =~ s/\s.*//;
$features += 1;
$flag = 1;
}
elsif ( /^DEFINITION\s*(.*)/ && $flag == 1 ) {
$head .= " $1";
$features += 2;
}
elsif ( /^ACCESSION/ && $flag == 1 ) {
/^ACCESSION\s*(.*\S)\s*$/;
$head .= " " . ( $acc = $1 );
$features += 4;
}
elsif ( /^FEATURES/ && $flag == 1 ){
$fflag = 1;
}
elsif ( /^ {5}CDS/ && $fflag == 1){
/^ {5}CDS\s*(.*\S)\s*$/;
$cds = $1;
}
elsif ( /^ORIGIN/ ) {
$seqflag = 1;
$features += 8;
}
elsif (/^\/\//) { # end of genbank file
$features += 16;
}
elsif ( $seqflag == 0 ) {
;
}
elsif ( $seqflag == 1 ) {
$seq .= $_;
}
if ( $features == 31 ) {
$seq =~ tr/a-zA-Z//cd;
$seq =~ tr/A-Z/a-z/;
my @pos = split /\.+/, $cds;
my ($begin,$end) = @pos;
$data->{_sequence} = $seq;
$data->{_header} = $head;
$data->{_id} = $id;
$data->{_accession} = $acc;
$data->{_cds} = substr($seq,$begin-1,$end-$begin+1);
return \$data;
}
}
return;
}
2009/4/23 guo <sczz520@126.com>
我在windows上下载安装了activeperl,想编写给程序把人类染色体序列上的编码区也就是CDS给提取出来。有朋友说用bioperl模块
来编程效率更高。可是我不会用。自己试着使用ppm来安装bioperl也没有成功。
现在主要是想请教一下,有没有人会bioperl,或者activperl的,能不能帮我写一个程序来把人类染色体序列中的CDS提取出来。
染色体序列我已经从NCBI上下载下来了。就是不知道如何编写程序。请教高手。
即使程序付费,我也愿意。
希望广大perl友帮忙。
--~--~---------~--~----~------------~-------~--~----~
您收到此信息是由于您订阅了 Google 论坛"PerlChina Mongers 讨论组"论坛。
要在此论坛发帖,请发电子邮件到 perlchina@googlegroups.com
要退订此论坛,请发邮件至 perlchina+unsubscribe@googlegroups.com
更多选项,请通过 http://groups.google.com/group/perlchina?hl=zh-CN 访问该论坛
-~----------~----~----~----~------~----~------~--~---
没有评论:
发表评论