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 warnings;
  32 use strict;
  33 use Data::Dumper;
  34 use Storable qw( dclone );
  35 use Maasha::Common;
  36 use Maasha::Filesys;
  37 use Maasha::Fasta;
  38 use Maasha::Seq;
  39 use Maasha::Biopieces;
  40 use vars qw ( @ISA @EXPORT );
  41
  42 use constant {
  43     SEQ_NAME => 0,
  44     SEQ      => 1,
  45 };
  46
  47 @ISA = qw( Exporter );
  48
  49
  50 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  51
  52
  53 sub match_mummer
  54 {
  55     # Martin A. Hansen, June 2007.
  56
  57     # Match sequences using MUMmer.
  58
  59     my ( $entries1,   # FASTA entries
  60          $entries2,   # FASTA entries
  61          $options,    # additional MUMmer options - OPTIONAL
  62          $tmp_dir,    # temporary directory
  63        ) = @_;
  64
  65     # Returns a list.
  66
  67     my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );
  68
  69     $tmp_dir ||= $ENV{ "BP_TMP" };
  70
  71     $options->{ "word_size" } ||= 20;
  72     $options->{ "direction" } ||= "both";
  73
  74     push @args, "-c";
  75     push @args, "-L";
  76     push @args, "-F";
  77     push @args, "-l $options->{ 'word_size' }";
  78     push @args, "-maxmatch";
  79     push @args, "-n" if not Maasha::Seq::seq_guess_type( $entries1->[ 0 ]->[ 1 ] ) eq "PROTEIN";
  80     push @args, "-b" if $options->{ "direction" } =~ /^b/;
  81     push @args, "-r" if $options->{ "direction" } =~ /^r/;
  82
  83     $arg = join " ", @args;
  84
  85     $file_in1 = "$tmp_dir/muscle1.tmp";
  86     $file_in2 = "$tmp_dir/muscle2.tmp";
  87     $file_out = "$tmp_dir/muscle3.tmp";
  88
  89     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
  90     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };
  91
  92     Maasha::Fasta::put_entries( $entries1, $file_in1 );
  93     Maasha::Fasta::put_entries( $entries2, $file_in2 );
  94
  95     Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );
  96
  97     $fh = Maasha::Filesys::file_read_open( $file_out );
  98
  99     while ( $line = <$fh> )
 100     {
 101         chomp $line;
 102
 103         if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
 104         {
 105             $result->{ "Q_ID" }  = $1;
 106             $result->{ "Q_LEN" } = $2;
 107             $result->{ "DIR" }   = "reverse";
 108         }
 109         elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
 110         {
 111             $result->{ "Q_ID" }  = $1;
 112             $result->{ "Q_LEN" } = $2;
 113             $result->{ "DIR" }   = "forward";
 114         }
 115         elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
 116         {
 117             $result->{ "S_ID" }    = $1;
 118             $result->{ "S_BEG" }   = $2 - 1;
 119             $result->{ "Q_BEG" }   = $3 - 1;
 120             $result->{ "HIT_LEN" } = $4;
 121             $result->{ "S_END" }   = $result->{ "S_BEG" } + $result->{ "HIT_LEN" } - 1;
 122             $result->{ "Q_END" }   = $result->{ "Q_BEG" } + $result->{ "HIT_LEN" } - 1;
 123
 124             push @results, dclone $result;
 125         }
 126     }
 127
 128     unlink $file_in1;
 129     unlink $file_in2;
 130     unlink $file_out;
 131
 132     return wantarray ? @results : \@results;
 133 }
 134
 135
 136 sub vmatch_index
 137 {
 138     # Martin A. Hansen, July 2008.
 139
 140     # Use mkvtree to create a vmatch index of a given file
 141
 142     my ( $file,      # FASTA file to index
 143          $src_dir,   # source directory with FASTA file
 144          $dst_dir,   # distination directory for Vmatch index
 145          $tmp_dir,   # temp directory - OPTIONAL
 146        ) = @_;
 147
 148     # Returns nothing.
 149
 150     my ( $fh, $entry, $tmp_file );
 151
 152     Maasha::Filesys::dir_create_if_not_exists( $dst_dir );
 153
 154 #    if ( Maasha::Common::file_size( $file ) < 200_000_000 )
 155 #    {
 156 #        &Maasha::Common::run( "mkvtree", "-db $src_dir/$file -dna -pl -allout -indexname $dst_dir/$file > /dev/null 3>&1" );
 157 #    }
 158 #    else
 159 #    {
 160         $fh = Maasha::Filesys::file_read_open( "$src_dir/$file" );
 161
 162         while ( $entry = Maasha::Fasta::get_entry( $fh ) )
 163         {
 164             $tmp_file = $entry->[ SEQ_NAME ] . ".fna";
 165
 166             Maasha::Fasta::put_entries( [ $entry ], "$tmp_dir/$tmp_file" );
 167
 168             Maasha::Common::run( "mkvtree", "-db $tmp_dir/$tmp_file -dna -pl -allout -indexname $dst_dir/$tmp_file > /dev/null 3>&1" );
 169
 170             unlink "$tmp_dir/$tmp_file";
 171         }
 172
 173         close $fh;
 174 }
 175
 176
 177 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 178
 179
 180 __END__