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::Filesys;
  36 use Maasha::Fasta;
  37 use Maasha::Seq;
  38 use Maasha::Biopieces;
  39 use vars qw ( @ISA @EXPORT );
  40
  41 use constant {
  42     SEQ_NAME => 0,
  43     SEQ      => 1,
  44 };
  45
  46 @ISA = qw( Exporter );
  47
  48
  49 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  50
  51
  52 sub match_mummer
  53 {
  54     # Martin A. Hansen, June 2007.
  55
  56     # Match sequences using MUMmer.
  57
  58     my ( $entries1,   # FASTA entries
  59          $entries2,   # FASTA entries
  60          $options,    # additional MUMmer options - OPTIONAL
  61          $tmp_dir,    # temporary directory
  62        ) = @_;
  63
  64     # Returns a list.
  65
  66     my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );
  67
  68     $tmp_dir ||= $ENV{ "BP_TMP" };
  69
  70     $options->{ "word_size" } ||= 20;
  71     $options->{ "direction" } ||= "both";
  72
  73     push @args, "-c";
  74     push @args, "-L";
  75     push @args, "-F";
  76     push @args, "-l $options->{ 'word_size' }";
  77     push @args, "-maxmatch";
  78     push @args, "-n" if not Maasha::Seq::seq_guess_type( $entries1->[ 0 ]->[ 1 ] ) eq "protein";
  79     push @args, "-b" if $options->{ "direction" } =~ /^b/;
  80     push @args, "-r" if $options->{ "direction" } =~ /^r/;
  81
  82     $arg = join " ", @args;
  83
  84     $file_in1 = "$tmp_dir/muscle1.tmp";
  85     $file_in2 = "$tmp_dir/muscle2.tmp";
  86     $file_out = "$tmp_dir/muscle3.tmp";
  87
  88     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
  89     map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };
  90
  91     Maasha::Fasta::put_entries( $entries1, $file_in1 );
  92     Maasha::Fasta::put_entries( $entries2, $file_in2 );
  93
  94     Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );
  95
  96     $fh = Maasha::Filesys::file_read_open( $file_out );
  97
  98     while ( $line = <$fh> )
  99     {
 100         chomp $line;
 101
 102         if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
 103         {
 104             $result->{ "Q_ID" }  = $1;
 105             $result->{ "Q_LEN" } = $2;
 106             $result->{ "DIR" }   = "reverse";
 107         }
 108         elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
 109         {
 110             $result->{ "Q_ID" }  = $1;
 111             $result->{ "Q_LEN" } = $2;
 112             $result->{ "DIR" }   = "forward";
 113         }
 114         elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
 115         {
 116             $result->{ "S_ID" }    = $1;
 117             $result->{ "S_BEG" }   = $2 - 1;
 118             $result->{ "Q_BEG" }   = $3 - 1;
 119             $result->{ "HIT_LEN" } = $4;
 120             $result->{ "S_END" }   = $result->{ "S_BEG" } + $result->{ "HIT_LEN" } - 1;
 121             $result->{ "Q_END" }   = $result->{ "Q_BEG" } + $result->{ "HIT_LEN" } - 1;
 122
 123             push @results, dclone $result;
 124         }
 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__