3 # Copyright (C) 2007 Martin A. Hansen.
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.
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.
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.
19 # http://www.gnu.org/copyleft/gpl.html
22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
25 # Routines to match sequences
28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
33 use Storable qw( dclone );
37 use Maasha::Biopieces;
38 use vars qw ( @ISA @EXPORT );
45 @ISA = qw( Exporter );
48 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
53 # Martin A. Hansen, June 2007.
55 # Match sequences using MUMmer.
57 my ( $entries1, # FASTA entries
58 $entries2, # FASTA entries
59 $options, # additional MUMmer options - OPTIONAL
60 $tmp_dir, # temporary directory
65 my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );
67 $tmp_dir ||= $ENV{ "BP_TMP" };
69 $options->{ "word_size" } ||= 20;
70 $options->{ "direction" } ||= "both";
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/;
81 $arg = join " ", @args;
83 $file_in1 = "$tmp_dir/muscle1.tmp";
84 $file_in2 = "$tmp_dir/muscle2.tmp";
85 $file_out = "$tmp_dir/muscle3.tmp";
87 map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
88 map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };
90 Maasha::Fasta::put_entries( $entries1, $file_in1 );
91 Maasha::Fasta::put_entries( $entries2, $file_in2 );
93 Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );
95 $fh = Maasha::Common::read_open( $file_out );
97 while ( $line = <$fh> )
101 if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
103 $result->{ "Q_ID" } = $1;
104 $result->{ "Q_LEN" } = $2;
105 $result->{ "DIR" } = "reverse";
107 elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
109 $result->{ "Q_ID" } = $1;
110 $result->{ "Q_LEN" } = $2;
111 $result->{ "DIR" } = "forward";
113 elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
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;
122 push @results, dclone $result;
131 return wantarray ? @results : \@results;
137 # Martin A. Hansen, April 2008.
139 # Vmatches a list of records against a list of index files and the full
140 # path to the result file is returned.
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
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 );
152 $query_file = "$tmp_dir/query.seq";
153 $result_file = "$tmp_dir/vmatch.out";
155 $fh_out = Maasha::Common::write_open( $query_file );
157 foreach $record ( @{ $records } )
159 if ( $entry = Maasha::Biopieces::record2fasta( $record ) )
161 next if length $entry->[ SEQ ] < 12; # assuming that the index is created for 12 as minimum length
163 push @seq_names, $entry->[ SEQ_NAME ];
165 Maasha::Fasta::put_entry( $entry, $fh_out, 100 );
171 if ( $options->{ 'genome' } ) {
172 $vmatch_args = "-complete -d -p -q $query_file";
174 $vmatch_args = "-complete -d -p -showdesc 100 -q $query_file";
177 $vmatch_args .= " -h " . $options->{ "hamming_dist" } if $options->{ "hamming_dist" };
178 $vmatch_args .= " -e " . $options->{ "edit_dist" } if $options->{ "edit_dist" };
180 for ( $i = 0; $i < @{ $index_files }; $i++ )
182 Maasha::Common::run( "vmatch", "$vmatch_args $index_files->[ $i ] > $result_file.$i" );
184 push @result_files, "$result_file.$i";
189 $count_list = vmatch_count_hits( \@result_files ) if ( $options->{ "count" } );
191 $fh_out = Maasha::Common::write_open( $result_file );
193 for ( $i = 0; $i < @{ $index_files }; $i++ )
195 $index_files->[ $i ] =~ s/.+\/(.+)\.fna$/$1/ if $options->{ 'genome' };
197 $fh_in = Maasha::Common::read_open( "$result_file.$i" );
199 while ( $line = <$fh_in> )
203 next if $line =~ /^#/;
205 @fields = split " ", $line;
207 next if $options->{ "max_hits" } and $count_list->[ $fields[ 5 ] ] > $options->{ 'max_hits' };
209 $fields[ 1 ] = $index_files->[ $i ]; # S_ID
210 $fields[ 9 ] = $count_list->[ $fields[ 5 ] ] if $options->{ "count" }; # SCORE
211 $fields[ 5 ] = $seq_names[ $fields[ 5 ] ]; # Q_ID
213 print $fh_out join( "\t", @fields ), "\n";
218 unlink "$result_file.$i";
227 sub vmatch_count_hits
229 # Martin A. Hansen, April 2008.
231 # Given a list of Vmatch result file, count duplications based
232 # on q_id. The counts are returned in a list where the list index
233 # corresponds to the q_id index in the query file.
235 my ( $files, # vmatch result files
240 my ( $file, $fh_in, $line, @fields, @count_list );
242 foreach $file ( @{ $files } )
244 $fh_in = Maasha::Common::read_open( $file );
246 while ( $line = <$fh_in> )
250 next if $line =~ /^#/;
252 @fields = split " ", $line;
254 $count_list[ $fields[ 5 ] ]++;
260 return wantarray ? @count_list : \@count_list;
266 # Martin A. Hansen, January 2008.
268 # Parses vmatch output records.
270 my ( $fh, # file handle to vmatch result file.
275 my ( $line, @fields, %record );
277 while ( $line = <$fh> )
281 next if $line =~ /^#/;
283 @fields = split "\t", $line;
285 $record{ "REC_TYPE" } = "VMATCH";
287 $record{ "S_LEN" } = $fields[ 0 ];
288 $record{ "S_ID" } = $fields[ 1 ];
289 $record{ "S_BEG" } = $fields[ 2 ];
291 if ( $fields[ 3 ] eq "D" ) {
292 $record{ "STRAND" } = "+";
294 $record{ "STRAND" } = "-";
297 $record{ "Q_LEN" } = $fields[ 4 ];
298 $record{ "Q_ID" } = $fields[ 5 ];
299 $record{ "Q_BEG" } = $fields[ 6 ];
300 $record{ "MATCH_DIST" } = $fields[ 7 ];
301 $record{ "E_VAL" } = $fields[ 8 ];
302 $record{ "SCORE" } = $fields[ 9 ];
303 $record{ "IDENT" } = $fields[ 10 ];
305 $record{ "Q_END" } = $record{ "Q_BEG" } + $record{ "Q_LEN" } - 1;
306 $record{ "S_END" } = $record{ "S_BEG" } + $record{ "S_LEN" } - 1;
308 return wantarray ? %record : \%record;
315 # Martin A. Hansen, July 2008.
317 # Use mkvtree to create a vmatch index of a given file
319 my ( $file, # FASTA file to index
320 $src_dir, # source directory with FASTA file
321 $dst_dir, # distination directory for Vmatch index
322 $tmp_dir, # temp directory - OPTIONAL
327 my ( $fh, $entry, $tmp_file );
329 Maasha::Common::dir_create_if_not_exists( $dst_dir );
331 # if ( Maasha::Common::file_size( $file ) < 200_000_000 )
333 # &Maasha::Common::run( "mkvtree", "-db $src_dir/$file -dna -pl -allout -indexname $dst_dir/$file > /dev/null 3>&1" );
337 $fh = Maasha::Common::read_open( "$src_dir/$file" );
339 while ( $entry = Maasha::Fasta::get_entry( $fh ) )
341 $tmp_file = $entry->[ SEQ_NAME ] . ".fna";
343 Maasha::Fasta::put_entries( [ $entry ], "$tmp_dir/$tmp_file" )
345 &Maasha::Common::run( "mkvtree", "-db $tmp_dir/$tmp_file -dna -pl -allout -indexname $dst_dir/$tmp_file > /dev/null 3>&1" );
347 unlink "$tmp_dir/$tmp_file";
355 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<