Reverse Complement A Genbank File

Question

[Lee Katz]

Hi all, is there an easy way to reverse-complement a genbank file? I tried it with bioperl, and all the features disappeared.

I am trying to revcom certain contigs to visualize a region in a genome, across several genomes, so that I can view them all at once in MAUVE, and some contigs need to be revcom'd.

[Zev.Kronenberg]

You don't want to write a script yourself right? This pdf http://www.bsi.umn.edu/resources/perl2.pdf has a nice demo, if you are willing to script it up.

[Lee Katz]

I can script it up if there's help! I'm not sure if the tutorial you attached addresses features though.

[Zev.Kronenberg]

No it doesn't. Were you throwing any errors in bioperl? In my limited bioperl experience it is usually something trivial that screws up the script.

[Lee Katz]

Doing a richseq->revcom removes the features. Or maybe it's the truncation step. But, I think I got the answer on my own.

Answer 1

[Lee Katz]

I think I basically got it right, after I thought about the math. Here is the relevant snippet of code. There are certain variables that are defined somewhere else like $start and $end (start/end of the contig), but the idea is there.

  my @feat=$targetSeq->get_SeqFeatures();
  my $truncSeq=$targetSeq->trunc($start,$end);
  $truncSeq=$truncSeq->revcom if($revcom);

  for my $feat(@feat){
    my $subseqFeatStart=$feat->location->start-$start+1;
    my $subseqFeatEnd=$feat->location->end-$start+1;
    my $strand=$feat->strand;
    next if( ($subseqFeatStart<1 && $subseqFeatEnd<1) || ($subseqFeatStart>$seqLength && $subseqFeatEnd>$seqLength) );

    if($revcom){
      $subseqFeatStart=$end-($subseqFeatStart-$start)+1;
      $subseqFeatEnd=$end-($subseqFeatEnd-$start)+1;
      ($subseqFeatStart,$subseqFeatEnd)=($subseqFeatEnd,$subseqFeatStart);
      $strand=$strand*-1;
    }

    next if(ref($feat->location) ne "Bio::Location::Simple");
    $subseqFeatStart=1 if($subseqFeatStart<1);
    $subseqFeatEnd=$seqLength if($subseqFeatEnd>$seqLength);

    $feat->location->start($subseqFeatStart);
    $feat->location->end($subseqFeatEnd);
    $feat->strand($strand);
    $truncSeq->add_SeqFeature($feat);
  }

[Zev.Kronenberg]

up vote for self answer :)