Compare sequence & maps headers in fasta file

Question

This is the perl code which compares the sequence in fasta file & maps the header. Though the code is working well, I still would like to make it more efficient. Since the files I compare has >100000 sequences it is taking huge time despite using hashes. Can you please suggest more efficient way ?

Example:

File1:
>Seq1
ABCDEFGH
>Seq2
RNADIMSEQ
>Seq6
XYZ

File2:
>Seq3
ABCDEFGH
>Seq4
RNADIMSEQ

Output:
>Seq1  >Seq3
>Seq2  >Seq4
>Seq6  Not found

The code: my $start_run = time();

%hash=();
open(out, ">Output.txt");
open(sbjt, "File1.fasta") or die "File not found"; #Bigger file
$count =0;
while(<sbjt>)
    {
    chomp;
    if($_ =~ m/^\w+/) 
        {
            $hash{$previous} = $_  ;
            #print "$previous\n";
        }
    else
        {
        $previous = $_;
        }
    $count++ ;  
    }
    close sbjt;
    #print "$hash{$previous}";

%dash=();
$previous = undef;
open(query, "File2.fasta") or die "File not found"; #smaller
while(<query>)
    {
    chomp;
    if($_ =~ m/^\w+/) 
        {
            $dash{$previous} = $_  ;
            #print "$previous\n";
        }
    else
        {
        $previous = $_;
        }
    }
    close query;

foreach $key (keys %dash)
{
         foreach $temp (keys %hash)
         {
            if($hash{$temp} eq $dash{$key})
            {
                push(@new_array, $temp);
            }
            next; #added
        }
        if(scalar @new_array>0)
        {       
            print out "$key\t@new_array\n";
        }
        else
        {
        print  out "$key\tNot found\n";
        }
        @new_array=();
 }

 my $end_run = time();
my $run_time = $end_run - $start_run;
print "Job took $run_time seconds\n";

Answer 1

Obvious things: use strict and warnings. Use three argument form of open, use lexical file handles, check the success of open or use autodie.

When I run the program, Seq6 is not reported in the output. That's because the outer loop exhausts File2, not File1.

#!/usr/bin/perl
use warnings;
use strict;

my $start_run = time;
my %hash;
my $previous;

open my $SBJ, '<', 'File1.fasta' or die 'File not found';
while (<$SBJ>) {
    chomp;
    if (/^\w+/) {
        $hash{$previous} = $_;
    } else {
        $previous = $_;
    }
}
close $SBJ;

my %dash;
undef $previous;
open my $QUERY, '<', 'File2.fasta' or die 'File not found';
while (<$QUERY>) {
    chomp;
    if (/^\w+/) {
        $dash{$previous} = $_;
    } else {
        $previous = $_;
    }
}
close $QUERY;

open my $OUT, '>', 'Output.txt' or die $!;
my @new_array;
for my $temp (keys %hash) {
    for my $key (keys %dash) {
        if($hash{$temp} eq $dash{$key}) {
            push @new_array, $key;
        }
    }
    if (@new_array > 0) {
        print {$OUT} "$temp\t@new_array\n";
    } else {
        print {$OUT} "$temp\tNot found\n";
    }
    @new_array = ();
}
close $OUT or die $!;

my $end_run = time;
my $run_time = $end_run - $start_run;
print "Job took $run_time seconds\n";

If performance is a problem, you can try trading memory for speed: Hash by the sequences, not by their names:

#!/usr/bin/perl
use warnings;
use strict;

my $start_run = time;
my %hash;

open my $SBJ, '<', 'File1.fasta' or die 'File not found';
my $id;
while (<$SBJ>) {
    chomp;
    if (/^>/) {
        $id = $_;
    } else {
        $hash{$_} = $id;
    }
}
close $SBJ;

my %dash;
undef $id;
open my $QUERY, '<', 'File2.fasta' or die 'File not found';
while (<$QUERY>) {
    chomp;
    if (/^>/) {
        $id = $_;
    } else {
        $dash{$_} = $id;
    }
}
close $QUERY;

open my $OUT, '>', 'Output.txt' or die $!;
for my $seq1 (keys %hash) {
    if (exists $dash{$seq1}) {
        print {$OUT} "$hash{$seq1}\t$dash{$seq1}\n";
    } else {
        print {$OUT} "$hash{$seq1}\tNot found\n";
    }
}
close $OUT or die $!;

my $end_run = time;
my $run_time = $end_run - $start_run;
print "Job took $run_time seconds\n";

Reading the files is exactly the same code with just the hash name different. I'd refactor it into a subroutine to keep the code DRY (Don't Repeat Yourself):

sub read_fasta {
    my $filename = shift;
    my ($id, %hash);
    open my $FH, '<', $filename or die "Can't open $filename: $!";
    while (<$FH>) {
        chomp;
        if (/^>/) {
            $id = $_;
        } else {
            $hash{$_} = $id;
        }
    }
    return %hash
}


my %hash = read_fasta('File1.fasta');
my %dash = read_fasta('File2.fasta');