http://qs321.pair.com?node_id=757871

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

dear all, i am having a weird problem. i am comparing a specific column between two large files. file1 is a database text file. and file2 is a query text file. i want to compare a particular column (which has alphabets) from file2 to file1. final results should have all the lines in file2, with an extra added column from file1. for example.
file1.db: zm811463 1190 31050 - A/G/T/C zm811462 1190 31051 - C/T zm1829427 1190 31789 + A/G zm445312 1190 31883 - A/G zm5377419 1190 32419 + A/C zm1052506 1190 32829 + C/G zm1052507 1190 32886 + A/C/T zm9115338 1190 33832 + A/G/CC file2.query 1190 31277 A > T 1 0 0 1190 31607 C > A 0 3 1 1190 31629 C > T 0 2 0 1190 31789 A > G 1 2 5 1190 31882 A > C 0 4 0 1190 31883 T > A 0 4 0 1190 31883 T > C 2 2 5 1190 32199 C > T 0 1 1 1190 32487 T > C 0 1 1 1190 32496 A > G 0 3 0 output which i am getting now: 1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 - Total number of HITS: 1 BUT, i want my output look like: 1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 31883 T > A 0 4 0 - 1190 31883 T > C 2 2 5 - 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 - Total number of HITS: 1
if you notice keenly, these two lines are missing in my output:
1190 31883 T > A 0 4 0 - 1190 31883 T > C 2 2 5 -
why is this happening? the program so far looks like this:
use strict; use warnings; use DB_File; my $myhashfile = "hash.$$"; tie my %hash1, "DB_File", $myhashfile, O_RDWR|O_CREAT, 0666, $DB_HASH or die "cannot open file $myhashfile: $!"; open(OUT,">output.out")or die "can not open"; open(my $fh1, "<", "file1.db") or die "file1.db: $!"; foreach (<$fh1>){ chomp; my @dbinfo = split(/\s+/); $hash1{"$dbinfo[1]#$dbinfo[2]"} = "$dbinfo[4]##$dbinfo[0]"; } close($fh1); my $c=0; open(my $fh2, "<", "file2.query") or die "file2.query: $!"; my @snplist; foreach (<$fh2>) { chomp($_);@snplist=(); my @queryinfo = split(/[\s>]+/); my $values = $hash1{"$queryinfo[0]#$queryinfo[1]"}; my ($variant,$rs)=split("##",$values); my $flag_left=0;my $flag_right=0; @snplist=split("/",$variant); if(defined($variant)){ foreach my $lis(@snplist){ if($queryinfo[2] eq $lis){$flag_left=1;} elsif($queryinfo[3] eq $lis){$flag_right=1;} } if(($flag_left == 1) && ($flag_right == 1)){ print OUT "$_\t$rs\n"; $c++; $flag_left=0;$flag_right=0; } } else{ print OUT "$_\t-\n"; } } print OUT "\nTotal number of HITS: $c\n"; close($fh2); untie %hash1; unlink($myhashfile); $myhashfile=();
using hash tie, because am dealing with large files.using @snplist because, there can be more alphabets to compare and only if both(seperated by ">") the alphabets in file2 third column is present in file1 5th column(which has many seperated by "/"), its considered as hit. plz help.thank you very much.

Replies are listed 'Best First'.
Re: duplicates getting omitted while comparing values inside foreach.
by ELISHEVA (Prior) on Apr 16, 2009 at 07:48 UTC
    Btree (or lack thereof) isn't the reason for this bug. The bug is in the placement of your print OUT statements. If you clean up the indenting you will notice that lines only print out if

    • ! defined($variant)
    • defined($variant) && ($flag_left==1) && ($flag_right==1)

    On the rows where it fails to print out, I presume that $variant is defined, but either the left or right flag is something other than 1.

    To fix this I recommend that you move the code for calculating the final column into a separate sub that returns either the file name or "-". Then your loop to read in file 2 should (pseudo code) look something like this:

    while (my $line = <DATA>) { my @aFields=split(/\s+/,$line); my $kDb= $aFields[0] . '#' . $aFields[1]; my ($variant,$rs) = split('##', $hash1{$kDb}); #code to calculate my $result=calcFinalColumnValue($variant,$rs,\@aFields); #print row print "@aFields $result\n"; }

    Also, you might want to consider two changes to make this code more efficient.

    • use while (<$fh2>) rather than foreach. I believe (someone correct me if I'm wrong), that foreach will cause the entire file to be slurped into an array before the loop starts. This could be a problem in a very large file and defeats the purpose of using a tied hash.
    • for the values of %hash1 use an array reference [$variant, $rs] rather than the string "foo##baz". This will eliminate the need to concatenate hash values when reading in file1 and the need to parse them when reading in file2.

    Best, beth

      i had missed an else statement after if.
      if(defined($variant)){ foreach my $lis(@snplist){ if($queryinfo[2] eq $lis){$flag_left=1;} elsif($queryinfo[3] eq $lis){$flag_right=1;} } if(($flag_left == 1) && ($flag_right == 1)){ print OUT "$_\t$rs\n"; $c++; $flag_left=0;$flag_right=0; } else{ print OUT "$_\t-\n"; } } else{ print OUT "$_\t-\n"; }
      thank u :)
      thank you so much for your valuable suggestions :)
