]> git.donarmstrong.com Git - biopieces.git/blobdiff - code_perl/Maasha/Biopieces.pm
correcting bug in Match.pm
[biopieces.git] / code_perl / Maasha / Biopieces.pm
index ee38e282d4392a56f181713480cfe66b59d40576..9e2d81a2eff6665952b2529d3e5951f433de28d1 100644 (file)
@@ -236,6 +236,7 @@ sub run_script
     elsif ( $script eq "rename_keys" )              { script_rename_keys(               $in, $out, $options ) }
     elsif ( $script eq "uniq_vals" )                { script_uniq_vals(                 $in, $out, $options ) }
     elsif ( $script eq "merge_vals" )               { script_merge_vals(                $in, $out, $options ) }
+    elsif ( $script eq "merge_records" )            { script_merge_records(             $in, $out, $options ) }
     elsif ( $script eq "grab" )                     { script_grab(                      $in, $out, $options ) }
     elsif ( $script eq "compute" )                  { script_compute(                   $in, $out, $options ) }
     elsif ( $script eq "flip_tab" )                 { script_flip_tab(                  $in, $out, $options ) }
@@ -574,6 +575,7 @@ sub get_options
         @options = qw(
             in_file|i=s
             genome|g=s
+            seed_size|s=s
             mismatches|m=s
             gap_size|G=s
             cpus|c=s
@@ -748,6 +750,13 @@ sub get_options
             delimit|d=s
         );
     }
+    elsif ( $script eq "merge_records" )
+    {
+        @options = qw(
+            keys|k=s
+            merge|m=s
+        );
+    }
     elsif ( $script eq "grab" )
     {
         @options = qw(
@@ -1008,6 +1017,10 @@ sub get_options
         {
             Maasha::Common::error( qq(Character '$options{ $opt }' is not allowed in table names) );
         }
+        elsif ( $opt eq "merge" and $options{ $opt } !~ /^(AandB|AorB|BorA|AnotB|BnotA)$/ )
+        {
+            Maasha::Common::error( qq(Argument to --$opt must be AandB, AorB, BorA, AnotB, or BnotA - not "$options{ $opt }") );
+        }
     }
 
     Maasha::Common::error( qq(no --database specified) )                if $script eq "create_blast_db"     and not $options{ "database" };
@@ -3580,23 +3593,31 @@ sub script_soap_seq
 
     # Returns nothing.
 
-    my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry );
+    my ( $genome, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $line, @fields, $entry, $count, $args );
 
+    $options->{ "seed_size" }  ||= 10;
     $options->{ "mismatches" } ||= 2;
     $options->{ "gap_size" }   ||= 0;
     $options->{ "cpus" }       ||= 1;
 
-    $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna" if $options->{ 'genome' };
+    if ( $options->{ "genome" } ) {
+        $options->{ "in_file" } = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/fasta/$options->{ 'genome' }.fna";
+    }
 
     $tmp_in  = "$BP_TMP/soap_query.seq";
     $tmp_out = "$BP_TMP/soap.result";
 
     $fh_out = Maasha::Common::write_open( $tmp_in );
