This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Reverse Complement A Genbank File

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.

genbank
• 6,822 views • link
written by Lee Katz Lee Katz

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

• 0 views • link
written by Lee Katz Lee Katz

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.

• 0 views • link
written by Zev.Kronenberg Zev.Kronenberg

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

• 0 views • link
written by Lee Katz Lee Katz

1 answer

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);
  }
• 0 views • link
written by Lee Katz Lee Katz

Log in to answer this question.

More posts like this