Re: duplicates getting omitted while comparing values inside foreach.
by ambrus (Abbot) on Apr 16, 2009 at 08:33 UTC

    There's no real need to write such a complicated script for this though, simple use of unix command line tools are enough.

    Firstly, verify that the files are indeed sorted on the third and second column resp with these commands.

    sort -ck3 file1.db sort -ck2 file2.query

    We get no error, so they are sorted. (Note that we're doing a textual comparision, not a numerical one, because join can't handle the latter.)

    Then let's use the join command to join the two files.

    join -a 2 -1 3 -2 2 -e - -o 2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,1.1 file1. +db file2.query | tee result

    Admittedly the command-line looks complicated because we use just about any option join has, but we're lucky because we don't need any extra features join doesn't have. (If you have numbers in the join column that aren't of the same width, then you're out of luck here and likely need at least a simple perl script.)

    The output is this.

    1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 31883 T > A 0 4 0 zm445312 1190 31883 T > C 2 2 5 zm445312 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 -

    The whitespaces separating the fields are destroyed, but if that bothers you you can fix it up easily with this command.

    (sed 's/^\([^ ]*\) \([^ ]*\) \([^ ]* [^ ]* [^ ]*\) \([^ ]*\) \([^ ]*\) + \([^ ]*\)\(.*\)/\1 \2 \3 \4 \5 \6 \7/' result; ech +o; echo -n "Total number of HITS: "; grep -cv ' -$' result) |tee resu +lt.fmt

    The result is this:

    1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 31883 T > A 0 4 0 zm445312 1190 31883 T > C 2 2 5 zm445312 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 - Total number of HITS: 3
    Update 2009 sep 2.

    See Re^2: Joining two files on common field for a list of other nodes where unix textutils is suggested to merge files.

Re: duplicates getting omitted while comparing values inside foreach.
by ambrus (Abbot) on Apr 16, 2009 at 08:51 UTC

    And to prove that what the join utility does (see solution above) is not trivial but also not very complicated, let's try to emulate what it does in a short perl program.

    perl -we 'open $_, "<", shift or die for my($D, $Q); my $p = -1e9999; +my $s = -1e9999; my $x; while (my $q = <$Q>) { chomp($q); my $k = (sp +lit " ", $q)[1]; $p <= $k or die "query not sorted"; $p = $k; while ( +$s < $p and defined(my $d = <$D>)) { (my($l), $x) = (split " ", $d)[2 +,0]; $s <= $l or die "db not sorted"; $s = $l; } my $z = $p == $s ? $ +x : "-"; print $q, " ", $z, "\n"; }' file1.db file2.query

    The output is this.

    1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 31883 T > A 0 4 0 zm445312 1190 31883 T > C 2 2 5 zm445312 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 -

    It is an exercise to the reader to add the "Number of HITS" message.

      this is correct... so instead of 3 hits, it should be 1. final result should look like :
      1190 31277 A > T 1 0 0 - 1190 31607 C > A 0 3 1 - 1190 31629 C > T 0 2 0 - 1190 31789 A > G 1 2 5 zm1829427 1190 31882 A > C 0 4 0 - 1190 31883 T > A 0 4 0 - 1190 31883 T > C 2 2 5 - 1190 32199 C > T 0 1 1 - 1190 32487 T > C 0 1 1 - 1190 32496 A > G 0 3 0 - Total number of Hits: 1
      sorry for the confusion in between :(

        Please explain why this should be the answer. What's your criterion for joining actually?

      sorry for bothering you all, i made a silly mistake and i rectified it :) thanks for all your suggestions. now that there are other methods also for me to look at from all your suggestions :)
      the corrections i did was: if(defined($variant)){ foreach my $lis(@snplist){ if($queryinfo[2] eq $lis){$flag_left=1;} elsif($queryinfo[3] eq $lis){$flag_right=1;} } if(($flag_left == 1) && ($flag_right == 1)){ print OUT "$_\t$rs\n"; $c++; $flag_left=0;$flag_right=0; } else{ print OUT "$_\t-\n"; } } else{ print OUT "$_\t-\n"; }
      i should just add an else statement under if.. thats all :)
      sorry to say this, but
      1190 31883 T > A 0 4 0 zm445312 1190 31883 T > C 2 2 5 zm445312
      this code "zm445312" doesnt come under hit category. only "zm1829427" has to be reported.
      the number of hits i should get is 1, not 3.
Re: duplicates getting omitted while comparing values inside foreach.
by Anonymous Monk on Apr 16, 2009 at 07:13 UTC
      tried with DB_BTREE now, but still getting the same results(wrong one). :( :( :(

        I'm not sure what your code is supposed to do. But if you want to check whether one list is a subset of the other, possibly with duplicates, I recommend using a hash of array references and "colliding" the elements out:

        my %elements; # open left file while (<$left>) { chomp; # extract key my $key = $_; $elements{ $key } ||= []; push @{$elements{ $key }}, $_; }; # open right file while (<$right>) { chomp; # extract key my $key = $_; if (! exists $elements{ $key }) { print "Missing totally from left file: $_\n"; } else { my $match = shift @{$elements{ $key }}; print "Matched: $_ against $match\n"; if (0 == @{$elements{ $key }}) { delete $elements{ $key }; }; }; }; for (keys %elements) { print "Only in left file: @{ $elements{ $_ }}" };