+    $count = 0;
 
     while ( $record = get_record( $in ) ) 
     {
-        if ( $entry = record2fasta( $record ) ) {
+        if ( $entry = record2fasta( $record ) )
+        {
             Maasha::Fasta::put_entry( $entry, $fh_out );
+
+            $count++;
         }
 
         put_record( $record, $out );
@@ -3604,32 +3625,48 @@ sub script_soap_seq
 
     close $fh_out;
 
-    Maasha::Common::run( "soap", "-r 2 -a $tmp_in -v $options->{ 'mismatches' } -g $options->{ 'gap_size' } -p $options->{ 'cpus' } -d $options->{ 'in_file' } -o $tmp_out > /dev/null 2>&1", 1 );
+    if ( $count > 0 )
+    {
+        $args = join( " ",
+            "-s $options->{ 'seed_size' }",
+            "-r 2",
+            "-a $tmp_in",
+            "-v $options->{ 'mismatches' }",
+            "-g $options->{ 'gap_size' }",
+            "-p $options->{ 'cpus' }",
+            "-d $options->{ 'in_file' }",
+            "-o $tmp_out",
+        );
 
-    unlink $tmp_in;
+        $args .= " > /dev/null 2>&1" if not $options->{ 'verbose' };
 
-    $fh_out = Maasha::Common::read_open( $tmp_out );
+        Maasha::Common::run( "soap", $args, 1 );
 
-    undef $record;
+        unlink $tmp_in;
 
-    while ( $line = <$fh_out> )
-    {
-        chomp $line;
+        $fh_out = Maasha::Common::read_open( $tmp_out );
 
-        @fields = split /\t/, $line;
+        undef $record;
 
-        $record->{ "REC_TYPE" }   = "SOAP";
-        $record->{ "Q_ID" }       = $fields[ 0 ];
-        $record->{ "SCORE" }      = $fields[ 3 ];
-        $record->{ "STRAND" }     = $fields[ 6 ];
-        $record->{ "S_ID" }       = $fields[ 7 ];
-        $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
-        $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
+        while ( $line = <$fh_out> )
+        {
+            chomp $line;
 
-        put_record( $record, $out );
-    }
+            @fields = split /\t/, $line;
 
-    close $fh_out;
+            $record->{ "REC_TYPE" }   = "SOAP";
+            $record->{ "Q_ID" }       = $fields[ 0 ];
+            $record->{ "SCORE" }      = $fields[ 3 ];
+            $record->{ "STRAND" }     = $fields[ 6 ];
+            $record->{ "S_ID" }       = $fields[ 7 ];
+            $record->{ "S_BEG" }      = $fields[ 8 ] - 1; # soap is 1-based
+            $record->{ "S_END" }      = $fields[ 8 ] + $fields[ 5 ] - 2;
+
+            put_record( $record, $out );
+        }
+
+        close $fh_out;
+    }
 
     unlink $tmp_out;
 }
@@ -4570,6 +4607,177 @@ sub script_merge_vals
 }
 
 
