bioperl 格式化genebank的输出

Posted 2020-09-26 庐州月光

代码如下：

use Bio::SeqIO;
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;

my $db_obj  = Bio::DB::GenBank->new;
my $seq_obj = $db_obj->get_Seq_by_acc(‘JN093905‘);
my $id      = $seq_obj->display_id();
my $organ;
foreach my $feat_object ($seq_obj->get_SeqFeatures) {

        if ($feat_object->primary_tag eq "source") {
            ($organ) = $feat_object->get_tag_values(‘organism‘);
        }

        if ($feat_object->primary_tag eq "CDS") {
            my ($gene_name) =  $feat_object->get_tag_values(‘gene‘);
            next if $gene_name ne ‘nifH‘;
            my ($seq)  =  $feat_object->get_tag_values(‘translation‘);
            my ($pro)  =  $feat_object->get_tag_values(‘product‘);
            $seq = lc($seq);
            my $len = (length($seq) + 1) * 3;
            print qq{>$id coded_by=<1..>$len,organism=$organ,definition=$pro\n$seq\n};
        }
}

运行结果如下：

>JN093905 coded_by=<1..>330,organism=uncultured Trichodesmium sp.,definition=dinitrogenase reductase
rlilnakaqttvlhvaaergavedveldevlkpgfggikcvesggpepgvgcagrgiitainfleeegaytdldfvsydvlgdvvcggfampirenkaqeiyivcsgem