Improved BLAST parameters and mini-map output

bin/abricate

@@ -27,7 +27,7 @@ my @BLAST_FIELDS = qw(
 );
 
 my @COLUMNS = qw(
- FILE SEQUENCE START END GENE COVERAGE COVERAGE_MAP %COVERAGE %IDENTITY
+ FILE SEQUENCE START END GENE COVERAGE COVERAGE_MAP GAPS %COVERAGE %IDENTITY
 );
 $COLUMNS[0] = '#'.$COLUMNS[0];
 
@@ -65,10 +65,10 @@ for my $fasta (@ARGV) {
   my %seen;
   my @hit;
   print STDERR "Processing: $fasta\n";
-  my $cmd = "blastn -ungapped -query \Q$fasta\E -db $datadir/resfinder "
-          . " -evalue 1E-99 -dust no -outfmt '$format'"
-#          . " -culling_limit 1"
-          . " -best_hit_overhang 0.49 -best_hit_score_edge 0.0001"
+  my $cmd = "blastn -query \Q$fasta\E -db $datadir/resfinder "
+          . " -task blastn -evalue 1E-99 -dust no -outfmt '$format'"
+          . " -culling_limit 1"
+#          . " -best_hit_overhang 0.49 -best_hit_score_edge 0.0001"
           ;
   print STDERR "Running: $cmd\n" if $verbose;
 
@@ -82,14 +82,16 @@ for my $fasta (@ARGV) {
     next if $seen{ join('~', @hit{qw(qseqid qstart qend)}) }++;
     $hit{sseqid} =~ s/_.*$//;
 #    print STDERR Dumper(\%hit);
-    my $minimap = minimap(@hit{qw(sstart send slen)});
+    my $minimap = minimap( @hit{qw(sstart send slen gaps)} );
+#    $minimap .= '!-'.$hit{gaps} if $hit{gaps} > 0;
     push @hit, [ 
       $fasta, 
       $hit{qseqid}, $hit{qstart}, $hit{qend},
       $hit{sseqid}, 
       $hit{sstart}.'-'.$hit{send}.'/'.$hit{slen}, 
       $minimap, 
-      sprintf("%.2f", 100 * $hit{length} / $hit{slen}),
+      $hit{gaps},
+      sprintf("%.2f", 100 * ($hit{length}-$hit{gaps}) / $hit{slen}),
       $hit{pident},
     ];
   }
@@ -106,8 +108,9 @@ for my $fasta (@ARGV) {
 #----------------------------------------------------------------------
 
 sub minimap {
-  my($x, $y, $L, $scale, $on, $off) = @_;
-  my $WIDTH = 10;
+  my($x, $y, $L, $broken, $scale, $on, $off) = @_;
+  my $WIDTH = 15 - ($broken ? 1 : 0);
+  $broken ||= 0;
   $scale ||= ($L/$WIDTH);
   $on  ||= '=';
   $off ||= '.';
@@ -118,6 +121,7 @@ sub minimap {
   for my $i (0 .. $WIDTH-1) {
 #    print STDERR "$i $x $y $L $scale\n";
     $map .= ($i >= $x and $i <= $y) ? $on : $off;
+    $map .= '/' if $broken and $i==int($WIDTH/2);
   }
   return $map;
 }