Reading file into a hash

PerlSufi has asked for the wisdom of the Perl Monks concerning the following question:

UPDATE: Solved and removed if(defined $seq) line
Hello monks, I have been trying a bio perl challenge site- so meta hints are welcome, too. I am having trouble reading a file into a hash properly. My data will be a fasta file in the format of:

>sequence_5849
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>sequence_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>sequence_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
[download]

my code is:

#! /usr/bin/perl

use strict;
use warnings;
use Data::Dumper;
use feature qw(say);


my $file = 'file.txt';
open (my $fh, '<', $file) or die "Could not open file '$file' $!";

my (%sequence_hash, $header, $seq, $count);

while ( my $line = <$fh> ) {
    chomp($line);
    if ( $line =~ m/^>(.*)/ ) {
        if ( $seq ) {
            say $seq;
            $sequence_hash{$header} = $seq;
        }
        $header = $1;
        $seq = '';
    }
    else {
        $seq .= $line;
    }
}

close $fh;

print Dumper(\%sequence_hash);
[download]

My problem is that I am not getting the header and sequence. I have a feeling that its because of clearing out the $seq variable, but I am not sure how else to get the header and sequences in the hash. Any insight would be highly appreciated :)

Comment on Reading file into a hash Select or Download Code

Replies are listed 'Best First'.
Re: Reading file into a hash by davido (Cardinal) on May 28, 2014 at 16:56 UTC
Your records always start with a '`>`'. But you don't have to look at that character as the start of a record, you can look at it as a record separator. Then if you set `$/ = '>';`, your file reads will not read lines, but records. That (in my mind) simplifies the process. Here's one way to do it: use strict; use warnings; use Data::Dumper; local $/ = '>'; my %hash; while( <DATA> ) { if( my($k,$v) = / ^(\N+)\n # Capture the key. ([CTGA]+)$ # Capture the value. /mx ) { $hash{$k}=$v; } } print Dumper \%hash; __DATA__ >sequence_5849 CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAAT +AATTCTGAGG >sequence_5959 CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTT +GTCAGCAGACACGC >sequence_0808 CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTG +CGGGCAGTAGGTGGAAT [download] Update: My solution assumes that there aren't newlines embedded in the CTG.... sequences. If there are, your regex should probably accept them, and then use `tr///` to filter them out again. Also, since you're modifying the input record separator, it's wise to constrain its effects to as narrow a scope as possible. The following does that: `print Dumper sub { local $/ = '>'; my( $fh, %recs ) = shift; while( <$fh> ) { if( my($k,$v) = / ^(\N+)\n # Capture the key. ([CTGA]+)$ # Capture the value. /mx ) { $recs{$k}=$v; } } return \%recs; }->(\*DATA);` [download] Dave	[reply] [d/l] [select]
Re^2: Reading file into a hash by PerlSufi (Friar) on May 29, 2014 at 19:45 UTC
Very Nice, davido. I have not used $/ much so I will keep that in mind for future reference :)	[reply]
Re: Reading file into a hash by hippo (Archbishop) on May 28, 2014 at 14:51 UTC
I think that the problem is that you are setting the value of $seq to be '' but then testing it via `defined`, which (other than the first time through) it will be. Change your test from `if ( defined $seq )` to just `if ($seq)` and see what difference that makes.	[reply] [d/l] [select]
Re^2: Reading file into a hash by PerlSufi (Friar) on May 28, 2014 at 14:53 UTC
Hi hippo, Thanks for the response. I changed it and still did not get the last header and sequence	[reply]
Re^3: Reading file into a hash (last record problem) by LanX (Saint) on May 28, 2014 at 15:14 UTC
Thats a very common problem, you are trying to add a record only when the successor is to be parsed, but the last record has no successor (sic ;)! Most people try to solve by repeating code to add the last record after the loop. But it's much cleaner this way (avoiding a posteriori state logic) use strict; use warnings; use Data::Dump; my $header; my %sequence; while ( my $line = <DATA> ){ chomp $line; if ( $line =~ /^>(.)/ ) { $header = $1; } else { $sequence{$header} .= $line; } } dd \%sequence; __DATA__ >sequence_5849 CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC TCCCACTAATAATTCTGAGG >sequence_5959 CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT ATATCCATTTGTCAGCAGACACGC >sequence_0808 CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC TGGGAACCTGCGGGCAGTAGGTGGAAT [download] output `{ sequence_0808 => "CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGAC +CAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT", sequence_5849 => "CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTC +CGGCCTTCCCTCCCACTAATAATTCTGAGG", sequence_5959 => "CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCG +CCGAAGGTCTATATCCATTTGTCAGCAGACACGC", }` [download] Cheers Rolf ( addicted to the Perl Programming Language) update a general pattern* to solve such problems while staying DRY is to use references `if ( $line =~ /^(HEAD_PATTERN)/ ) { $data = \ $deeply{nested}{structure}{$1}; # reference data + } else { $$data .= $line; # derefrence data }` [download] like this you don't need to repeat the path of a deeply nested data structure, which might vary in multiple dimensions update added some explanation	[reply] [d/l] [select]
Re: Reading file into a hash by LanX (Saint) on May 28, 2014 at 14:53 UTC
at first glanceš this doesn't seem to make much sense: `if ( defined $seq ) { # ... $seq = '';` [download] empty strings are defined, you may either want to assign `undef` or just use that '' is false. Cheers Rolf ( addicted to the Perl Programming Language) š) that is w/o trying to understand the full intended logic	[reply] [d/l] [select]
Re^2: Reading file into a hash by PerlSufi (Friar) on May 28, 2014 at 14:56 UTC
Hi Lanx, I tried that, too. I also updated my post to reflect not using `if ( defined $seq)` anymore	[reply] [d/l]

update

update