2
\$\begingroup\$

I have written a program to split a Multi-FASTA file into individual fasta files.

use strict;
use autodie;

open( my $IN, '<', $ARGV[0] );

my ( $id, $seq_obj );

sub dump_seq {

    my $sid = shift;
    my $seq = shift;

    my @parts = split( /\h+/, $sid );
    my $f_name = $parts[0];

    $f_name =~ s/\|/_/g;
    $f_name =~ s/^>//;
    $f_name =~ m/\_$/
        ? ( $f_name =~ s/.$/\.fasta/ )
        : ( $f_name .= ".fasta" );  # Controlling extension
    my @lines = $seq =~ m/.{1,60}/g; # controlling line length

    open( my $out, '>', $f_name );

    print {$out} "$sid\n";
    print {$out} "$_\n" for (@lines);

    return 0;

}

while (<$IN>) {

    chomp;

    if (/^>/) {

        if ( $. > 1 ) { # Else there would be no-output for first fasta
            dump_seq( $id, $seq_obj );
        }

        $id = $_;
        next;
    }

    $seq_obj = $seq_obj . $_;

    dump_seq( $id, $seq_obj ) if (eof); # For last fasta only
}

In this code I made calls to dump_seq() two times, I want to minimize it to one and place the call appropriately.

I can also do it by storing all the $id and $seq_obj in a hash and call the function once, which takes up lot of memory, which I don't want.

Sample Input:

>gi|129367|sp|P13481.1|P53_CHLAE RecName: Full=Cellular tumor antigen p53; AltName: Full=Tumor suppressor p53 [Chlorocebus aethiops] 
MEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEAAPHMAPTPAAP
KPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPAGSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSD
>gi|635087583|ref|XP_008008385.1| PREDICTED: cellular tumor antigen p53 isoform X2 [Chlorocebus sabaeus]
MEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEAAPPMAPTPAAP
SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENFRKKGEPCHELPPGSTKRALPNNTSSSPQPKK
KPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPAGSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSD
>gi|402898625|ref|XP_003912321.1| PREDICTED: cellular tumor antigen p53 [Papio anubis]
MEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEATPPMAPTPAAP
KPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPAGSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSD
\$\endgroup\$
  • \$\begingroup\$ Do you need to control the line-length like you do, or is the input-file line length good enough? \$\endgroup\$ – rolfl Oct 22 '14 at 14:34
  • \$\begingroup\$ Controlling line length makes the file look pretty and is also standard. \$\endgroup\$ – VAR121 Oct 22 '14 at 14:49
  • \$\begingroup\$ Somewhat frustrated by your question... your code produces the wrong output. \$\endgroup\$ – rolfl Oct 22 '14 at 15:14
  • \$\begingroup\$ Where are you facing problem ? It was working fine with me. \$\endgroup\$ – VAR121 Oct 22 '14 at 15:29
  • 2
    \$\begingroup\$ Sorry mistake while copying the code. \$\endgroup\$ – VAR121 Oct 22 '14 at 15:32
2
\$\begingroup\$

There are a number of things in here which make the code somewhat hard to read. Also, there's a big bug.

Bug

You never reset the $seq_obj variable. As a result, each time you add a new file you are actually adding the new data to the end of the old data too. Thus, your files get bigger, and bigger, etc.

Thus, the third sequence in your input:

>gi|402898625|ref|XP_003912321.1| PREDICTED: cellular tumor antigen p53 [Papio anubis]
MEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEATPPMAPTPAAP
KPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPAGSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSD

actually produces the output:

>gi|402898625|ref|XP_003912321.1| PREDICTED: cellular tumor antigen p53 [Papio anubis]
MEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGP
DEAPRMSEAAPHMAPTPAAPKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPA
GSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSDMEEPQSDPSIEPPLSQETFSDLWKLLP
ENNVLSPLPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEAAPPMAPTPAAPSCMGGMN
RRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENFRKKGEPCHELPPGSTKRA
LPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPAGSRAHSS
HLKSKKGQSTSRHKKFMFKTEGPDSDMEEPQSDPSIEPPLSQETFSDLWKLLPENNVLSP
LPSQAVDDLMLSPDDLAQWLTEDPGPDEAPRMSEATPPMAPTPAAPKPLDGEYFTLQIRG
RERFEMFRELNEALELKDAQAGKEPAGSRAHSSHLKSKKGQSTSRHKKFMFKTEGPDSD

which is the accumulation of all the previous inputs too. I have spent 20 minutes trying to figure out why my code is wrong.... but it's your code.


Better system

Additionally, the process you have 'batches' the data in one variable, and then dumps the variable when all the data is read.

