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