Hi monks, I have a list of 121,000 fasta headers in a .txt file and I want to compare this to a multiple FASTA file containing 132,000 sequences. I have written a script to do this but it is too slow!!! The script puts all the lines from the .txt file in a hash and then compares the header of each FASTA sequence in the multiple FASTA file with the keys in the hash and if they match prints out a new FASTA sequence with the renamed header from the .txt file. Is there anything sensible I can do to speed this up?

Input txt file looks like: 

F0Z7V0F01A03EB_210
F0Z7V0F01A03EB_180
F0Z7V0F01A03EB_136
F0Z7V0F01A03EB_362
F0Z7V0F01A03GP_298
F0Z7V0F01A03GP_205
etc...........

[download]

Input FASTA looks like:

>GJKKTUG01DYDGC
GGGTATTCCTTCTCCACCTTGCAGCTAACATCAGTGTTTCGTCTACTCAAGCACGCCAAC
ACGCCCTAGAGCGCCCTGTCCAGGGGATGGCAACCAACTCTGACCCTGCAAGTGCAGCAG
ACATGAGGAATACAAACTACAATCTTTTACTTGATGATGCAATGCCGGACAAACTCTAGA
>F0Z7V0F01EDB3V
AAGGCGAGNGGTATCACGCAGTAAGTTACGGTTTTCGGGTAACGCGTCNGNGGNACTAAC
CCACGGNGGGTAACCCGTCNCTACCGGTATAGGACTAAGGTTACCGGAACGTCGTGGGGT
ACCCCCCGGACGGGGACCGTCCCCTCATANAGTCAACNGTNTGAGATGGACTAACTCAAA
CCTAGTTTCAAGTACTATTTAACTTACTTACGTTACCCGTAATTTCGGCGTTTAGAGGCG

etc....................

[download]

My script which is working VERY slowly is :

#!/usr/bin/perl

use strict;
use warnings;
use Bio::SeqIO;
use Data::Dumper;
my %seq_id;
my $fasta_id;
        open HEADER , "<FASTA.headers" or die $!;
        while (<HEADER>){
            chomp $_;
            #print $_."\n";
            $fasta_id = $_; 
            $fasta_id =~ s/_.*//g ;
            #print $fasta_id."\n";                
            %seq_id = ("$fasta_id" => "$_");
#            print Dumper (\%seq_id);


            my $infile=$ARGV[0] || die ("Please give me an input fasta
+ file\n");
            my $inseq = new Bio::SeqIO(-format => 'fasta',
                        -file => $infile);
            while (my $seq_obj = $inseq->next_seq ) {
             my $id = $seq_obj->id ;
            chomp $id;
#            print $id."\n";
            my $seq = $seq_obj->seq ;
                if (exists ($seq_id{$id})) {
                print ">";
                print $seq_id{$fasta_id};
                print "\n".$seq."\n";
                } 
                

            }
}

[download]

In reply to renaming 1000's of FASTA files by garyboyd

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post, it's "PerlMonks-approved HTML":


  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.