[Bioperl-l] some contigs do not work for sequence retrievel

Hermann Norpois hnorpois at googlemail.com
Tue May 8 16:16:28 UTC 2012


Hello,

for getting a sequence 5 prime upstream of TTS I wrote a script that works
for some geneids but not for all. I always get a contig and coordinates. I
do not have an idea why I do not get a sequence ( I only get fasta
headers). Actually the sequence ID should be out of importance if I see
that a contig is detected. Has anybody an idea?

Thanks
Hermann Norpois


#!/bin/perl -w
use strict;
use Bio::DB::EntrezGene;
use Bio::SeqIO;
use Bio::DB::GenBank;

my $id = "12064"; #Works with geneid 18619 (Penk1) but not with 54161
(copg) or 12064 (bdnf)

my $seqio_obj = Bio::SeqIO->new(-file => ">bdnf.fasta", -format => 'fasta'
);

my $db = new Bio::DB::EntrezGene;

my $seq = $db->get_Seq_by_id($id);

my $ac = $seq->annotation;

for my $ann ($ac->get_Annotations('dblink')) {
    if ($ann->database eq "Evidence Viewer") {
                # get the sequence identifier, the start, and the stop
        my ($contig,$from,$to) = $ann->url =~
          /contig=([^&]+).+from=(\d+)&to=(\d+)/;
                 my $chr_start = $from-700;
                 my $chr_stop = $from;
               #  my $strand = 1;
                print "CONTIG:\t$contig\tFROM\t$from\tTO\t$to\n\tFETCHING
SEQUENCE FROM\t$chr_start\tTO\t$chr_stop\n"; # Control that something was
detected.
                my $gb  = Bio::DB::GenBank->new(-format     => 'fasta',
                                -seq_start  => $chr_start,
                                -seq_stop   => $chr_stop,
                #               -strand     => $strand
                #               -complexity => 1
                                );
              #  $gb->request_format('fasta');
                my $obj = $gb->get_Seq_by_id($contig);

                $seqio_obj->write_seq($obj);

    }
}



More information about the Bioperl-l mailing list