trimming FASTA file

[baika]

Hi All
I have a multi-fasta file with approximately 2000 sequences of varying length that have different start and end regions. I need to trim all my sequences in a way that all the sequences start with “ATAGCCGGCACCCTGGT” and ends with “GGCCATATGAGTGGGCC”. Any script would be really helpful to remove bases upstream of ATAGCCGGCACCCTGGT and downstream of GGCCATATGAGTGGGCC so that all the sequences become of equal length, and have the same start and end.

Thanks for your help

Baika

[krobison]

#!/usr/bin/perl
use strict;
use Bio::SeqIO;

## debugging/tuning left as exercise for student :-)

my $reader=new Bio::SeqIO(-format=>'fasta',-file=>$ARGV[0]);
my $writer=new Bio::SeqIO(-format=>'fasta',-file=>$ARGV[0]."trimmed");
while (my $rec=$reader->next_seq)
{
  ## works only on forward strand & 0 mismatches!!!!
   if ($seq->seq=~/(ATAGCCGGCACCCTGGT.*GGCCATATGAGTGGGCC)/i)
   {
   $rec->seq($1);  
   $writer->write_seq($rec);
  } 
  else
  {
     print STDERR "Could not find head and/or tail sequences for ",$seq->id,"\n";
   }
}

[baika]

Thanks Krobison for the perl script. It gives an error-
Global symbol "$seq" requires explicit package name at ../../Scripts/trim_fastaseq.pl line 12.
Global symbol "$seq" requires explicit package name at ../../Scripts/trim_fastaseq.pl line 19.

Baika

[kmcarr]

Sorry baika; I was going to answer until I saw this line in Keith's code:

## debugging/tuning left as exercise for student :-)

Keith, did you stick the bug in there on purpose?

[A_Morozov]

No fun, guys. Baika might be (and claims to be in introduction section) a non-bioniformatitian stuck with what is not his/her area of expertise. If you went to genomics/bioinformatics lab, you should at least learn some Perl/Python, but for now - seems like it should be $rec, not $seq under regexp (the scary thingy with // and lots of uppercase).

[kmcarr]

Originally Posted by A_Morozov

No fun, guys. Baika might be (and claims to be in introduction section) a non-bioniformatitian stuck with what is not his/her area of expertise. If you went to genomics/bioinformatics lab, you should at least learn some Perl/Python, but for now - seems like it should be $rec, not $seq under regexp (the scary thingy with // and lots of uppercase).

Oh, it was late and I was feeling a tad impish; I wasn't going to leave baika hanging long. Anyway the better solution is to change line #9, naming the first object $seq.

Change

while (my $rec=$reader->next_seq)

to

while (my $seq=$reader->next_seq)

[d1antho]

Not that it should make a difference, because it is unlikely to match in the middle of a sequence, but the the match operator should be bounded to the start and the end of the sequence read. But honestly I doubt it should make any difference at all

So this

if ($seq->seq=~/(ATAGCCGGCACCCTGGT.*GGCCATATGAGTGGGCC)/i)
   {
   $rec->seq($1);  
   $writer->write_seq($rec);
  }

Should then be

if ($seq->seq=~/(^ATAGCCGGCACCCTGGT.*GGCCATATGAGTGGGCC$)/i)
   {
   $rec->seq($1);  
   $writer->write_seq($rec);
  }

[baika]

Thanks krobison for writing this script, and kmcarr for pointing out the error. After incorporating all your suggestions and help from my friend Robert, finally it is working.

Thank you all

baika

#!/usr/bin/perl -w

#Usage: trim_fasta.pl YOUR_FASTA_FILE.fasta OUT_FILE_TRIMMED.FASTA

use strict;
use Bio::SeqIO;

my $reader=new Bio::SeqIO(-format=>'fasta',-file=>$ARGV[0]);
my $writer=new Bio::SeqIO(-format=>'fasta',-file=>">$ARGV[1]");
while (my $seq=$reader->next_seq)
{
  ## works only on forward strand & 0 mismatches!!!!
   if ($seq->seq=~/(CCAGTATTTGGTA.*AGTTGATAACTGGGAA)/i)
   {
   $seq->seq($1);  
   $writer->write_seq($seq);
  } 
  else
  {
     print STDERR "Could not find head and/or tail sequences for ",$seq->id,"\n";
   }
}