If you reverse the algorithm, and open a file when the header arrives, and then write all the data (that's not a header) to that file, then the process would be simpler.

Additionally, you should use the 'diamond' magic filehandle <> for ease of use:

The null filehandle <> is special: it can be used to emulate the behavior of sed and awk, and any other Unix filter program that takes a list of filenames, doing the same to each line of input from all of them. Input from <> comes either from standard input, or from each file listed on the command line.

So, if I break your code in to three sections. One section closes an open file, the next section starts a new file, and then a simple loop that checks to see whether it should start a new file, the code is much simpler.

Note also the general lack of global variables. Also, the 'id' line is not printed to the file when the file is created, but rather after the file is created as part of the 'normal' loop.

Finally, the process below copies the input lines unchanged from the input to the output. Your code 'reformats' the lines to 56 chars wide. This code does not do the reformat. I consider this to be a feature because of two things:

  1. preserving input format is generally a good idea.
  2. it saves having to chomp every line
  3. it allows the data to be processed 1 line at a time, instead of batching it

If that does not work for you, then the code will have to be changed significantly.

#!/usr/bin/perl

use strict;
use autodie;

sub flushfile {
    my $file = shift;
    return undef unless $file;
    close $file;
    return undef;
}

sub makefile {
    my $control = shift;
    chomp $control;


    my @parts = split( /\h+/, $control, 2);
    my $f_name = $parts[0];
    $f_name =~ s/\|/_/g;
    $f_name =~ s/^>//;
    $f_name =~ s/\_$//;
    $f_name .= ".fasta";

    print "Processing $f_name from $control\n";

    open (my $ofile, '>', $f_name);
    return $ofile;
}

my $outfile;
while (<>) {

    if (/^>/) {
        $outfile = flushfile($outfile);
        $outfile = makefile($_);
    }

    $outfile or die "Unexpected input: $_";
    print {$outfile} $_;
}
flushfile($outfile);

Line Trimming version

I have written it again using a system that does the line sizing too.

The following code is different. It accumulates partial lines and the outputs them when they are full. It also needs to print the control line first because that cannot be line-trimmed.

#!/usr/bin/perl

use strict;
use autodie;

my $linesize = 60;

sub flushfile {
    my $file = shift;
    my $partial = shift;
    return undef unless $file;
    print {$file} $partial . "\n";
    close $file;
    return undef;
}

sub makefile {
    my $control = shift;
    chomp $control;

    my @parts = split( /\h+/, $control, 2);
    my $f_name = $parts[0];
    $f_name =~ s/\|/_/g;
    $f_name =~ s/^>//;
    $f_name =~ s/\_$//;
    $f_name .= ".fasta";

    print "Processing $f_name from $control\n";
    open (my $ofile, '>', $f_name);
    print {$ofile} $control . "\n";
    return $ofile;
}

sub addline {
    my $file = shift;
    my $part = shift;
    my $nxt = shift;
    $file or die "Unexpected input: $nxt";

    $part .= $nxt;
    while (length($part) >= $linesize) {
        my $front = substr ($part, 0, $linesize);
        #print "Printing $front\n";
        print {$file} ($front , "\n");
        $part = substr ($part, $linesize);
    }
    return $part;
}

my ($partial, $outfile);
$partial = "";
while (<>) {

    chomp;

    if (/^>/) {
        $outfile = flushfile($outfile, $partial);
        $outfile = makefile($_);
    } else {
        $partial = addline($outfile, $partial, $_);
    }
}
flushfile($outfile, $partial);
\$\endgroup\$
  • \$\begingroup\$ Batching of all lines of a record was specially done to add a control filter later. Like if I need sequences if(length > N). \$\endgroup\$ – VAR121 Oct 22 '14 at 15:14
  • \$\begingroup\$ Editing my answer to list the bugs in your code... give me a sec. \$\endgroup\$ – rolfl Oct 22 '14 at 15:14
  • \$\begingroup\$ My Mistake, It was a copy error, I was copying through x2vnc from one computer to other and the source was vim, I missed one line while scrolling. \$\endgroup\$ – VAR121 Oct 22 '14 at 15:38
1
\$\begingroup\$

I'd probably drop some newlines and parentheses, but otherwise, looks ok.

I often encounter a similar situation: after processing input line by line, I have to call a processing sub once more to process the last entry. To simplify the syntax, you can remove the if eof part and move the line after the loop.

I also don't quite like the $. > 1 condition. Wouldn't if $id work the same but give more clues?

\$\endgroup\$
  • \$\begingroup\$ $. > 1 can be substituted by ($id), and might some sense in to the construct. \$\endgroup\$ – VAR121 Oct 22 '14 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.