+sub script_merge_records
+{
+    # Martin A. Hansen, July 2008.
+
+    # Merges records in the stream based on identical values of two given keys.
+
+    my ( $in,        # handle to in stream
+         $out,       # handle to out stream
+         $options,   # options hash
+       ) = @_;
+
+    # Returns nothing.
+
+    my ( $merge, $record, $file1, $file2, $fh1, $fh2, $key1, $key2, @keys1, @keys2, @vals1, @vals2,
+         $num1, $num2, $num, $cmp, $i );
+
+    $merge = $options->{ "merge" } || "AandB";
+
+    $file1 = "$BP_TMP/merge_records1.tmp";
+    $file2 = "$BP_TMP/merge_records2.tmp";
+
+    $fh1   = Maasha::Common::write_open( $file1 );
+    $fh2   = Maasha::Common::write_open( $file2 );
+
+    $key1  = $options->{ "keys" }->[ 0 ];
+    $key2  = $options->{ "keys" }->[ 1 ];
+
+    $num   = $key2 =~ s/n$//;
+    $num1  = 0;
+    $num2  = 0;
+
+    while ( $record = get_record( $in ) ) 
+    {
+        if ( exists $record->{ $key1 } )
+        {
+            @keys1 = $key1;
+            @vals1 = $record->{ $key1 };
+
+            delete $record->{ $key1 };
+
+            map { push @keys1, $_; push @vals1, $record->{ $_ } } keys %{ $record };
+
+            print $fh1 join( "\t", @vals1 ), "\n";
+
+            $num1++;
+        }
+        elsif ( exists $record->{ $key2 } )
+        {
+            @keys2 = $key2;
+            @vals2 = $record->{ $key2 };
+
+            delete $record->{ $key2 };
+
+            map { push @keys2, $_; push @vals2, $record->{ $_ } } keys %{ $record };
+
+            print $fh2 join( "\t", @vals2 ), "\n";
+
+            $num2++;
+        }
+    }
+
+    close $fh1;
+    close $fh2;
+
+    if ( $num )
+    {
+        Maasha::Common::run( "sort", "-k 1,1n $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
+        Maasha::Common::run( "sort", "-k 1,1n $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
+    }
+    else
+    {
+        Maasha::Common::run( "sort", "-k 1,1 $file1 > $file1.sort" ) and rename "$file1.sort", $file1;
+        Maasha::Common::run( "sort", "-k 1,1 $file2 > $file2.sort" ) and rename "$file2.sort", $file2;
+    }
+
+    $fh1 = Maasha::Common::read_open( $file1 );
+    $fh2 = Maasha::Common::read_open( $file2 );
+
+    @vals1 = Maasha::Common::get_fields( $fh1 );
+    @vals2 = Maasha::Common::get_fields( $fh2 );
+
+    while ( $num1 > 0 and $num2 > 0 )
+    {
+        undef $record;
+
+        if ( $num ) {
+            $cmp = $vals1[ 0 ] <=> $vals2[ 0 ];
+        } else {
+            $cmp = $vals1[ 0 ] cmp $vals2[ 0 ];
+        }
+
+        if ( $cmp < 0 )
+        {
+            if ( $merge =~ /^(AorB|AnotB)$/ )
+            {
+                for ( $i = 0; $i < @keys1; $i++ ) {
+                    $record->{ $keys1[ $i ] } = $vals1[ $i ];
+                }
+
+                put_record( $record, $out );
+            }
+
+            @vals1 = Maasha::Common::get_fields( $fh1 );
+            $num1--;
+        }
+        elsif ( $cmp > 0 )
+        {
+            if ( $merge =~ /^(BorA|BnotA)$/ )
+            {
+                for ( $i = 0; $i < @keys2; $i++ ) {
+                    $record->{ $keys2[ $i ] } = $vals2[ $i ];
+                }
+
+                put_record( $record, $out );
+            }
+
+            @vals2 = Maasha::Common::get_fields( $fh2 );
+            $num2--;
+        }
+        else
+        {
+            if ( $merge =~ /^(AandB|AorB|BorA)$/ )
+            {
+                for ( $i = 0; $i < @keys1; $i++ ) {
+                    $record->{ $keys1[ $i ] } = $vals1[ $i ];
+                }
+
+                for ( $i = 1; $i < @keys2; $i++ ) {
+                    $record->{ $keys2[ $i ] } = $vals2[ $i ];
+                }
+            
+                put_record( $record, $out );
+            }
+
+            @vals1 = Maasha::Common::get_fields( $fh1 );
+            @vals2 = Maasha::Common::get_fields( $fh2 );
+            $num1--;
+            $num2--;
+        }
+    }
+
+    close $fh1;
+    close $fh2;
+
+    unlink $file1;
+    unlink $file2;
+
+    if ( $num1 > 0 and $merge =~ /^(AorB|AnotB)$/ )
+    {
+        undef $record;
+
+        for ( $i = 0; $i < @keys1; $i++ ) {
+            $record->{ $keys1[ $i ] } = $vals1[ $i ];
+        }
+
+        put_record( $record, $out );
+    }
+
+    if ( $num2 > 0 and $merge =~ /^(BorA|BnotA)$/ )
+    {
+        undef $record;
+
+        for ( $i = 0; $i < @keys2; $i++ ) {
+            $record->{ $keys2[ $i ] } = $vals2[ $i ];
+        }
+
+        put_record( $record, $out );
+    }
+}
+
+
 sub script_grab
 {
     # Martin A. Hansen, August 2007.
@@ -5065,6 +5273,7 @@ sub script_count_vals
 
     $fh_out   = Maasha::Common::write_open( $tmp_file );
 
+    $cache    = 0;
     $num      = 0;
 
     while ( $record = get_record( $in ) ) 
@@ -5091,7 +5300,7 @@ sub script_count_vals
 
     if ( $cache )
     {
-        $num      = 0;
+        $num   = 0;
 
         $fh_in = Maasha::Common::read_open( $tmp_file );
 
@@ -6079,21 +6288,24 @@ sub clean_tmp
 
     foreach $dir ( @dirs )
     {
-        if ( $dir =~ /(.+)_(\d+)_(\d+)_bp_tmp/ )
+        if ( $dir =~ /^$tmpdir\/(.+)_(\d+)_(\d+)_bp_tmp$/ )
         {
             $user = $1;
             $sid  = $2;
             $pid  = $3;
 
-            if ( not Maasha::Common::process_running( $pid ) )
-            {
-                # print STDERR "Removing stale dir: $dir\n";
-                Maasha::Common::dir_remove( $dir );
-            }
-            elsif ( $pid == $curr_pid )
+            if ( $user eq Maasha::Common::get_user() )
             {
-                # print STDERR "Removing current dir: $dir\n";
-                Maasha::Common::dir_remove( $dir );
+                if ( not Maasha::Common::process_running( $pid ) )
+                {
+                    # print STDERR "Removing stale dir: $dir\n";
+                    Maasha::Common::dir_remove( $dir );
+                }
+                elsif ( $pid == $curr_pid )
+                {
+                    # print STDERR "Removing current dir: $dir\n";
+                    Maasha::Common::dir_remove( $dir );
+                }
             }
         }
     }