Support for HAMAP.hmm.gz

Author: tseemann

bin/prokka

@@ -46,7 +46,7 @@ use Bio::Tools::GuessSeqFormat;
 my @CMDLINE = ($0, @ARGV);
 my $OPSYS = $^O;
 my $EXE = $FindBin::RealScript;
-my $VERSION = "1.15.5";
+my $VERSION = "1.15.6";
 my $AUTHOR = 'Torsten Seemann <torsten.seemann@gmail.com>';
 my $URL = 'https://github.com/tseemann/prokka';
 my $PROKKA_PMID = '24642063';
@@ -1659,8 +1659,15 @@ sub setup_db {
   }
 
   check_tool('hmmpress');
-  for my $hmm (<$dbdir/hmm/*.hmm>) {
+  for my $hmm (<$dbdir/hmm/*.{hmm,hmm.gz}>) {
     msg("Pressing HMM database: $hmm");
+    # We may need to uncompress first
+    if ($hmm =~ m/^(.*?)\.gz$/) {
+      msg("Uncompressing: $hmm");
+      my $hmm_u = $1;
+      runcmd("gzip -dc \Q$hmm\E > \Q$hmm_u\E");
+      $hmm = $hmm_u;
+    }
     runcmd("hmmpress \Q$hmm\E");
   }
 

bin/prokka-genbank_to_fasta_db

@@ -41,7 +41,7 @@ while (my $seq = $in->next_seq) {
 
     next unless $f->length >= $minlen;
 
-    next if !$pseudo and $f->has_tag('pseudo');
+    next if !$pseudo and ($f->has_tag('pseudo') or $f->has_tag('pseudogene'));
 
     my $prod = TAG($f, 'product') || $blank;
 
@@ -54,7 +54,7 @@ while (my $seq = $in->next_seq) {
       $id = TAG($f, $idtag);
       last if $id;
     }
-    $id or die "\nFeature #$counter does not have any of these tags: @TRYIDTAG";
+    $id or die "\nFeature #$counter does not have any of these tags: @TRYIDTAG".Dumper($f);
     $counter++;
     print STDERR "[$counter] ", $seq->display_id, " | $id | $prod\n" if $verbose;
 

bin/prokka-hamap_to_hmm

@@ -6,32 +6,44 @@ use Bio::SeqIO;
 use Data::Dumper;
 use File::Temp qw(tempdir);
 
-my(@Options, $verbose, $datadir, $hypo, $sep, $blank, $pseudo, $minlen);
+my(@Options, $verbose, $datadir, $hypo, $sep, $blank, 
+             $pseudo, $minlen, $seeds, $index, $rules);
 setOptions();
 
 $datadir or die "please specify HAMAP --datadir";
 -d $datadir or die "--datadir $datadir is not a valid folder";
+$seeds or die "need --seeds .tgz file";
+$index or die "need --index .dat file";
+$rules or die "need --rules .dat file";
 
 my $dir = tempdir(CLEANUP=>1);
 print STDERR "Using $dir as temp folder.\n";
 
+
 my %DESC;
-my $ali_fh = get_file("alignment_id.dat", "alignment /gene /product info");
+my $ali_fh = get_file($index, "alignment /gene /product info");
 while (<$ali_fh>) {
+  #0	MF_00001
+  #1	MF_00001
+  #2	28
+  #3	PYRB
+  #4	pyrB
+  #5	Aspartate carbamoyltransferase [pyrB]
   chomp;
   my @x = split m/\t/;
-  next if $x[5] =~ m/UPF\d+|homolog|[@<\[]/;
+  $x[5] =~ s/\s\[\w+\]$//;
+  next if $x[5] =~ m/UPF\d+|homolog/;
   $x[5] =~ s/^probable\s+//gi;
   $DESC{ $x[0] } = "$sep$x[4]$sep$x[5]";
 }
-printf STDERR "Accepted descriptions for %d families\n", scalar keys %DESC;
+printf STDERR "Accepted descriptions for %d families\n", scalar keys %DESC;;
 
-my $fam_fh = get_file("hamap_families.dat", "alignment /EC_number");
+my $fam_fh = get_file($rules, "alignment /EC_number");
 my($AC,$EC,$Bac);
 while (<$fam_fh>) {
   chomp;
   if (m{^//}) {
-    #print STDERR "AC=$AC EC=$EC Bac=$Bac\n";
+    print STDERR "AC=$AC EC=$EC Bac=$Bac\n";
     if ($AC and $EC and exists $DESC{$AC}) {
       $DESC{$AC} = $EC.$DESC{$AC}; 
     }
@@ -53,7 +65,7 @@ while (<$fam_fh>) {
 }
 printf STDERR "%d models made it through filtering\n", scalar keys %DESC;
 
-my $aln_file = "$datadir/hamap_seed_alignments.tar.gz";
+my $aln_file = $seeds;
 if (-r $aln_file) {
   print STDERR "Extracting $aln_file.\n";
   system("tar -C $dir --strip-components 1 -zxf $aln_file")==0 or die "ERROR: $?";
@@ -83,23 +95,24 @@ else {
 }
 
 my $dest = 'HAMAP.hmm';
-print STDERR "Combining into $dest\n";
-system("cat $dir/*.hmm > $dest");
-
+print STDERR "Combining into $dest.ascii\n";
+system("cat $dir/*.hmm > $dest.ascii");
+print STDERR "Converting to binary $dest\m";
+system("hmmconvert -b $dest.ascii > $dest.hmm");
 
 print STDERR "Done.\n";
 
 #----------------------------------------------------------------------
 
 sub get_file {
   my($pattern, $desc) = @_;
-  my($first) = glob("$datadir/$pattern");
+  my($first) = glob("$pattern");
   if ($first and -r $first) {
     print STDERR "Found $desc: $first\n";
     open my $fh, '<', $first;
     return $fh;
   }
-  die "Problem finding $pattern in $datadir ($desc)";
+  die "Problem finding '$pattern'";
 }
 
 #----------------------------------------------------------------------
@@ -111,7 +124,10 @@ sub setOptions {
   @Options = (
     {OPT=>"help",      VAR=>\&usage,                   DESC=>"This help"},
     {OPT=>"verbose!",  VAR=>\$verbose, DEFAULT=>0,     DESC=>"Verbose progress"},
-    {OPT=>"datadir=s", VAR=>\$datadir, DEFAULT=>'',    DESC=>"Path to downloaded HAMAP folder (ftp://ftp.expasy.org/databases/hamap/)" },
+    {OPT=>"datadir=s", VAR=>\$datadir, DEFAULT=>'hamap_seed_alignments',    DESC=>"Path to downloaded HAMAP folder (ftp://ftp.expasy.org/databases/hamap/)" },
+    {OPT=>"seeds=s", VAR=>\$seeds, DEFAULT=>'hamap_seed_alignments.tar.gz',    DESC=>"Seed alignments .tgz" },
+    {OPT=>"index=s", VAR=>\$index, DEFAULT=>'rules_index.dat',    DESC=>"Index file" },
+    {OPT=>"rules=s", VAR=>\$rules, DEFAULT=>'hamap_rules.dat',    DESC=>"Rules file" },
     {OPT=>"sep=s",     VAR=>\$sep,     DEFAULT=>'~~~', DESC=>"Separator between EC/gene/product" },
     {OPT=>"blank=s",   VAR=>\$blank,   DEFAULT=>'',    DESC=>"Replace empty EC/gene/product with this"},
 #    {OPT=>"pseudo!",   VAR=>\$pseudo,  DEFAULT=>0,     DESC=>"Include /pseudo genes"},

db/hmm/HAMAP.hmm

@@ -1,6 +1,7 @@
 
 
 DIR := ../db/kingdom
+DIR_HMM := ../db/hmm
 DIR_B := $(DIR)/Bacteria
 DIR_A := $(DIR)/Archaea
 DIR_V := $(DIR)/Viruses
@@ -14,19 +15,32 @@ CP := cp -vf
 
 .PHONY: install all
 
-
-ALL := AMR sprot_b sprot_v sprot_a
+ALL := AMR sprot_b sprot_v sprot_a HAMAP.hmm.gz
 
 all: $(ALL)
 
 install : $(ALL)
-	$(CP) AMR $(DIR_B)/AMR
+	$(CP) HAMAP.hmm.gz $(DIR_HMM)/
+	$(CP) AMR $(DIR_B)/
 	$(CP) sprot_b $(DIR_B)/sprot
 	$(CP) sprot_a $(DIR_A)/sprot
 	$(CP) sprot_v $(DIR_V)/sprot
 #	$(CP) sprot_m $(DIR_M)/sprot
-	
 
+HAMAP.hmm.gz : HAMAP.hmm
+	pigz -9 < $< > $@
+
+HAMAP.hmm : hamap_rules.dat rules_index.dat hamap_seed_alignments.tar.gz
+	../bin/prokka-hamap_to_hmm 
+
+hamap_rules.dat :
+	wget https://ftp.expasy.org/databases/hamap/hamap_rules.dat
+
+rules_index.dat :
+	wget https://ftp.expasy.org/databases/hamap/rules_index.dat
+	
+hamap_seed_alignments.tar.gz :
+	wget https://ftp.expasy.org/databases/hamap/hamap_seed_alignments.tar.gz
 
 AMR: AMRProt.fa
 	./AMRprot-to-AMR.pl $< > $@
@@ -54,5 +68,3 @@ uniprot_sprot_virsues.dat.gz :
 
 clean :
 	rm -f $(ALL)
-	
-