code_perl/Maasha/Match.pm

   1 package Maasha::Match;
   2
   3 # Copyright (C) 2007 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines to match sequences
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use strict;
  32 use Data::Dumper;
  33 use Storable qw( dclone );
  34 use Maasha::Common;
  35 use Maasha::Fasta;
  36 use Maasha::Seq;
  37 use vars qw ( @ISA @EXPORT );
  38
  39 use constant {
  40     SEQ_NAME => 0,
  41     SEQ      => 1,
  42 };
  43
  44 @ISA = qw( Exporter );
  45
  46
  47 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  48
  49
  50 sub match_mummer
  51 {
  52     # Martin A. Hansen, June 2007.
  53
  54     # Match sequences using MUMmer.
  55
  56     my ( $entries1,   # FASTA entries
  57          $entries2,   # FASTA entries
  58          $options,    # additional MUMmer options - OPTIONAL
  59          $tmp_dir,    # temporary directory
  60        ) = @_;
  61
  62     # Returns a list.
  63
  64     my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );
  65
  66     $tmp_dir ||= $ENV{ "BP_TMP" };
  67
  68     $options->{ "word_size" } ||= 20;
  69     $options->{ "direction" } ||= "both";
  70
  71     push @args, "-c";
  72     push @args, "-L";
  73     push @args, "-F";
  74     push @args, "-l $options->{ 'word_size' }";
  75     push @args, "-maxmatch";
  76     push @args, "-n" if not Maasha::Seq::seq_guess_type( $entries1->[ 0 ]->[ 1 ] ) eq "protein";
  77     push @args, "-b" if $options->{ "direction" } =~ /^b/;
  78     push @args, "-r" if $options->{ "direction" } =~ /^r/;
  79
  80     $arg = join " ", @args;
  81
  82     $file_in1 = "$tmp_dir/muscle1.tmp";
  83     $file_in2 = "$tmp_dir/muscle2.tmp";
  84     $file_out = "$tmp_dir/muscle3.tmp";
  85
  86     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
  87     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };
  88
  89     Maasha::Fasta::put_entries( $entries1, $file_in1 );
  90     Maasha::Fasta::put_entries( $entries2, $file_in2 );
  91
  92     Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );
  93
  94     $fh = Maasha::Common::read_open( $file_out );
  95
  96     while ( $line = <$fh> )
  97     {
  98         chomp $line;
  99
 100         if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
 101         {
 102             $result->{ "Q_ID" }  = $1;
 103             $result->{ "Q_LEN" } = $2;
 104             $result->{ "DIR" }   = "reverse";
 105         }
 106         elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
 107         {
 108             $result->{ "Q_ID" }  = $1;
 109             $result->{ "Q_LEN" } = $2;
 110             $result->{ "DIR" }   = "forward";
 111         }
 112         elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
 113         {
 114             $result->{ "S_ID" }    = $1;
 115             $result->{ "S_BEG" }   = $2 - 1;
 116             $result->{ "Q_BEG" }   = $3 - 1;
 117             $result->{ "HIT_LEN" } = $4;
 118             $result->{ "S_END" }   = $result->{ "S_BEG" } + $result->{ "HIT_LEN" } - 1;
 119             $result->{ "Q_END" }   = $result->{ "Q_BEG" } + $result->{ "HIT_LEN" } - 1;
 120
 121             push @results, dclone $result;
 122         }
 123
 124     }
 125
 126     unlink $file_in1;
 127     unlink $file_in2;
 128     unlink $file_out;
 129
 130     return wantarray ? @results : \@results;
 131 }
 132
 133
 134 sub match_vmatch
 135 {
 136     # Martin A. Hansen, April 2008.
 137
 138     # Vmatches a list of records against a list of index files and the full
 139     # path to the result file is returned.
 140
 141     my ( $tmp_dir,       # directory in where to save temp files
 142          $records,       # list of records
 143          $index_files,   # list of index files
 144          $options,       # argument hash
 145        ) = @_;
 146
 147     # Returns a string.
 148
 149     my ( $query_file, $result_file, @result_files, $fh_in, $fh_out, $line, @fields, $i, $record, $vmatch_args, @index_names, @seq_names, $count_list );
 150
 151     $query_file  = "$tmp_dir/query.seq";
 152     $result_file = "$tmp_dir/vmatch.out";
 153
 154     $fh_out = Maasha::Common::write_open( $query_file );
 155
 156     foreach $record ( @{ $records } )
 157     {
 158         if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
 159         {
 160             next if length $record->{ "SEQ" } < 12; # assuming that the index is created for 12 as minimum length
 161
 162             push @seq_names, $record->{ "SEQ_NAME" };
 163
 164             Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_out, 80 );
 165         }
 166     }
 167
 168     close $fh_out;
 169
 170     if ( $options->{ 'genome' } ) {
 171         $vmatch_args  = "-complete -d -p -q $query_file";
 172     } else {
 173         $vmatch_args  = "-complete -d -p -showdesc 100 -q $query_file";
 174     }
 175
 176     $vmatch_args .= " -h " . $options->{ "hamming_dist" } if $options->{ "hamming_dist" };
 177     $vmatch_args .= " -e " . $options->{ "edit_dist" }    if $options->{ "edit_dist" };
 178
 179     for ( $i = 0; $i < @{ $index_files }; $i++ )
 180     {
 181         Maasha::Common::run( "vmatch", "$vmatch_args $index_files->[ $i ] > $result_file.$i" );
 182
 183         push @result_files, "$result_file.$i";
 184     }
 185
 186     unlink $query_file;
 187
 188     $count_list = vmatch_count_hits( \@result_files ) if ( $options->{ "count" } );
 189
 190     $fh_out = Maasha::Common::write_open( $result_file );
 191
 192     for ( $i = 0; $i < @{ $index_files }; $i++ )
 193     {
 194         $index_files->[ $i ] =~ s/.+\/(.+)\.fna$/$1/ if $options->{ 'genome' };
 195
 196         $fh_in = Maasha::Common::read_open( "$result_file.$i" );
 197
 198         while ( $line = <$fh_in> )
 199         {
 200             chomp $line;
 201
 202             next if $line =~ /^#/;
 203
 204             @fields = split " ", $line;
 205
 206             next if $options->{ "max_hits" } and $count_list->[ $fields[ 5 ] ] > $options->{ 'max_hits' };
 207
 208             $fields[ 1 ] = $index_files->[ $i ];                                     # S_ID
 209             $fields[ 9 ] = $count_list->[ $fields[ 5 ] ] if $options->{ "count" };   # SCORE
 210             $fields[ 5 ] = $seq_names[ $fields[ 5 ] ];                               # Q_ID
 211
 212             print $fh_out join( "\t", @fields ), "\n";
 213         }
 214
 215         close $fh_in;
 216
 217         unlink "$result_file.$i";
 218     }
 219
 220     close $fh_out;
 221
 222     return $result_file;
 223 }
 224
 225
 226 sub vmatch_count_hits
 227 {
 228     # Martin A. Hansen, April 2008.
 229
 230     # Given a list of Vmatch result file, count duplications based
 231     # on q_id. The counts are returned in a list where the list index
 232     # corresponds to the q_id index in the query file.
 233
 234     my ( $files,   # vmatch result files
 235        ) = @_;
 236
 237     # Returns a list.
 238
 239     my ( $file, $fh_in, $line, @fields, @count_list );
 240
 241     foreach $file ( @{ $files } )
 242     {
 243         $fh_in = Maasha::Common::read_open( $file );
 244
 245         while ( $line = <$fh_in> )
 246         {
 247             chomp $line;
 248
 249             next if $line =~ /^#/;
 250
 251             @fields = split " ", $line;
 252
 253             $count_list[ $fields[ 5 ] ]++;
 254         }
 255
 256         close $fh_in;
 257     }
 258
 259     return wantarray ? @count_list : \@count_list;
 260 }
 261
 262
 263 sub vmatch_get_entry
 264 {
 265     # Martin A. Hansen, January 2008.
 266
 267     # Parses vmatch output records.
 268
 269     my ( $fh,   # file handle to vmatch result file.
 270        ) = @_;
 271
 272     # Returns a hash.
 273
 274     my ( $line, @fields, %record );
 275
 276     while ( $line = <$fh> )
 277     {
 278         chomp $line;
 279
 280         next if $line =~ /^#/;
 281
 282         @fields = split "\t", $line;
 283
 284         $record{ "REC_TYPE" } = "VMATCH";
 285
 286         $record{ "S_LEN" }      = $fields[ 0 ];
 287         $record{ "S_ID" }       = $fields[ 1 ];
 288         $record{ "S_BEG" }      = $fields[ 2 ];
 289
 290         if ( $fields[ 3 ] eq "D" ) {
 291             $record{ "STRAND" } = "+";
 292         } else {
 293             $record{ "STRAND" } = "-";
 294         }
 295
 296         $record{ "Q_LEN" }      = $fields[ 4 ];
 297         $record{ "Q_ID" }       = $fields[ 5 ];
 298         $record{ "Q_BEG" }      = $fields[ 6 ];
 299         $record{ "MATCH_DIST" } = $fields[ 7 ];
 300         $record{ "E_VAL" }      = $fields[ 8 ];
 301         $record{ "SCORE" }      = $fields[ 9 ];
 302         $record{ "IDENT" }      = $fields[ 10 ];
 303
 304         $record{ "Q_END" }      = $record{ "Q_BEG" } + $record{ "Q_LEN" } - 1;
 305         $record{ "S_END" }      = $record{ "S_BEG" } + $record{ "S_LEN" } - 1;
 306
 307         return wantarray ? %record : \%record;
 308     }
 309 }
 310
 311
 312 sub vmatch_index
 313 {
 314     # Martin A. Hansen, July 2008.
 315
 316     # Use mkvtree to create a vmatch index of a given file
 317
 318     my ( $file,      # FASTA file to index
 319          $src_dir,   # source directory with FASTA file
 320          $dst_dir,   # distination directory for Vmatch index
 321          $tmp_dir,   # temp directory - OPTIONAL
 322        ) = @_;
 323
 324     # Returns nothing.
 325
 326     my ( $fh, $entry, $tmp_file );
 327
 328     Maasha::Common::dir_create_if_not_exists( $dst_dir );
 329
 330     if ( Maasha::Common::file_size( $file ) < 200_000_000 )
 331     {
 332         &Maasha::Common::run( "mkvtree", "-db $src_dir/$file -dna -pl -allout -indexname $dst_dir/$file > /dev/null 3>&1" );
 333     }
 334     else
 335     {
 336         $fh = Maasha::Common::read_open( "$src_dir/$file" );
 337
 338         while ( $entry = Maasha::Fasta::get_entry( $fh ) )
 339         {
 340             $tmp_file = $entry->[ SEQ_NAME ] . ".fna";
 341
 342             Maasha::Fasta::put_entries( [ $entry ], "$tmp_dir/$tmp_file" )
 343
 344             &Maasha::Common::run( "mkvtree", "-db $tmp_dir/$tmp_file -dna -pl -allout -indexname $dst_dir/$tmp_file > /dev/null 3>&1" );
 345
 346             unlink "$tmp_dir/$tmp_file";
 347         }
 348
 349         close $fh;
 350     }
 351 }
 352
 353
 354 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 355
 356
 357 __END__