[Bioperl-l] finding sequence 2000 bp upstream in bioperl?

Jason Stajich jason at bioperl.org
Wed Jul 2 19:57:15 UTC 2008


Right on Kevin - except that you have to call trunc() not subseq() if  
you want to get back a sequence object that can be written out with  
Bio::SeqIO and you'll have to do it slightly differently if the gene  
is on the opposite strand operating on the $f->end, $f->end + 2000  
and then doing a revcom.

I *think* this will get you what you want:
my $promotor;
if( $f->strand < 0 ) {
  $promotor = $entry->trunc($f->end+1, $f->end + 2000)->revcom;
} else {
  $promotor = $entry->trunc($f->start - 2000, $f->start -1);
}

-jason
On Jul 2, 2008, at 3:36 PM, Kevin Brown wrote:

> Get the start point for a Sequence Feature and request a subseq of the
> main object that starts 2000 above that.
>
> my $gb  = new Bio::DB::GenBank;
> $entry = $gb->get_Seq_by_id($id);
> foreach my $f ($entry->all_SeqFeatures())
> {
> 	if ($f->primary_tag eq 'CDS')
> 	{
> 		my $seq = $entry->subseq($f->start - 2000,$f->start -
> 1);
> 		$out->write_seq($seq);
> 	}
> }
>
>> -----Original Message-----
>> From: bioperl-l-bounces at lists.open-bio.org
>> [mailto:bioperl-l-bounces at lists.open-bio.org] On Behalf Of
>> RICHARD GREEN
>> Sent: Wednesday, July 02, 2008 11:10 AM
>> To: bioperl-l at lists.open-bio.org
>> Subject: [Bioperl-l] finding sequence 2000 bp upstream in bioperl?
>>
>> Howdy Bioperlers,
>>
>> Quick question, I have scripts that pull out a sequence from
>> genbank, but I am looking for a command that would extract
>> 2000 base pairs upstream of the sequence. Is there a way to
>> extract just the Gene promoter region from genbank? Any
>> advice that folks can give me a is muchly appreciated.
>>
>> Thanks again
>>
>> -Rich Green
>>
>>
>>
>>
>> _______________________________________________
>> Bioperl-l mailing list
>> Bioperl-l at lists.open-bio.org
>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>>
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l




More information about the Bioperl-l mailing list