simple perl script to trim a FASTA file

[RabidMortal]

hello,

i need help finding/writing a simple script to edit a very large FASTA (text) file.

the text format of the FASTA file is simple:

>HWI-EAS158_40_3_1_46_535
GTGAATGCGTGATACAGGAATGTTCGTTGTGACCAT
>HWI-EAS158_40_3_1_47_579
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
>HWI-EAS158_40_3_1_46_731
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA

each file has 6000000 lines, all with the exact same format. i need a script that will go through and trim off x nucleotides from the beginning, and y nucleotides from the end, of each and every sequence in the file. so, the script should not touch anything on the lines beginning with ">"

i would really really appreciate any help.

thank you.

[BrowserUk]

Generally, FASTA files contain multiple sequence lines (wrapped at n? chars), per header line, which means your example could be misleading. If that is the case, then something like this may be needed:

#! perl -slw
use strict;

my( $x, $y ) = ( 3, 5 );

while( <DATA> ) {
    chomp;
    print ">$_";
    local $/ = '>';
    $_ = <DATA>;
    chomp;
    tr[\n][]d;
    print for unpack '(A36)*', substr $_, $x, length() - $x - $y;
}

__DATA__
>HWI-EAS158_40_3_1_46_535
GTGAATGCGTGATACAGGAATGTTCGTTGTGACCAT
>HWI-EAS158_40_3_1_47_579
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
>HWI-EAS158_40_3_1_46_731
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA

Produces:

c:\test>junk46
>>HWI-EAS158_40_3_1_46_535
AATGCGTGATACAGGAATGTTCGTTGTG
>HWI-EAS158_40_3_1_47_579
GTGAATGCGTGATACAGGAATGTTCGTTGTGACAAA
GTGAATGCGTGATACAGGAATGTTCGTTGTGACAAA
GTGAATGCGTGATACAGGAATGTTCGTTGTGACAAA
GTGAATGCGTGATACAGGAATGTTCGTT
>HWI-EAS158_40_3_1_46_731
TCATGCGTGATACAGGAATGTTCGTTGTGAAAAGTG
TCATGCGTGATACAGGAATGTTCGTTGTGAAAAGTG
TCATGCGTGATACAGGAATGTTCGTTGT

[almut]

trim off x nucleotides from the beginning, and y nucleotides from the end

#!/usr/bin/perl

while (<DATA>) {
    s/^(?!>)...(.*)..../$1/;
    #       ^^^    ^^^^  
    #       x=3    y=4

    print;
}

__DATA__
>HWI-EAS158_40_3_1_46_535
GTGAATGCGTGATACAGGAATGTTCGTTGTGACCAT
>HWI-EAS158_40_3_1_47_579
AAAGTGAATGCGTGATACAGGAATGTTCGTTGTGAC
>HWI-EAS158_40_3_1_46_731
GTGTCATGCGTGATACAGGAATGTTCGTTGTGAAAA

For longer stretches you can also write .{N} (with N being the number of nucleotides to trim).

[RabidMortal]

Thank you thank you.

So far, I have not been able to get either script to run. I keep getting a missing semicolon error--I am not sure what that means but that seems like something I can figure out on my own.

What I cannot figure out is if there an easy way I can get these scripts to run on a specified FASTA file? Again, my files are very very large (gigabytes) and I'd prefer not to edit them manually...

Again, I really appreciate the help.

[erix]

You should indeed not edit the sequence files.

The previous answers (by almut and BrowserUK) just use __DATA__ because that makes it easy to show as compact example. And to use the examples they'll need a little editing.

For instance, in almut's script, replace the line:

while( <DATA> ) {

with

while( <> ) {

and invoke the script as:

perl almuts_answer.pl sequences.fa > sequences_cleaned.fa

Replacing <DATA> with <> means it will not read the __DATA__ section of the program-file itself, but instead read from STDIN (aka standard input).

hth

[RabidMortal]

Thank you thank you thank you!!!

I really appreciate the help and both scripts work great.

This is a great way to end the week!

:)

[silver_steve]

Could someone please explain what the (A36) does in BrowserUk's code? Thank you.

[toolic]

perlpacktut