Skip empty FASTA sequences, and allow to force taxids

.gitignore

@@ -0,0 +1,7 @@
+/install/
+/Debug/
+/tests/dbs
+/tests/data
+/tests/install
+*.dSYM
+\.idea/

README.md

@@ -1,10 +1,45 @@
-Kraken taxonomic sequence classification system
+KrakenHLL taxonomic sequence classification system with unique k-mer counting
 ===============================================
 
-Please see the [Kraken webpage] or the [Kraken manual]
-for information on installing and operating Kraken.
-A local copy of the [Kraken manual] is also present here
-in the `docs/` directory (`MANUAL.html` and `MANUAL.markdown`).
+[Kraken](https://github.com/DerrickWood/kraken) is a fast taxonomic classifier for metagenomics data. This project, kraken-hll, adds some additional functionality - most notably a unique k-mer count using the HyperLogLog algorithm. Spurious identifications due to sequence contamination in the dataset or database often leads to many reads, however they usually cover only a small portion of the genome. 
 
-[Kraken webpage]:   http://ccb.jhu.edu/software/kraken/
-[Kraken manual]:    http://ccb.jhu.edu/software/kraken/MANUAL.html
+KrakenHLL adds two additional columns to the Kraken report - total number of k-mers observed for taxon, and the total number of unique k-mers observed for taxon (columns 3 and 4, resp.). 
+
+Here's a small example of a classification against a viral database with k=25. There are three species identified by just one read - Enterobacteria phage BP-4795, Salmonella phage SEN22, Sulfolobus monocaudavirus SMV1. Out of those, the identification of Salmonella phage SEN22 is the strongest, as there read was matched with 116 k-mers that are unique to the sequence, while the match to Sulfolobus monocaudavirus SMV1 is only based on a single 25-mer.
+
+```
+99.0958 2192    2192    255510  272869  no rank 0   unclassified
+0.904159    20  0   2361    2318    no rank 1   root
+0.904159    20  0   2361    2318    superkingdom    10239     Viruses
+0.904159    20  0   2361    2318    no rank 35237       dsDNA viruses, no RNA stage
+0.768535    17  0   2074    2063    order   548681        Herpesvirales
+0.768535    17  0   2074    2063    family  10292           Herpesviridae
+0.768535    17  0   2074    2063    subfamily   10374             Gammaherpesvirinae
+0.768535    17  0   2074    2063    genus   10375               Lymphocryptovirus
+0.768535    17  16  2001    1987    species 10376                 Human gammaherpesvirus 4
+0.045208    1   1   4   4   sequence    1000041143                  KC207814.1 Human herpesvirus 4 strain Mutu, complete genome
+0.0904159   2   0   254 254 order   28883         Caudovirales
+0.045208    1   0   28  28  family  10699           Siphoviridae
+0.045208    1   0   28  28  genus   186765            Lambdavirus
+0.045208    1   0   28  28  no rank 335795              unclassified Lambda-like viruses
+0.045208    1   1   28  28  species 196242                Enterobacteria phage BP-4795
+0.045208    1   0   116 116 family  10744           Podoviridae
+0.045208    1   0   116 116 no rank 196895            unclassified Podoviridae
+0.045208    1   0   116 116 no rank 1758253             Escherichia phage phi191 sensu lato
+0.045208    1   1   116 116 species 1647458               Salmonella phage SEN22
+0.045208    1   0   1   1   no rank 51368         unclassified dsDNA viruses
+0.045208    1   1   1   1   species 1351702         Sulfolobus monocaudavirus SMV1
+```
+
+## Usage
+
+For usage, see `krakenhll --help`. Note that you can use the same database as Kraken with one difference - instead of the files `DB_DIR/taxonomy/nodes.dmp` and `DB_DIR/taxonomy/names.dmp` than kraken relies upon, `kraken-hll` needs the file `DB_DIR/taxDB`. This can be generated with the script `build_taxdb`: `KRAKEN_DIR/build_taxdb DB_DIR/taxonomy/names.dmp DB_DIR/taxonomy/nodes.dmp > DB_DIR/taxDB`. The code behind the taxDB is based on [k-SLAM](https://github.com/aindj/k-SLAM).
+
+### Differences to `kraken`
+ - Use `krakenhll --report-file FILENAME ...` to write the kraken report to `FILENAME`.
+ - Use `krakenhll --db DB1 --db DB2 --db DB3 ...` to first attempt, for each k-mer, to assign it based on DB1, then DB2, then DB3. You can use this to prefer identifications based on DB1 (e.g. human and contaminant sequences), then DB2 (e.g. completed bacterial genomes), then DB3, etc. Note that this option is incompatible with `krakenhll-build --generate-taxonomy-ids-for-sequences` since the taxDB between the databases has to be absolutely the same.
+ - Add a suffix `.gz` to output files to generate gzipped output files
+
+### Differences to `kraken-build`
+ - Use `krakenhll-build --generate-taxonomy-ids-for-sequences ...` to add pseudo-taxonomy IDs for each sequence header. An example for the result using this is in the ouput above - one read has been assigned specifically to `KC207814.1 Human herpesvirus 4 strain Mutu, complete genome`.
+ - `seqid2taxid.map` mapping sequence IDs to taxonomy IDs does NOT parse or require `>gi|`, but rather the sequence ID is the header up to just before the first space

VERSION

@@ -0,0 +1 @@
+0.1

install_kraken.sh → install_krakenhll.sh

@@ -1,9 +1,8 @@
 #!/bin/bash
 
+# Portions (c) 2017, Florian Breitwieser <[email protected]>
 # Copyright 2013-2015, Derrick Wood <[email protected]>
 #
-# This file is part of the Kraken taxonomic classification system.
-#
 # Kraken is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
@@ -19,11 +18,22 @@
 
 set -e
 
-VERSION="0.10.6-unreleased"
+DIR=$(dirname $0)
+VERSION=`cat $(dirname $0)/VERSION`
+
+if [ "$1" == "--install-jellyfish" ]; then
+ INSTALL_JELLYFISH=1;
+ shift;
+fi
 
 if [ -z "$1" ] || [ -n "$2" ]
 then
-  echo "Usage: $(basename $0) KRAKEN_DIR"
+  echo "Usage: $(basename $0) [--install-jellyfish] KRAKEN_DIR
+
+If --install-jellyfish is specified, the source code for version 1.1
+is downloaded from http://www.cbcb.umd.edu/software/jellyfish and installed 
+in KRAKEN_DIR. Note that this may overwrite other jellyfish installation in 
+the same path."
   exit 64
 fi
 
@@ -34,13 +44,32 @@ then
   exit 1
 fi
 
+
 # Perl cmd used to canonicalize dirname - "readlink -f" doesn't work
 # on OS X.
 export KRAKEN_DIR=$(perl -MCwd=abs_path -le 'print abs_path(shift)' "$1")
 
 mkdir -p "$KRAKEN_DIR"
-make -C src install
-for file in scripts/*
+if [ "$INSTALL_JELLYFISH" == "1" ]; then
+  WD=`pwd`
+  cd $KRAKEN_DIR
+  if [[ ! -d jellyfish ]]; then
+    wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.11.tar.gz
+    tar xvvf jellyfish-1.1.11.tar.gz
+    mv jellyfish-1.1.11 jellyfish
+  fi
+  cd jellyfish
+  [[ -f Makefile ]] || ./configure
+  make
+  #make install ## doest not work for me on OSX
+  #cp $KRAKEN_DIR/jellyfish-install/bin/jellyfish $KRAKEN_DIR
+  #rm -r jellyfish-1.1.11.tar.gz jellyfish-1.1.11
+  cd $WD
+fi
+
+#make -C src clean
+make -C $DIR/src install
+for file in $DIR/scripts/*
 do
   perl -pl -e 'BEGIN { while (@ARGV) { $_ = shift; ($k,$v) = split /=/, $_, 2; $H{$k} = $v } }'\
            -e 's/#####=(\w+)=#####/$H{$1}/g' \
@@ -52,12 +81,16 @@ do
   fi
 done
 
-echo
-echo "Kraken installation complete."
-echo
-echo "To make things easier for you, you may want to copy/symlink the following"
-echo "files into a directory in your PATH:"
-for file in $KRAKEN_DIR/kraken*
+echo -n "
+Kraken installation complete.
+
+To make things easier for you, you may want to copy/symlink the following
+files into a directory in your PATH:
+
+ln -s"
+for file in $KRAKEN_DIR/krakenhll*
 do
-  [ -x "$file" ] && echo "  $file"
+  [ -x "$file" ] && echo -n " $file"
 done
+echo " DEST_DIR"
+exit 0