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 vars qw ( @ISA @EXPORT );
44 @ISA = qw( Exporter );
47 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
52 # Martin A. Hansen, June 2007.
54 # Match sequences using MUMmer.
56 my ( $entries1, # FASTA entries
57 $entries2, # FASTA entries
58 $options, # additional MUMmer options - OPTIONAL
59 $tmp_dir, # temporary directory
64 my ( @args, $arg, $file_in1, $file_in2, $cmd, $file_out, $fh, $line, $result, @results );
66 $tmp_dir ||= $ENV{ "BP_TMP" };
68 $options->{ "word_size" } ||= 20;
69 $options->{ "direction" } ||= "both";
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/;
80 $arg = join " ", @args;
82 $file_in1 = "$tmp_dir/muscle1.tmp";
83 $file_in2 = "$tmp_dir/muscle2.tmp";
84 $file_out = "$tmp_dir/muscle3.tmp";
86 map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries1 };
87 map { $_->[ 0 ] =~ tr/ /_/ } @{ $entries2 };
89 Maasha::Fasta::put_entries( $entries1, $file_in1 );
90 Maasha::Fasta::put_entries( $entries2, $file_in2 );
92 Maasha::Common::run( "mummer", "$arg $file_in1 $file_in2 > $file_out 2>/dev/null" );
94 $fh = Maasha::Common::read_open( $file_out );
96 while ( $line = <$fh> )
100 if ( $line =~ /^> (.+)Reverse\s+Len = (\d+)$/ )
102 $result->{ "Q_ID" } = $1;
103 $result->{ "Q_LEN" } = $2;
104 $result->{ "DIR" } = "reverse";
106 elsif ( $line =~ /^> (.+)Len = (\d+)$/ )
108 $result->{ "Q_ID" } = $1;
109 $result->{ "Q_LEN" } = $2;
110 $result->{ "DIR" } = "forward";
112 elsif ( $line =~ /^\s*(.\S+)\s+(\d+)\s+(\d+)\s+(\d+)$/ )
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;
121 push @results, dclone $result;
130 return wantarray ? @results : \@results;
136 # Martin A. Hansen, April 2008.
138 # Vmatches a list of records against a list of index files and the full
139 # path to the result file is returned.
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
149 my ( $query_file, $result_file, @result_files, $fh_in, $fh_out, $line, @fields, $i, $record, $vmatch_args, @index_names, @seq_names, $count_list );
151 $query_file = "$tmp_dir/query.seq";
152 $result_file = "$tmp_dir/vmatch.out";
154 $fh_out = Maasha::Common::write_open( $query_file );
156 foreach $record ( @{ $records } )
158 if ( $record->{ "SEQ_NAME" } and $record->{ "SEQ" } )
160 next if length $record->{ "SEQ" } < 12; # assuming that the index is created for 12 as minimum length
162 push @seq_names, $record->{ "SEQ_NAME" };
164 Maasha::Fasta::put_entry( [ $record->{ "SEQ_NAME" }, $record->{ "SEQ" } ], $fh_out, 80 );
170 if ( $options->{ 'genome' } ) {
171 $vmatch_args = "-complete -d -p -q $query_file";
173 $vmatch_args = "-complete -d -p -showdesc 100 -q $query_file";
176 $vmatch_args .= " -h " . $options->{ "hamming_dist" } if $options->{ "hamming_dist" };
177 $vmatch_args .= " -e " . $options->{ "edit_dist" } if $options->{ "edit_dist" };
179 for ( $i = 0; $i < @{ $index_files }; $i++ )
181 Maasha::Common::run( "vmatch", "$vmatch_args $index_files->[ $i ] > $result_file.$i" );
183 push @result_files, "$result_file.$i";
188 $count_list = vmatch_count_hits( \@result_files ) if ( $options->{ "count" } );
190 $fh_out = Maasha::Common::write_open( $result_file );
192 for ( $i = 0; $i < @{ $index_files }; $i++ )
194 $index_files->[ $i ] =~ s/.+\/(.+)\.fna$/$1/ if $options->{ 'genome' };
196 $fh_in = Maasha::Common::read_open( "$result_file.$i" );
198 while ( $line = <$fh_in> )
202 next if $line =~ /^#/;
204 @fields = split " ", $line;
206 next if $options->{ "max_hits" } and $count_list->[ $fields[ 5 ] ] > $options->{ 'max_hits' };
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
212 print $fh_out join( "\t", @fields ), "\n";
217 unlink "$result_file.$i";
226 sub vmatch_count_hits
228 # Martin A. Hansen, April 2008.
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.
234 my ( $files, # vmatch result files
239 my ( $file, $fh_in, $line, @fields, @count_list );
241 foreach $file ( @{ $files } )
243 $fh_in = Maasha::Common::read_open( $file );
245 while ( $line = <$fh_in> )
249 next if $line =~ /^#/;
251 @fields = split " ", $line;
253 $count_list[ $fields[ 5 ] ]++;
259 return wantarray ? @count_list : \@count_list;
265 # Martin A. Hansen, January 2008.
267 # Parses vmatch output records.
269 my ( $fh, # file handle to vmatch result file.
274 my ( $line, @fields, %record );
276 while ( $line = <$fh> )
280 next if $line =~ /^#/;
282 @fields = split "\t", $line;
284 $record{ "REC_TYPE" } = "VMATCH";
286 $record{ "S_LEN" } = $fields[ 0 ];
287 $record{ "S_ID" } = $fields[ 1 ];
288 $record{ "S_BEG" } = $fields[ 2 ];
290 if ( $fields[ 3 ] eq "D" ) {
291 $record{ "STRAND" } = "+";
293 $record{ "STRAND" } = "-";
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 ];
304 $record{ "Q_END" } = $record{ "Q_BEG" } + $record{ "Q_LEN" } - 1;
305 $record{ "S_END" } = $record{ "S_BEG" } + $record{ "S_LEN" } - 1;
307 return wantarray ? %record : \%record;
314 # Martin A. Hansen, July 2008.
316 # Use mkvtree to create a vmatch index of a given file
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
326 my ( $fh, $entry, $tmp_file );
328 Maasha::Common::dir_create_if_not_exists( $dst_dir );
330 if ( Maasha::Common::file_size( $file ) < 250_000_000 )
332 &Maasha::Common::run( "mkvtree", "-db $src_dir/$file -dna -pl -allout -indexname $dst_dir/$file > /dev/null 3>&1" );
336 $fh = Maasha::Common::read_open( "$src_dir/$file" );
338 while ( $entry = Maasha::Fasta::get_entry( $fh ) )
340 $tmp_file = $entry->[ SEQ_NAME ] . ".fna";
342 Maasha::Fasta::put_entries( [ $entry ], "$tmp_dir/$tmp_file" )
344 &Maasha::Common::run( "mkvtree", "-db $tmp_dir/$tmp_file -dna -pl -allout -indexname $dst_dir/$tmp_file > /dev/null 3>&1" );
346 unlink "$tmp_dir/$tmp_file";
354 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<