Repository: qiyunzhu/woltka (language: python)

Issue:
Further improved naive solution - right place
#147 should be pushing into this branch...


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index dc82b1f..7d1bf90 100644
--- a/README.md
+++ b/README.md
@@ -27,10 +27,11 @@ Woltka ships with a **QIIME 2 plugin**. [See here for instructions](woltka/q2).
 - Profile tools
   - [Collapse](doc/collapse.md), [Coverage](doc/coverage.md), [Normalize](doc/normalize.md), [Filter](doc/filter.md), [Merge](doc/merge.md)
 - Tutorials
+  - [Sequence alignment](doc/align.md)
   - [Working with WoL](doc/wol.md)
   - [OGU analysis](doc/ogu.md)
 - For users of
-  - [QIIME 2](woltka/q2), [Qiita](doc/qiita.md), [SHOGUN](doc/wol.md#sequence-alignment), [GTDB](doc/gtdb.md), [MetaCyc](doc/metacyc.md), [KEGG](doc/kegg.md)
+  - [QIIME 2](woltka/q2), [Qiita](doc/qiita.md), [SHOGUN](doc/align.md#the-shogun-protocol), [GTDB](doc/gtdb.md), [MetaCyc](doc/metacyc.md), [KEGG](doc/kegg.md)
 - References
   - [Command-line interface](doc/cli.md)
   - [Computational efficiency](doc/perform.md)
@@ -82,9 +83,9 @@ Woltka provides several small test datasets under [woltka/tests/data](woltka/tes
 
 One can execute the following commands to make sure that Woltka functions correctly, and to get an impression of the basic usage of Woltka.
 
-(Note: a more complete list of commands at provided [here](woltka/tests/data). Alternatively, you can skip this test dataset check out the [instructions](doc/wol.md) for working with WoL.)
+(Note: a more complete list of commands at provided [here](woltka/tests/data). Alternatively, you can skip this test dataset and check out the [instructions](doc/wol.md) for working with WoL.)
 
-1\. OGU table generation ([details](doc/ogu.md)):
+1\. OGU (operational genomic unit) table generation ([details](doc/ogu.md)):
 
 ```bash
 woltka classify -i align/bowtie2 -o ogu.biom
@@ -106,11 +107,11 @@ woltka classify \
   -o output_dir
 ```
 
-The mapping file (`taxid.map`) translates genome IDs to taxonomic IDs, which then allows Woltka to classify query sequences based on the NCBI taxonomy (`nodes.dmp` and `names.dmp`).
+The mapping file (`taxid.map`) translates genome IDs to taxonomy IDs, which then allows Woltka to classify query sequences based on the NCBI taxonomy (`nodes.dmp` and `names.dmp`).
 
 The output directory (`output_dir`) will contain three feature tables: `phylum.biom`, `genus.biom` and `species.biom`, each representing a taxonomic profile at one of the three ranks.
 
-3\. Functional profiling by UniRef entries then by GO terms (molecular process):
+3\. Functional profiling by UniRef entries, then by GO molecular processes ([details](doc/ordinal.md)):
 
 ```bash
 woltka classify \
@@ -122,7 +123,7 @@ woltka classify \
   -o output_dir
 ```
 
-Here, the input files are still read-to-genome alignments, instead of read-to-gene ones, but Woltka matches reads to genes based on their coordinates on the genomes (as indicated by the file `coords.txt`). This ensures consistency between taxonomic and functional classifications.
+Here, the input files are still read-to-genome alignments, instead of read-to-gene ones, but Woltka matches reads with genes based on their coordinates on the genomes, using an efficient [algorithm](doc/ordinal.md). The gene coordinates are given by the database file `coords.txt`. The read coordinates are extracted from the alignment files. This ensures consistency between taxonomic and functional classifications.
 
 Subsequently, Woltka is able to assign query sequences to functional units, as defined in mapping files (`uniref.map` and `go/process.tsv`). As you can see, compressed files are supported and auto-detected.
 
diff --git a/doc/align.md b/doc/align.md
new file mode 100644
index 0000000..c93669f
--- /dev/null
+++ b/doc/align.md
@@ -0,0 +1,171 @@
+# Sequence alignment
+
+Woltka takes sequence alignment files as input and generate feature tables. This means that you need to prepare alignments BEFORE running Woltka. The way you want your sequences to be aligned is highly flexible, and can be customized according to your research goals. This page provides general guidelines and detailed examples.
+
+## Contents
+
+- [Sequence data](#sequence-data)
+- [Reference database](#reference-database)
+- [Alignment with Bowtie2](#alignment-with-bowtie2)
+- [Other alignment tools](#other-alignment-tools)
+  - [SHOGUN](#the-shogun-protocol), [Minimap2](#minimap2), [BWA](#bwa), [BURST](#burst), [BLASTn](#blastn)
+
+## Sequence data
+
+This will be your sequence data generated from microbiome samples. They can be DNA (metagenome), RNA (metatranscriptome), or other types of data. They are usually in the file format of FastQ or Fasta, but other formats may also apply.
+
+- **Note**: Woltka is designed for whole-(meta)genome shotgun sequencing data. Technically, you can analyze 16S rRNA data using Woltka, but I don't see the rationale. Please instead learn about [QIIME 2](https://docs.qiime2.org/).
+
+Proper quality control should be performed prior to sequence alignment. Here is an [introduction](https://www.metagenomics.wiki/tools/short-read/quality-control). Also see our latest benchmarks: [Armstrong et al. (2022)](https://journals.asm.org/doi/full/10.1128/msystems.01378-21).
+
+## Reference database
+
+The sequence data should be aligned to a reference database. It can be a database of **reference genomes**. Or, it can be a database of reference genes (nucleotides or proteins). The latter scenario is less typical for Woltka, but it works.
+
+Our [WoL](wol.md) database can be downloaded from this [Dropbox site](https://www.dropbox.com/sh/rwhetdmwomi9k34/AAAKRLIJ2nFniJESgb8lHR-na?dl=0). `wolr1_seqs` is the genome sequences (recommended). `wolr1_orfs` is the re-annotated ORF sequences (amino acids). This database is generally useful for microbiome studies.
+
+We have also tested multiple other databases, including [GEBA-I](https://genome.jgi.doe.gov/portal/geba1003/geba1003.info.html), [UHGG](https://www.ebi.ac.uk/metagenomics/genome-catalogues/human-gut-v2-0), [MIDAS](https://github.com/snayfach/MIDAS), [GTDB](https://gtdb.ecogenomic.org/) (we have a dedicated [tutorial](gtdb.md) for GTDB), and [GEM](https://genome.jgi.doe.gov/portal/GEMs/GEMs.home.html). The results were reported in [Zhu et al. (2022)](https://journals.asm.org/doi/10.1128/msystems.00167-22).
+
+In summary, Woltka works with any database. The choice of database should depend on your research scope, and technical feasibility.
+
+## Alignment with Bowtie2
+
+Woltka has been most extensively tested and validated on sequence alignments generated by [**Bowtie2**](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) ([Langmead and Salzberg, 2012](https://www.nature.com/articles/nmeth.1923)), an ultrafast short sequence aligner that has been widely adopted in many applications.
+
+### Database indexing
+
+Let's say your reference genome sequences are in a Fasta-formatted file `db.fna`. Before running Bowtie2 alignment, you need to index the reference genome database (8 is the number of CPU cores you plan to allocate, same below):
+
+```bash
+bowtie2-build --threads 8 db.fna db
+```
+
+This takes quite a while, and it could cost a substantial amount of memory, especially for large databases (many reference genomes). So do it on a supercomputer. The good news is that you just need to do it once.
+
+For example, our [WoL](wol.md) database, which has 10,575 reference genomes, took nearly 9 hours to build and had a peak memory consumption of 190 GB. The indexed database costs 54.5 GB disk space.
+
+Below is the benchmark we generated for your reference (disclaimer: we make no warranty to its accuracy).
+
+![Bowtie2 indexing](img/bowtie2_build.png)
+
+### Sequence alignment
+
+Then you can align your sequencing data against the database. Let's assume you have paired-end sequences in FastQ files `R1.fq` and `R2.fq` for each sample, or for the entire multiplexed dataset. The basic alignment command is as simple as:
+
+```bash
+bowtie2 -p 8 -x db -1 R1.fq -2 R2.fq -S output.sam
+```
+
+The output file `output.sam` is an alignment file in [SAM](https://en.wikipedia.org/wiki/SAM_(file_format)) format. This will be the [input](input.md) file for Woltka.
+
+As a rule of thumb, the `--very-sensitive` flag is recommended, because we want to maximize discovery of diverse microbes in the sample (instead for pursuing accurate alignment to the human genome). 
+
+Meanwhile, we can suppress SAM header (`--no-head`) and unaligned sequences (`--no-unal`) to reduce the size of output files. Woltka won't need these. But in case you need them for other applications, you may skip these flags.
+
+So the command becomes (**recommended**):
+
+```bash
+bowtie2 -p 8 -x db -1 R1.fq -2 R2.fq -S output.sam --very-sensitive --no-head --no-unal
+```
+
+The alignment step is fast and it costs much less memory compared with the database indexing step. Below is our benchmark on the [HMP](https://www.hmpdacc.org/hmp/) metagenome dataset against our [WoL](wol.md) database (disclaimer: we make no warranty to its accuracy).
+
+- From the benchmark we learned that Bowtie2 scales well with several CPU cores (8 is a good number). But it perhaps makes less sense to use more than 16 CPU cores for a single task.
+
+![Bowtie2 alignment](img/bowtie2_align.png)
+
+It maybe worth checking out the [Bowtie2 manual](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) and this [performing tuning guide](https://community.arm.com/developer/tools-software/hpc/b/hpc-blog/posts/tuning-bowtie2-better-performance) on how to optimize the Bowtie2 framework.
+
+### The SHOGUN protocol
+
+[SHOGUN](https://github.com/knights-lab/SHOGUN) ([Hillmann et al., 2018](https://journals.asm.org/doi/full/10.1128/mSystems.00069-18), [2020](https://academic.oup.com/bioinformatics/article/36/13/4088/5828930)) is a metagenomics package developed by our collaborators and us. It wraps up sequence aligners such as Bowtie2. The team did benchmarks and found that the following additional Bowtie2 parameters are suitable for the task of metagenomics. Our team also validated this in a separate work ([Zhu et al., 2022](https://journals.asm.org/doi/10.1128/msystems.00167-22)).
+
+- In brief, this means that Bowtie2 will return up to 16 hits that have 95% or higher sequence identity with each query sequence. You may consider adding them to the command above.
+
+```bash
+-k 16 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.05"
+```
+
+So the entire command becomes:
+
+```bash
+bowtie2 -p 8 -x db -1 R1.fq -2 R2.fq -S output.sam --very-sensitive --no-head --no-unal -k 16 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.05"
+```
+
+This is equivalent to the SHOGUN command:
+
+```bash
+shogun align -d db -a bowtie2 -t 8 -p 0.95 -i input.fa -o .
+```
+
+
+## Other alignment tools
+
+Woltka has no limitation on what alignment tool should be used. Beyond Bowtie2, it is totally okay to use whatever you have already been using, for consistency with other parts of your research project (which is important!). Below are a few tools which we have tested.
+
+
+### Minimap2
+
+[Minimap2](https://github.com/lh3/minimap2) ([Li, 2018](https://academic.oup.com/bioinformatics/article/34/18/3094/4994778)) is a rapid aligner for short and long DNA and mRNA sequences. The default setting is for long sequences (such as PacBio and Nanopore). To make it work the best with short reads (such as Illumina), you need to add the `-ax sr` parameter.
+
+Database indexing:
+
+```bash
+minimap2 -ax sr db.fna -d db.mmi
+```
+
+Sequence alignment:
+
+```bash
+minimap2 -t 8 -ax sr db R1.fq R2.fq -o output.sam
+```
+
+### BWA
+
+[BWA](https://github.com/lh3/bwa) ([Li and Durbin, 2009](https://academic.oup.com/bioinformatics/article/25/14/1754/225615)) (Burrows-Wheeler Aligner) is a classical short sequence aligner. It has been used many. Researchers of modern time may want to check out [Minimap2](https://github.com/lh3/minimap2) or [BWA-MEM2](https://github.com/bwa-mem2/bwa-mem2). Anyway, here is the usage:
+
+Database indexing:
+
+```bash
+bwa index -p db db.fna
+```
+
+Sequence alignment:
+
+```bash
+bwa mem -t 8 db R1.fq R2.fq > output.sam
+```
+
+### BURST
+
+[**BURST**](https://github.com/knights-lab/BURST) ([Al-Ghalith and Knights, 2020](https://www.biorxiv.org/content/10.1101/2020.09.08.287128v1)) is a fast and optimal sequence aligner that scales well with modern metagenomic datasets while enjoying high accuracy. It is **exhaustive**, as it guarantees to return the best hits from the entire database (unlike tools such as Bowtie2 which is heuristic). You may consider it as a BLASTn alternative. Not as fast as Bowtie2, though, but a good choice when alignment accuracy is the priority.
+
+Database building (these are the parameters we tested):
+
+```bash
+burst -t 8 -r db.fna -o db.edx -a db.acx -d DNA 320 -i 0.95 -s 1500
+```
+
+Sequence alignment:
+
+```bash
+burst -t 8 -r db.edx -a db.acx -q input.fa -o output.b6o -sa -fr -i 0.95
+```
+
+The output file is in [BLAST](https://www.metagenomics.wiki/tools/blast/blastn-output-format-6) format.
+
+### BLASTn
+
+The old-school [**BLAST**](https://en.wikipedia.org/wiki/BLAST_(biotechnology)) algorithm ([Altschul et al., 1990](https://pubmed.ncbi.nlm.nih.gov/2231712/)) is still viable as long as your database / sequence data are of limited size. Despite being slower than modern algorithms, its alignment accuracy is probably the best.
+
+Database building:
+
+```bash
+makeblastdb -in db.fa -dbtype nucl -out db -title db
+```
+
+Sequence alignment (parameters are for reference only):
+
+```bash
+blastn -num_threads 8 -db db -query input.fa -out output.b6o -outfmt 6 -evalue 1e-10 -max_target_seqs 25
+```
diff --git a/doc/classify.md b/doc/classify.md
index 937e83c..d34bbd0 100644
--- a/doc/classify.md
+++ b/doc/classify.md
@@ -55,11 +55,11 @@ In many cases a query sequence has matches in multiple reference sequences, and
 
 ### 1. Keep them all, and divide (default)
 
-In the resulting profile, each subject receives 1 / _k_ count, where _k_ is the total number of subjects of the current query.
+In the resulting profile, each subject receives 1 / _k_ count, where _k_ is the total number of subjects of the current query. In this way, the total feature frequency per sample is consistent with the number of reads in the input alignment.
 
-For example, sequence A was aligned to five genomes: two under genus [_Escherichia_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=561&lvl=3&lin=f&keep=1&srchmode=1&unlock), one under each of genera [_Salmonella_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=590&lvl=3&lin=f&keep=1&srchmode=1&unlock), [_Klebsiella_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=570&lvl=3&lin=f&keep=1&srchmode=1&unlock), and [_Enterobacter_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=547&lvl=3&lin=f&keep=1&srchmode=1&unlock). So at the **genus** level, _Escherichia_ receives 2/5 count, and each of the other three receives 1/5 each.
+For example, sequence A was aligned to five genomes: two under genus [_Escherichia_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=561&lvl=3&lin=f&keep=1&srchmode=1&unlock), one under each of genera [_Salmonella_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=590&lvl=3&lin=f&keep=1&srchmode=1&unlock), [_Klebsiella_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=570&lvl=3&lin=f&keep=1&srchmode=1&unlock), and [_Enterobacter_](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=547&lvl=3&lin=f&keep=1&srchmode=1&unlock). In the **genus** level classification by Woltka, _Escherichia_ receives 2/5 count, and each of the other three receives 1/5 each.
 
-In the read-to-feature maps (`--outmap`), a multi-assignment will be reported as:
+In the read-to-feature maps (`--outmap`), a multi-assignment is reported as:
 
 ```
 A <tab> Escherichia:2 <tab> Salmonella:1 <tab> Klebsiella:1 <tab> Enterobacter:1
diff --git a/doc/faq.md b/doc/faq.md
index 2220300..936cace 100644
--- a/doc/faq.md
+++ b/doc/faq.md
@@ -1,20 +1,100 @@
 # Frequently asked questions
 
+
+## Computing
+
+### How many CPU cores does Woltka use?
+
+Woltka works the best with two CPU cores (threads): one for file decompression and the other for classification. This happens automatically.
+
+### Is Woltka deterministic or stochastic?
+
+Woltka is **deterministic**. Given the same input files and parameters, it always produces the identical output files. There is no "seed" parameter.
+
+### Is Woltka exhaustive/exact/optimal? Or is it approximate/heuristic?
+
+The former. Woltka **exhaustively** captures all valid matches from the alignment file(s).
+
+### Are Woltka results consistent across versions?
+
+To date, all Woltka versions (0.1.0 to 0.1.4) generate **identical** output files given the same setting. Later versions are more efficient and have more features, though.
+
+
 ## Input files
 
+### I ran `woltka classify -i input.fastq ...`, and got an error saying it cannot determine alignment file format. Why?
+
+Woltka takes alignment files as input, NOT original sequencing data (FASTQ, FASTA, etc.). You need to perform alignment on the sequencing data by yourself, such as:
+
+```bash
+bowtie2 -x db -f input.fastq -S output.sam
+```
+
+Then feed the resulting alignment(s) into Woltka.
+
+```bash
+woltka classify -i output.sam ... -o output.tsv
+```
+
+See [here](align.md) for details.
+
+### I ran `woltka classify -i S01.sam ...`. The output feature table has a single column with an empty header. Why?
+
+Woltka is designed to deal with multiple samples at once. If the input is a single alignment file, Woltka automatically treats it as a multiplexed file and attempts to extract individual samples out of it. But in case it is actually not -- Woltka will leave the sample ID empty.
+
+If you really wants Woltka to process one sample at a time, you can do this (aside from manually adding sample ID to the output table):
+
+```bash
+woltka classify --no-demux -i S01.sam ... -o S01.tsv
+```
+
+The `--no-demux` flag will tell Woltka not to try to demultiplex the alignment file. Instead, it will use the filename `S01` as the only sample ID.
+
+See [here](input.md#demultiplexing) for details.
+
 ### Can Woltka parse compressed files?
 
-Yes. All input files for Woltka (alignments and databases) can be supplied as compressed in gzip, bzip2 or xz formates. Woltka will automatically recognize and parse them.
+Yes. All input files for Woltka (alignments and databases) can be supplied as compressed in gzip, bzip2 or xz formats. Woltka will automatically recognize and parse them.
+
+
+## Output files
+
+### Are cell values in the output feature table absolute or relative?
+
+By default, cell values are feature _frequencies_, i.e., the numbers of sequencing reads assigned to individual features. Therefore they are **absolute** abundance.
+
+Woltka provides multiple normalization features. For example, if you want to get **relative** abundance (fraction) instead of absolute abundance, you can do:
+
+```bash
+woltka classify ... --frac
+```
+
+See [here](normalize.md) for details.
+
+### I ran Woltka separately on multiple subsets of data. Can I merge the results?
+
+Yes. The `woltka tools merge` command is for you. See [here](merge.md) for details.
+
 
 ## Hierarchies
 
+### What if a read is aligned to multiple reference sequences?
+
+By default, Woltka will keep all matches and divide them by the number of matches. Meanwhile, Woltka lets the user choose from multiple alternative behaviors.
+
+See [here](classify.md#ambiguous-assignment) for details.
+
 ### What if some subjects are not found in the classification system?
 
-It is okay. Tolerating them and reporting them as "unassigned" is by design and not a bug in Woltka.
+They will not be counted in the output feature table, unless you specify the `--unassigned` flag, in which case an extra feature `unassigned` will be appear.
 
+See [here](classify.md#unassigned-sequences) for details.
 
-## Computing
 
-### How many CPU cores can Woltka use?
+## Stratification
+
+### I attempted to collapse a stratified feature table, and the output is empty?
+
+In this case, you will need to specify which field of a stratified feature should be collapsed, using parameter `--field` followed by field index (starting from 1), otherwise the program will try to collapse the entire feature.
 
-Woltka works the best with two CPU cores (threads): one for decompression and the other for classification. This is automatically determined.
\ No newline at end of file
+See [here](collapse.md#stratification) for details.
diff --git a/doc/gtdb.md b/doc/gtdb.md
index 8d565f0..65c0b4d 100644
--- a/doc/gtdb.md
+++ b/doc/gtdb.md
@@ -1,6 +1,6 @@
 # Working with GTDB
 
-[**GTDB**](https://gtdb.ecogenomic.org/) (Genome Taxonomy Database) ([Parks et al., 2018](https://www.nature.com/articles/nbt.4229)) is a standardized taxonomy system for bacteria and archaea. It is based on the phylogenetic trees of genomes, hence reflecting the evolutionary relationships among microorganisms better, compared to conventional taxonomy. Our research shows that GTDB shares high consistency with the [WoL](https://biocore.github.io/wol/) phylogeny (which was built using more marker genes and more expensive methods) ([Zhu et al., 2019](https://www.nature.com/articles/s41467-019-13443-4)).
+**GTDB** (https://gtdb.ecogenomic.org/) ([Parks et al., 2018](https://www.nature.com/articles/nbt.4229)) is a standardized taxonomy system for bacteria and archaea. It is based on the phylogenetic trees of genomes, hence reflecting the evolutionary relationships among microorganisms better, compared to conventional taxonomy. Our research shows that GTDB shares high consistency with the [WoL](https://biocore.github.io/wol/) phylogeny (which was built using more marker genes and more expensive methods) ([Zhu et al., 2019](https://www.nature.com/articles/s41467-019-13443-4)).
 
 This tutorial will be based on GTDB [**R04-RS89**](https://data.ace.uq.edu.au/public/gtdb/data/releases/release89/) (version 4, indexed to RefSeq release 89), which was the current release at the time of writing. This release contains 143,512 bacterial and 2,392 archaeal organisms in its taxonomy, or 23,458 bacterial and 1,248 archaeal genomes in its collection of actual genome sequences. These 24,706 "**species clusters**" cover the complete GTDB taxonomic framework ([Parks et al., 2020](https://www.nature.com/articles/s41587-020-0501-8)).
 
diff --git a/doc/img/bowtie2_align.png b/doc/img/bowtie2_align.png
new file mode 100644
index 0000000..643193c
Binary files /dev/null and b/doc/img/bowtie2_align.png differ
diff --git a/doc/img/bowtie2_build.png b/doc/img/bowtie2_build.png
new file mode 100644
index 0000000..771b95d
Binary files /dev/null and b/doc/img/bowtie2_build.png differ
diff --git a/doc/input.md b/doc/input.md
index d0efe20..e03836b 100644
--- a/doc/input.md
+++ b/doc/input.md
@@ -33,8 +33,8 @@ S03 S03.sam
 Woltka supports the following alignment formats (specified by parameter `--format` or `-f`):
 
 - `map`: A **simple map** in the format of query \<tab\> subject.
-- `sam`: [**SAM**](https://en.wikipedia.org/wiki/SAM_(file_format)) format. Supported by multiple tools such as Bowtie2 and BWA.
-- `b6o`: **BLAST** tabular format (i.e., BLAST parameter `-outfmt 6`). Supported by multiple tools such as BLAST, DIAMOND, VSEARCH, BURST, etc.
+- `sam`: [**SAM**](https://en.wikipedia.org/wiki/SAM_(file_format)) format. Supported by multiple tools such as Bowtie2, BWA and Minimap2.
+- `b6o`: [**BLAST**](https://www.metagenomics.wiki/tools/blast/blastn-output-format-6) tabular format (i.e., BLAST parameter `-outfmt 6`). Supported by multiple tools such as BLAST, DIAMOND, VSEARCH, BURST, etc.
 
   If not specified, Woltka will automatically infer the format of input alignment files.
 
diff --git a/doc/install.md b/doc/install.md
index d5e52c9..033cbdb 100644
--- a/doc/install.md
+++ b/doc/install.md
@@ -73,7 +73,7 @@ conda install -c conda-forge cython=0.29.6 biom-format=2.1.8
 
 ## Test
 
-You may test whether Woltka functions correctly by running it on some small test datasets. You will need to download the [repository](https://github.com/qiyunzhu/woltka/archive/master.zip) to get them. See [example usage](../README.md#example-usage) for how to run these small tests manually. Alternatively, you may execute the following command in the package directory:
+You may test whether Woltka functions correctly by running it on some small test datasets. You will need to download the [repository](https://github.com/qiyunzhu/woltka/archive/master.zip) to get them. See [example usage](../README.md#example-usage) for how to run these small tests manually. Alternatively, you may execute the following command in the repository directory:
 
 ```bash
 python -m unittest
diff --git a/doc/kegg.md b/doc/kegg.md
index 30b81c1..acb678f 100644
--- a/doc/kegg.md
+++ b/doc/kegg.md
@@ -1,8 +1,8 @@
 # Working with KEGG
 
-**KEGG** (https://www.genome.jp/kegg/) ([Kanehisa et al., 2021](https://academic.oup.com/nar/article/49/D1/D545/5943834)) is the classical database of biological functions. It provides a well-organized hierarchical identification system, such as orthologies (K), modules (M), reactions (R), compounds (C), pathways, diseases and more.
+**KEGG** (https://www.genome.jp/kegg/) ([Kanehisa et al., 2021](https://academic.oup.com/nar/article/49/D1/D545/5943834)) is a classical database of biological functions. It provides a well-organized hierarchical identification system, such as orthologies (K), modules (M), reactions (R), compounds (C), pathways, diseases and more.
 
-Whereas the FTP access to KEGG is limited to subscribed users, the mapping of UniRef entries to KEGG orthology (KO) entries is freely available from the [UniProt](https://www.uniprot.org/downloads) data release. From this point on, we provide a Python script: [**kegg_query.py**](https://github.com/qiyunzhu/utils/blob/main/kegg_query.py) to automatically retrieve higher-level classification information of a given KO list or table from the KEGG server. This is made possible using the official [KEGG REST API](https://www.kegg.jp/kegg/rest/), which is freely available to academic users (however restrictions may apply; see official policy [here](https://www.kegg.jp/kegg/rest/)).
+Whereas the FTP access to KEGG is limited to users with a subscription, the mapping of UniRef entries to KEGG orthology (KO) entries is freely available from the [UniProt](https://www.uniprot.org/downloads) data release. From this point on, we provide a Python script: [**kegg_query.py**](https://github.com/qiyunzhu/utils/blob/main/kegg_query.py) to automatically retrieve higher-level classification information of a given KO list or table from the KEGG server. This is made possible using the official [KEGG REST API](https://www.kegg.jp/kegg/rest/), which is freely available to academic users (however restrictions may apply; see official policy [here](https://www.kegg.jp/kegg/rest/)).
 
 **Step 1**: Classify sequencing data to **KO** entries, bridged by a UniRef-to-KO mapping file:
 
@@ -24,7 +24,7 @@ python kegg_query.py ko.tsv
 
 Be patient, as the KEGG server limits the number of queries per time to 10 (see [policy](https://www.kegg.jp/kegg/rest/keggapi.html#list)).
 
-This will generate multiple of mapping files in the current directory. The file names are self-explanatory. For examples: `ko-to-reaction.txt` is a mapping of KOs to [reactions](https://www.genome.jp/kegg/reaction/) (R), `ko-to-module.txt` is a mapping to [modules](https://www.genome.jp/kegg/module.html) (M), `ko-to-pathway.txt` is to [pathways](https://www.genome.jp/kegg/pathway.html), etc.
+This will generate multiple mapping files in the current directory. The filenames are self-explanatory. For examples: `ko-to-reaction.txt` is a mapping of KOs to [reactions](https://www.genome.jp/kegg/reaction/) (R), `ko-to-module.txt` is a mapping to [modules](https://www.genome.jp/kegg/module.html) (M), `ko-to-pathway.txt` is to [pathways](https://www.genome.jp/kegg/pathway.html), etc.
 
 **Step 3**: Use Woltka's [**collapse**](collapse.md) command to convert KO's to higher-level classification units. This command supports many-to-many mapping, because that's the nature of the relationships between functional units (e.g., genes vs pathways). For example, the following command will generate a profile of **reactions**:
 
diff --git a/doc/merge.md b/doc/merge.md
index 19bf743..606e108 100644
--- a/doc/merge.md
+++ b/doc/merge.md
@@ -6,4 +6,10 @@ The **merge** command merges two or more profiles into one, while treating overl
 woltka tools merge -i input1.biom -i input2.biom -i input3.biom -o output.biom
 ```
 
-The output file from the merge command is **identical** or nearly identical to the output file generated by merging sequence alignment file prior to running Woltka. Small errors (differring by the count of **1**) could be introduced during the normalization of _multiple assignments_ due to floating point arithmetic issues, which is usually not troublesome. In addition to sticking to one-to-one alignments, one can use classification parameters `--rank free`, `--uniq`, `--major`, or `--above` to prevent small errors ([see details](classify.md#ambiguous-assignment)).
+This command is different from QIIME 2's `qiime feature-table merge` command, as the latter can only handle unique sample IDs.
+
+The output file from the merge command is **identical** or nearly identical to the output file generated by merging sequence alignment file prior to running Woltka.
+
+- Very small errors (differring by the count of **1**) could be introduced during the normalization of _multiple assignments_ due to floating point arithmetic issues, which is usually not troublesome. If you do want to avoid this, you may: 1) Use the `--digits` parameter to increase precision level of profiles, which mitigates these errors. 2) Let the upstream alignment step produce unique alignments only (e.g., force one hit per query). 3) Use classification parameters `--rank free`, `--uniq`, `--major`, or `--above` to prevent such small errors ([see details](classify.md#ambiguous-assignment)). 
+
+The **merge** command is generally useful for merging tables. These do not have to be Woltka outputs, nor do they have to be microbiome feature tables. You may supply any feature tables you want to merge.
diff --git a/doc/normalize.md b/doc/normalize.md
index 0e5247f..cf4353a 100644
--- a/doc/normalize.md
+++ b/doc/normalize.md
@@ -2,7 +2,7 @@
 
 By default, the cell values in a feature table (profile) are **counts** (**frequencies**) of each feature in each sample. This behavior maximizes the flexibility of downstream analyses. Meanwhile, Woltka provides multiple functions for normalizing cell values to alternative units and scales:
 
-- `--sizes`: **Divide counts by subject sizes** (e.g., genome / gene lengths).
+- `--sizes`: Divide counts by subject sizes (e.g., by genome or gene lengths).
 - `--frac`: Convert counts into fractions (relative abundances).
 - `--scale`: Scale up all values by a constant factor (e.g., "1k", "1M").
 - `--digits`: Keep certain number of digits after the decimal point.
@@ -26,20 +26,20 @@ Woltka has two modes of normalization, each providing all four functions mention
 
 - **Alternatively**, one can perform normalization on profiles that are already generated (`woltka tools normalize`). This saves the need for rerunning the lengthy classification process.
 
-There is a major advantage of normalization by subject size **during** classification, as detailed below. Other than that, the two modes produce mutually identical results.
+There is a major advantage of normalization by subject size **during** classification, as detailed below. Other than that, the two modes produce mutually identical results. See below for some examples.
 
 
 ## Normalization by subject size
 
-The number of reads assigned to each classification unit can be normalized against a **subject-specific** property, which is usually size / length, but can also be other metrics depending on the specific research.
+The number of reads assigned to each classification unit can be normalized against a **subject-specific** property, which is usually size (length), but can also be other metrics depending on the specific research.
 
-During classification:
+During classification (as part of main workflow):
 
 ```bash
 woltka classify --sizes size.map ...
 ```
 
-On existing profiles:
+Or post classification (on existing profiles):
 
 ```bash
 woltka tools normalize --sizes size.map ...
diff --git a/doc/ogu.md b/doc/ogu.md
index 1be7894..1507cca 100644
--- a/doc/ogu.md
+++ b/doc/ogu.md
@@ -1,8 +1,12 @@
 # OGU analysis
 
-The notion of “**OGU**” (operational genomic unit) is the minimal unit for community ecology studies based on shotgun metagenome or other forms of whole-genome microbiome data. OGUs are simply the reference genomes to which input sequences are aligned. There is no need to assign taxonomy to them. This is in constrast to conventional practices, in which analyses are based on taxonomic units such as genera or species. Therefore, OGU is analogous to ASV in 16S rRNA studies.
+The notion of “**OGU**” (operational genomic unit) is the minimal unit for community ecology studies based on shotgun metagenome or other forms of whole-genome microbiome data. OGUs are simply the reference genomes to which input sequences are aligned. There is no need to assign taxonomy to them. This is in constrast to conventional practices, in which analyses are based on taxonomic units such as genera or species. In this sense, OGU is analogous to ASV in 16S rRNA studies.
 
-The advantage of using OGUs includes 1) highest-possible resolution, 2) independent from taxonomy which is coarse and error-prone as a classification system. 3) allowing for phylogeny-based analysis such as Faith’s PD and UniFrac. The last part is enhanced by the "Web of Life" ([WoL](https://biocore.github.io/wol/)) reference phylogeny.
+The advantage of using OGUs includes:
+
+1. Highest-possible resolution.
+2. Independent from taxonomy, which is coarse and can be error-prone as a classification system.
+3. Allowing for phylogeny-aware analysis such as Faith’s PD and UniFrac. This is enhanced by the "Web of Life" ([WoL](https://biocore.github.io/wol/)) reference phylogeny, or any similar works.
 
 The OGU analysis was explained, benchmarked and discussed in:
 
@@ -20,16 +24,18 @@ The OGU analysis was explained, benchmarked and discussed in:
 
 ## OGU table generation
 
-To generate an OGU table, one needs a multiplexed alignment file, or a directory of per-sample alignment files. These files can be generated by aligning sequencing data against a reference genome database. We recommend using [**SHOGUN**](https://github.com/knights-lab/SHOGUN) with the "Web of Life" database (**WoL**, available for download at: https://biocore.github.io/wol/) (see [details](wol.md)). For example:
+To generate an OGU table, one needs a multiplexed alignment file, or a directory of per-sample alignment files. These files can be generated by aligning sequencing data against a reference genome database. We recommend using Bowtie2 with our "Web of Life" (**WoL**) database (see [details](wol.md)). For example:
 
 ```bash
-shogun align -a bowtie2 -d WoLr1 -i input.fasta -o .
+while read $id; do 
+  bowtie2 -p 8 -x databases/bowtie2/WoLr1 -1 $id.R1.fq -2 $id.R2.fq -S bt2out/$id.sam --very-sensitive --no-head --no-unal
+done < sample.list
 ```
 
 Then one can run Woltka to convert the alignment file(s) into an OGU table:
 
 ```bash
-woltka classify -i alignment.bowtie2.sam -o table.biom
+woltka classify -i bt2out -o table.biom
 ```
 
 The output file `table.biom` is a BIOM table with rows as genome IDs (OGUs), columns as sample IDs, and cell values as counts of OGUs in samples.
diff --git a/doc/perform.md b/doc/perform.md
index 1b87b4e..2a87fbe 100644
--- a/doc/perform.md
+++ b/doc/perform.md
@@ -1,6 +1,6 @@
 # Computational efficiency
 
-Woltka is for handling very large alignment files and complex classification systems. Therefore computational expense is an important consideration. This page explains how much time does Woltka spend on typical datasets and how much memory does it need, as well as multiple tips and tricks for improving performance.
+Woltka is for handling very large alignment files and complex classification systems. Therefore computational expense is an important consideration. This page explains how much time Woltka spends on typical datasets and how much memory it needs, as well as multiple tips and tricks for improving performance.
 
 
 ## Contents
@@ -34,18 +34,20 @@ For the Web of Life ([WoL](https://biocore.github.io/wol/)) database, which incl
 
 In this example, we started with the [CAMI](https://data.cami-challenge.org/) high complexity toy dataset, which contains 5 samples with 15 Gbp HiSeq sequencing data each. We aligned them using the [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) aligner through [SHOGUN](https://github.com/knights-lab/SHOGUN), against the Web of Life ([WoL](https://biocore.github.io/wol/)) database. This step generates up to 16 high-score alignments (matches) for each query sequence.
 
-The five resulting alignment files (SAM format, gzipped) total **38.6 GB** in size. They include **1.05 billion** alignments (lines), with 501 million unique query sequences. They were stored in a hard disk drive (59,000 rpm, SATA 3). The test was performed on with a Xeon E3-1230 v3 CPU (Haswell microarchitecture) and 32 GB DDR3 memory. The software environment was Ubuntu 18.04 plus Python 3.7.
+The five resulting alignment files (SAM format, gzipped) total **38.6 GB** in size. They include **1.05 billion** alignments (lines), with 501 million unique query sequences. They were stored in a hard disk drive (59,000 rpm, SATA 3). The test was performed on with a Xeon E3-1230 v3 CPU (Haswell microarchitecture) and 32 GB DDR3 memory. The software environment was Ubuntu 20.04 plus Python 3.10.
 
 Here are the benchmarks of multiple typical Woltka analyses:
 
 Analysis | Runtime | Memory
 --- | ---: | ---:
-OGU (no classification) | 50:13 | 101.6 MB
-Free-rank taxonomic classification | 1:00:19 | 107.0 MB
-Taxonomic classification at genus | 52:20 | 106.2 MB
-Taxonomic classification at genus, while writing gzipped read maps | 1:05:35 | 106.3 MB
-Taxonomic classification at three ranks: phylum, genus and species | 1:18:26 | 107.5 MB
-Coordinates-based functional classification | 4:24:10 | 19.7 GB
+Operational genomic unit ([OGU](ogu.md)) | 34:39 | 108.4 MB
+Free-rank taxonomic classification | 40:07 | 115.1 MB
+Taxonomic classification at genus | 35:49 | 113.4 MB
+Taxonomic classification at genus, while writing gzipped read maps | 51:59 | 113.0 MB
+Taxonomic classification at three ranks: phylum, genus and species | 52.01 | 115.4 MB
+Matching reads with genes | 4:08:53 | 7.745 GB
+Matching reads with genes, reporting RPK | 4:16:31 | 11.99 GB
+Functional profiling at UniRef entry, then GO process | 4:24:49 | 14.92 GB
 
 
 ## Tips for efficient computing
@@ -114,7 +116,7 @@ For the "simple map" format, Woltka will parse first two columns and query and s
 
 ### Trim SAM file columns
 
-A [SAM](https://en.wikipedia.org/wiki/SAM_(file_format)) file can be bulky. The last several columns may contain the entire sequence and quality scores (e.g., Bowtie2 and Minimap2 do this). These texts are not useful for Woltka. If you haven't opted out this during alignment (e.g., calling Bowtie2 with `--omit-sec-seq`), you can do so afterwards to save disk space:
+A [SAM](https://en.wikipedia.org/wiki/SAM_(file_format)) file can be bulky. The last several columns may contain the entire sequence and quality scores (e.g., Bowtie2 and Minimap2 do this). These texts are not useful for Woltka. You can trim them off to save disk space:
 
 ```bash
 zcat input.sam.gz | grep -v ^@ | cut -f1-9 | sed 's/$/\t*\t*/' | gzip > output.sam.gz
diff --git a/doc/qiita.md b/doc/qiita.md
index c7780c7..642630f 100644
--- a/doc/qiita.md
+++ b/doc/qiita.md
@@ -1,6 +1,6 @@
 # Working with Qiita
 
-[**Qiita**](https://qiita.ucsd.edu/) is a web platform for hosting and analyzing microbiome datasets. Woltka has been implemented in Qiita as part of the standard operating procedure for shotgun metagenomic data analysis. Please find below the instructions for using Woltka and retrieving Woltka results.
+**Qiita** (https://qiita.ucsd.edu/) ([Gonzalez et al., 2018](https://www.nature.com/articles/s41592-018-0141-9)) is a web platform for hosting and analyzing microbiome datasets. Woltka has been implemented in Qiita as part of the standard operating procedure for shotgun metagenomic data analysis. Please find below the instructions for using Woltka and retrieving Woltka results.
 
 ## Existing results
 
@@ -49,4 +49,4 @@ While multiple downstream analyses can be performed on the Woltka output files,
 7. Click "Add command", then click "Run". You job will be queued.
 8. Once the job is completed, you will see a number of output files. Enjoy browsing and interpreting the results!
 
-![Qiita4](img/qiita4.png)
\ No newline at end of file
+![Qiita4](img/qiita4.png)
diff --git a/doc/wol.md b/doc/wol.md
index 74c2fbe..e516334 100644
--- a/doc/wol.md
+++ b/doc/wol.md
@@ -6,71 +6,55 @@
 
 Meanwhile, a basic WoL/Woltka workflow is available from our web-based microbiome study platform: **Qiita** (https://qiita.ucsd.edu/) (see [details](qiita.md)).
 
-## The WoL resource
+## The WoL database
 
-The "Web of Life" (WoL) aims to reconstruct an accurate reference phylogeny for microbial genomes, and to build resources that can benefit microbiome researchers. In phase I ([Zhu et al., 2019](https://www.nature.com/articles/s41467-019-13443-4)), we built a reference [tree](https://biocore.github.io/wol/data/trees/tree.nwk) of 10,575 bacterial and archaeal genomes using 381 marker genes. The basic WoL data release, which contains all necessary files for performing microbiome data analysis, is available for download at the following FTP site:
+The "Web of Life" (WoL) aims to reconstruct an accurate reference phylogeny for microbial genomes, and to build resources that can benefit microbiome researchers. In phase I ([Zhu et al., 2019](https://www.nature.com/articles/s41467-019-13443-4)), we built a reference [tree](https://biocore.github.io/wol/data/trees/tree.nwk) of 10,575 bacterial and archaeal genomes using 381 marker genes. The basic WoL data release, which contains all necessary files for performing microbiome data analysis, is available for download from the following FTP site:
 
-- [ftp://ftp.microbio.me/pub/wol-20April2021](ftp://ftp.microbio.me/pub/wol-20April2021)
+```bash
+wget -r -np -nH --cut-dirs=2 ftp://ftp.microbio.me/pub/wol-20April2021/
+```
+
+Alternatively, you may download the database from the following Dropbox site. Specifically, `wolr1_seqs` is the genome sequences with which you can build a database for sequence alignment, which will happen prior to running Woltka; `wolr1_maps` contains all mapping files which you will need for doing taxonomic, functional and phylogenetic analyses using Woltka.
+
+- https://www.dropbox.com/sh/rwhetdmwomi9k34/AAAKRLIJ2nFniJESgb8lHR-na?dl=0
 
 Meanwhile, the full WoL data release, which also contains raw and processed sequence data, metadata, alternative trees, and pre-built databases for multiple metagenomics tools, are available from our Globus endpoint: [WebOfLife](https://app.globus.org/file-manager/collections/31acbeb8-c62f-11ea-bef9-0e716405a293) (see [instruction](https://biocore.github.io/wol/download#download-via-globus)). The project is detailed at our website: https://biocore.github.io/wol/, including documentation, code, protocols, a gallery and a visualizer.
 
 The following tutorial assumes that you have downloaded the basic or full WoL data release directory. The paths mentioned below are relative to this directory. Specifically, the following directories and files are relevant:
 
-- `databases/shogun/`
 - `databases/bowtie2/`
 - `proteins/coords.txt.xz`
 - `taxonomy/`
 - `function/`
 
-Note: If you downloaded the [basic](ftp://ftp.microbio.me/pub/wol-20April2021) release, you will need to build a Bowtie2 index under `databases/bowtie2` following the instruction provided in the `README` file, and you will need to skip the SHOGUN protocol and use the Bowtie2 protocol (see below).
-
 
-## Sequence alignment
-
-First, you need to align your sequences (namely your FastQ / Fast5 / BAM files) against the WoL database using an aligner of your choice. Let's take [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) for example. Our bioinformatics tool, [SHOGUN](https://github.com/knights-lab/SHOGUN), provides a Bowtie2 wrapper optimized for shotgun metagenomic datasets.
-
-The input file format for SHOGUN needs to be Fasta. If your sequencing data are in FastQ format, you may convert them into Fasta using [seqtk](https://github.com/lh3/seqtk):
-
-```bash
-seqtk seq -a input.fq.gz > input.fa
-```
+## Database indexing
 
-Paired-end FastQ files should be merged (interleaved) followed by conversion:
+The full WoL data release contains a pre-built [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) index which can be directly used (and you can skip section). The basic release only contains the original genome sequences (`concat.fna.xz`). You will need to build a Bowtie2 index by yourself (see [here](align.md#database-indexing) for details):
 
 ```bash
-seqtk mergepe R1.fq.gz R2.fq.gz | seqtk seq -A > input.fa
+mkdir -p databases/bowtie2
+xz -d -k concat.fna.xz
+bowtie2-build --threads 8 concat.fna databases/bowtie2/WoLr1
+rm concat.fna
 ```
 
-Then you can perform sequence alignment using SHOGUN's `align` command:
-
-```bash
-shogun align -d databases/shogun -a bowtie2 -t 16 -p 0.95 -i input.fa -o .
-```
--  The parameter `-p 0.95` represents a sequence identity threshold of 95%, which is the recommended value for shotgun metagenomic data.
+## Sequence alignment
 
-This command is equivalent to (in case you prefer to directly run Bowtie2 without SHOGUN):
+With an indexed database, you can align your sequencing data (e.g., paired-end sequences `R1.fq` and `R2.fq`) against reference genomes using [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) (see [here](align.md#sequence-alignment) for details):
 
 ```bash
-bowtie2 -x databases/bowtie2/WoLr1 -p 16 -f input.fa -S output.sam -k 16 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.05" --very-sensitive --no-head --no-unal
+bowtie2 -p 8 -x db -1 R1.fq -2 R2.fq -S output.sam --very-sensitive --no-head --no-unal
 ```
 
 This step will generate a [SAM](https://en.wikipedia.org/wiki/SAM_(file_format)) format alignment file.
 
-**Note**: You can also run Bowtie2 manually using your choice of parameters, or using other aligners and corresponding databases. Woltka is designed for flexibility. For instance, you can run Bowtie2 to handle paired-end reads as follows:
+**Note**: You can also run Bowtie2 manually using your choice of parameters, or using other aligners and corresponding databases. Woltka is designed for flexibility.
 
-```bash
-bowtie2 -x databases/bowtie2/WoLr1 -p 16 -1 forward.fastq -2 reverse.fastq -S output.sam -k 16 --np 1 --mp "1,1" --rdg "0,1" --rfg "0,1" --score-min "L,0,-0.05" --very-sensitive --no-head --no-unal
-```
-It maybe worth checking out the [Bowtie2 manual](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) and this [performing tuning guide](https://community.arm.com/developer/tools-software/hpc/b/hpc-blog/posts/tuning-bowtie2-better-performance) on how to optimize the Bowtie2 framework.
-
-**Another Note** : It is currently recommended to perform quality filtering before performing mapping. An example of a quality filtering command using [fastp](https://github.com/OpenGene/fastp) is given as follows:
-```bash
-fastp -l 100 -i forward.fastq -I reverse.fastq -o forward_trimmed.fastq -O reverse_trimmed.fastq
-```
 
 ## The OGU analysis
 
-**OGU** (operational genomic unit) ([Zhu et al., 2022](https://journals.asm.org/doi/10.1128/msystems.00167-22)) is a notion we proposed to define the minimum unit of microbiome composition allowed by shotgun metagenomic data. OGUs are reference genomes to which any input sequences have matches. This maximizes the resolution of microbiome composition, and allows for phylogeny-aware analyses using the WoL reference phylogeny. See [details](ogu.md).
+**OGU** (operational genomic unit) ([Zhu et al., 2022](https://journals.asm.org/doi/10.1128/msystems.00167-22)) is a notion we proposed to define the minimum unit of microbiome composition allowed by shotgun metagenomic data. OGUs are reference genomes to which any input sequences have matches. This maximizes the resolution of microbiome composition, and allows for phylogeny-aware analyses using the WoL reference phylogeny. See [here](ogu.md) for details.
 
 In Woltka, an OGU analysis is simply a classification process without a classification system. In such case, the output features will just be the subjects in the alignment file, namely reference genomes.
 
@@ -78,7 +62,7 @@ In Woltka, an OGU analysis is simply a classification process without a classifi
 woltka classify -i input.sam -o ogu.biom
 ```
 
-Woltka accepts a multiplexed alignment file, or a directory of multiple per-sample files, or a mapping of sample IDs to filepaths as input. See [details](input.md#input-filepath). Woltka automatically parses compressed alignment files (such as `input.sam.gz`) so you may [compress](perform.md#compress-alignment-files) them to save disk space.
+Woltka accepts a multiplexed alignment file, or a directory of multiple per-sample files, or a mapping of sample IDs to filepaths as input. See [here](input.md#input-filepath) for details. Woltka automatically parses compressed alignment files (such as `input.sam.gz`) so you may [compress](perform.md#compress-alignment-files) them to save disk space.
 
 The resulting file, `ogu.biom` is a **feature table** describing the composition of the microbiome using OGUs (reference genome IDs, such as `G000123456`). This file can be analyzed using standard microbiome data analysis packages such as [QIIME 2](https://qiime2.org/), in combination with the WoL tree. Here is an [example](ogu.md#ogu-analysis-using-qiime-2).
 
diff --git a/woltka/ordinal.py b/woltka/ordinal.py
index 38d99bb..34a46a6 100644
--- a/woltka/ordinal.py
+++ b/woltka/ordinal.py
@@ -13,6 +13,7 @@
 
 from collections import defaultdict
 from itertools import chain
+from bisect import bisect
 
 from .align import infer_align_format, assign_parser
 
@@ -78,16 +79,12 @@ def flush():
         # master read-to-gene(s) map
         res = defaultdict(set)
 
+        # iterate over nucleotides
         for nucl, locs in locmap.items():
 
-            # merge and sort coordinates
-            # question is to add unsorted read coordinates into already-sorted
-            # gene coordinates
-            # Python's timsort algorithm is efficient for this task
-            try:
-                queue = sorted(chain(coords[nucl], locs))
-
             # it's possible that no gene was annotated on the nucleotide
+            try:
+                glocs = coords[nucl]
             except KeyError:
                 continue
 
@@ -97,11 +94,26 @@ def flush():
             # append prefix if needed
             pfx = nucl + '_' if prefix else ''
 
-            # map reads to genes using the core algorithm
-            for read, gene in match_read_gene(queue):
+            # execute ordinal algorithm when reads are many
+            # 8 (5+ reads) is an empirically determined cutoff
+            if len(locs) > 8:
 
-                # merge read-gene pairs to the master map
-                res[rids[read]].add(pfx + gids[gene])
+                # merge and sort coordinates
+                # question is to add unsorted read coordinates into pre-sorted
+                # gene coordinates
+                # Python's Timsort algorithm is efficient for this task
+                queue = sorted(chain(glocs, locs))
+
+                # map reads to genes using the core algorithm
+                for read, gene in match_read_gene(queue):
+
+                    # add read-gene pairs to the master map
+                    res[rids[read]].add(pfx + gids[gene])
+
+            # execute naive algorithm when reads are few
+            else:
+                for read, gene in match_read_gene_quart(glocs, locs):
+                    res[rids[read]].add(pfx + gids[gene])
 
         # return matching read Ids and gene Ids
         return res.keys(), res.values()
@@ -118,7 +130,7 @@ def flush():
         except (TypeError, IndexError):
             continue
 
-        # skip if length is not available
+        # skip if length is not available or zero
         if not length:
             continue
 
@@ -140,7 +152,9 @@ def flush():
         idx = len(rids)
         rid_append(query)
 
+        # effective length = length * th
         # -int(-x // 1) is equivalent to math.ceil(x) but faster
+        # this value must be >= 1
         locmap[subject].extend((
             (beg << 48) + (-int(-length * th // 1) << 31) + idx,
             (end << 48) + idx))
@@ -362,7 +376,7 @@ def match_read_gene(queue):
             # if code >> 31 & 131071:
             if code & (1 << 31):
 
-                # add it to cache
+                # add its index and coordinate to cache
                 genes[code & (1 << 30) - 1] = code >> 48
 
             # when a gene ends,
@@ -375,6 +389,10 @@ def match_read_gene(queue):
                 for rid, rloc in reads_items():
 
                     # is a match if read/gene overlap is long enough
+                    #   code >> 48:     read end coordinate
+                    #   gloc:           gene start coordinate
+                    #   rloc >> 17`:    read start coordinate
+                    #   rloc & 131071`: effective length - 1
                     if (code >> 48) - max(gloc, rloc >> 17) >= rloc & 131071:
                         yield rid, code & (1 << 30) - 1
 
@@ -384,7 +402,8 @@ def match_read_gene(queue):
             # when a read begins,
             if code >> 31 & 131071:
 
-                # add it and its effective length to cache
+                # add its index, coordinate and effective length - 1 to cache
+                # the latter two are stored as a single integer
                 reads[code & (1 << 30) - 1] = (code >> 31) - 1
 
             # when a read ends,
@@ -394,6 +413,13 @@ def match_read_gene(queue):
                 rloc = reads_pop(code & (1 << 30) - 1)
 
                 # check cached genes
+                # a potential optimization is to pre-calculate `rloc >> 17`
+                #   and `(code >> 48) - (rloc & 131071)`, however, it is not
+                #   worth the overhead in real applications because there is
+                #   usually zero or one gene in cache
+                # another potential optimization is to replace `max` (which is
+                #   a function call with a ternary operator, but one needs the
+                #   first optimization prior to this
                 for gid, gloc in genes_items():
 
                     # same as above
@@ -401,6 +427,246 @@ def match_read_gene(queue):
                         yield code & (1 << 30) - 1, gid
 
 
+def match_read_gene_naive(geneque, readque):
+    """Associate reads with genes using a naive approach, which performs
+    nested iteration over genes and reads.
+
+    Parameters
+    ----------
+    geneque : list of int
+        Sorted queue of genes.
+    readque : list of int
+        Paired queue of reads.
+
+    Yields
+    ------
+    int
+        Read index.
+    int
+        Gene index.
+
+    See Also
+    --------
+    match_read_gene
+
+    Notes
+    -----
+    This is a reference implementation of the naive nested iteration method
+    with limited optimization. It is O(nm), where n and m are the numbers of
+    genes and reads, respectively. It should be much slower than the ordinal
+    method. However, when the number of reads is small, it may be efficient
+    because it saves the sorting step.
+    """
+    # only genes are cached
+    genes = {}
+    genes_pop = genes.pop
+
+    # iterate over reads (paired)
+    it = iter(readque)
+    for x, y in zip(it, it):
+
+        # pre-calculate read metrics
+        rid = x & (1 << 30) - 1     # index
+        beg = x >> 48               # start coordinate
+        end = y >> 48               # end coordinate
+        L = (x >> 31) - 1 & 131071  # effective length - 1
+
+        # iterate over genes (ordinal)
+        for code in geneque:
+
+            # add gene to cache
+            if code & (1 << 31):
+                genes[code & (1 << 30) - 1] = code >> 48
+
+            # check overlap while removing gene from cache:
+            # min(gene end, read end) - max(gene start, read start) >=
+            # effective length - 1
+            elif (min(code >> 48, end) -
+                  max(genes_pop(code & (1 << 30) - 1), beg)) >= L:
+                yield rid, code & (1 << 30) - 1
+
+
+def match_read_gene_quart(geneque, readque):
+    """Associate reads with genes by iterating between read region and nearer
+    end of gene queue.
+
+    Parameters
+    ----------
+    geneque : list of int
+        Sorted queue of genes.
+    readque : list of int
+        Paired queue of reads.
+
+    Yields
+    ------
+    int
+        Read index.
+    int
+        Gene index.
+
+    See Also
+    --------
+    match_read_gene_naive
+
+    Notes
+    -----
+    This optimized method is similar to the naive method, but it reduces the
+    search space from n to n / 4, thus it significantly improves efficiency.
+
+    The idea is that one only needs to iterate between the position where read
+    can fit it, and the nearer end of the gene queue.
+
+    Further reduction of search space is not possible, because genes may span
+    over half of the genome, and they cannot be detected within out iterating
+    to end of the gene queue.
+
+    This method uses bisection to find insertion points of read start in the
+    pre-sorted gene queue, which has O(logn) time.
+    """
+    n = len(geneque)  # entire search space
+    mid = n // 2      # mid point
+
+    # iterate over paired read starts and ends
+    it = iter(readque)
+    for x, y in zip(it, it):
+
+        # pre-calculate read metrics
+        rid = x & (1 << 30) - 1     # index
+        beg = x >> 48               # start coordinate
+        end = y >> 48               # end coordinate
+        L = (x >> 31) - 1 & 131071  # effective length - 1
+
+        # genes starting within read region
+        # will record their coordinates
+        within = {}
+        within_pop = within.pop
+
+        # read starts in left half of gene queue
+        if x < geneque[mid]:
+
+            # genes starting before read region
+            # no need to record coordinates as overlap is not possible
+            before = set()
+            before_add = before.add
+            before_remove = before.remove
+
+            # locate read start using bisection
+            # Python's `bisect` is more efficient than homebrew code
+            i = bisect(geneque, x, hi=mid)
+
+            # iterate over gene positions before read region
+            for code in geneque[:i]:
+
+                # add gene to cache when it starts
+                if code & (1 << 31):
+                    before_add(code & (1 << 30) - 1)
+
+                # drop gene from cache when it ends
+                else:
+                    before_remove(code & (1 << 30) - 1)
+
+            # iterate over gene positions after read start
+            for code in geneque[i:]:
+
+                # stop when exceeding read end
+                if code > y:
+                    break
+
+                # when gene starts, add it and its coordinate to cache
+                elif code & (1 << 31):
+                    within[code & (1 << 30) - 1] = code >> 48
+
+                # when gene ends, check overlap and remove it from cache
+                else:
+                    gid = code & (1 << 30) - 1
+
+                    # most genes should start before read region
+                    try:
+                        before_remove(gid)
+
+                    # if gene started within read region,
+                    # overlap = gene end - gene start
+                    except KeyError:
+                        if (code >> 48) - within_pop(gid) >= L:
+                            yield rid, gid
+
+                    # if gene started before read region,
+                    # overlap = gene end - read start
+                    else:
+                        if (code >> 48) - beg >= L:
+                            yield rid, gid
+
+            # yield all genes that started before read but not yet ended
+            for gid in before:
+                yield rid, gid
+
+            # yield genes that started after read but not yet ended,
+            # until overlap = read end - gene start is too short
+            for gid, gloc in within.items():
+                if end - gloc < L:
+                    break
+                yield rid, gid
+
+        # read starts in right half of gene queue
+        else:
+
+            # genes starting after read region
+            # no need to record coordinates as overlap is not possible
+            after = set()
+            after_add = after.add
+            after_remove = after.remove
+
+            # locate read start using bisection
+            i = bisect(geneque, x, lo=mid + 1, hi=n)
+
+            # iterate over gene positions within read region
+            # for i in range(bisect(geneque, x), n):
+            while i < n:
+                code = geneque[i]
+
+                # stop when exceeding read end
+                if code > y:
+                    break
+
+                # when gene starts, add it and its coordinate to cache
+                if code & (1 << 31):
+                    within[code & (1 << 30) - 1] = code >> 48
+
+                # when gene ends, check overlap and remove it from cache
+                # gene end must be <= read end
+                # if gene not found in cache (meaning gene started before read
+                # region), use read start, otherwise use gene start
+                elif (code >> 48) - within_pop(code & (1 << 30) - 1, beg) >= L:
+                    yield rid, code & (1 << 30) - 1
+
+                # move to next position
+                i += 1
+
+            # iterate over gene positions after read region
+            for code in geneque[i:]:
+
+                # add gene to cache when it starts
+                if code & (1 << 31):
+                    after_add(code & (1 << 30) - 1)
+
+                # when gene ends,
+                else:
+
+                    # most remaining genes should start after read region
+                    try:
+                        after_remove(code & (1 << 30) - 1)
+
+                    # otherwise, the gene spans over entire read region
+                    except KeyError:
+
+                        # check overlap while removing gene from cache
+                        # gene end must be >= read end
+                        # if gene not found in cache, use read start, otherwise
+                        # use gene start
+                        if end - within_pop(code & (1 << 30) - 1, beg) >= L:
+                            yield rid, code & (1 << 30) - 1
+
+
 def match_read_gene_dummy(queue, lens, th):
     """Associate reads with genes based on a sorted queue of coordinates.
     Parameters


Test patch:
diff --git a/woltka/tests/test_ordinal.py b/woltka/tests/test_ordinal.py
index 56f55da..f500bc7 100644
--- a/woltka/tests/test_ordinal.py
+++ b/woltka/tests/test_ordinal.py
@@ -9,7 +9,7 @@
 # ----------------------------------------------------------------------------
 
 """Unit tests for the ordinal system, as well as a demonstration of how the
-system works. See docstrings of `test_match_read_gene`.
+system works. See docstrings of `test_match_read_gene_dummy`.
 """
 
 from unittest import TestCase, main
@@ -22,8 +22,9 @@
 from woltka.file import openzip
 from woltka.align import parse_b6o_line, parse_sam_line
 from woltka.ordinal import (
-    match_read_gene_dummy, match_read_gene, ordinal_parser_dummy,
-    ordinal_mapper, load_gene_coords, calc_gene_lens)
+    match_read_gene_dummy, match_read_gene, match_read_gene_naive,
+    match_read_gene_quart, ordinal_parser_dummy, ordinal_mapper,
+    load_gene_coords, calc_gene_lens)
 
 
 class OrdinalTests(TestCase):
@@ -149,6 +150,91 @@ def test_match_read_gene(self):
                (8, 3)]
         self.assertListEqual(obs, exp)
 
+    def test_match_read_gene_naive(self):
+        """The naive solution should produce identical result compared to
+        the default (ordinal) solution.
+        """
+        genes = [(1,  5, 29),
+                 (2, 33, 61),
+                 (3, 65, 94)]
+        reads = [(1, 10, 29),
+                 (2, 16, 35),
+                 (3, 20, 39),
+                 (4, 22, 41),
+                 (5, 30, 49),
+                 (6, 30, 49),
+                 (7, 60, 79),
+                 (8, 65, 84),
+                 (9, 82, 95)]
+        genes = [x for idx, beg, end in genes for x in (
+            (beg << 48) + (1 << 31) + (1 << 30) + idx,
+            (end << 48) + (0 << 31) + (1 << 30) + idx)]
+        reads = [x for idx, beg, end in reads for x in (
+            (beg << 48) + (14 << 31) + (0 << 30) + idx,
+            (end << 48) + (0 << 31) + (0 << 30) + idx)]
+
+        # don't sort, but directly feed both queues
+        obs = list(match_read_gene_naive(genes, reads))
+        exp = [(1, 1),
+               (2, 1),
+               (5, 2),
+               (6, 2),
+               (7, 3),
+               (8, 3)]
+        self.assertListEqual(obs, exp)
+
+    def test_match_read_gene_quart(self):
+        """The naive solution should produce identical result compared to
+        the default (ordinal) solution.
+        """
+        genes = [(1,  5, 29),
+                 (2, 33, 61),
+                 (3, 65, 94),
+                 (4, 61, 76),  # added a small gene within a read
+                 (5, 68, 72)]  # added a tiny gene
+        reads = [(1, 10, 29),
+                 (2, 16, 35),
+                 (3, 20, 39),
+                 (4, 22, 41),
+                 (5, 30, 49),
+                 (6, 30, 49),
+                 (7, 60, 79),
+                 (8, 65, 84),
+                 (9, 82, 95)]
+        genes = [x for idx, beg, end in genes for x in (
+            (beg << 48) + (1 << 31) + (1 << 30) + idx,
+            (end << 48) + (0 << 31) + (1 << 30) + idx)]
+        genes.sort()
+        reads = [x for idx, beg, end in reads for x in (
+            (beg << 48) + (14 << 31) + (0 << 30) + idx,
+            (end << 48) + (0 << 31) + (0 << 30) + idx)]
+
+        obs = list(match_read_gene_quart(genes, reads))
+        exp = [(1, 1),
+               (2, 1),
+               (5, 2),
+               (6, 2),
+               (7, 4),  # match to added small gene
+               (7, 3),
+               (8, 3)]
+        self.assertListEqual(obs, exp)
+
+        # a special case with a giant read
+        genes = [(1, 1, 5),
+                 (2, 6, 7),
+                 (3, 7, 8)]
+        reads = [(1, 4, 9)]
+        genes = [x for idx, beg, end in genes for x in (
+            (beg << 48) + (1 << 31) + (1 << 30) + idx,
+            (end << 48) + (0 << 31) + (1 << 30) + idx)]
+        genes.sort()
+        reads = [x for idx, beg, end in reads for x in (
+            (beg << 48) + (5 << 31) + (0 << 30) + idx,
+            (end << 48) + (0 << 31) + (0 << 30) + idx)]
+        obs = list(match_read_gene_quart(genes, reads))
+        exp = []
+        self.assertListEqual(obs, exp)
+
     def test_ordinal_parser_dummy(self):
         # b6o (BLAST, DIAMOND, BURST, etc.)
         b6o = (