diff --git a/metagraph/external-libraries/caches b/metagraph/external-libraries/caches
index b969baa3e7..9843c53447 160000
--- a/metagraph/external-libraries/caches
+++ b/metagraph/external-libraries/caches
@@ -1 +1 @@
-Subproject commit b969baa3e774fac964487e14bb0122708b09114b
+Subproject commit 9843c534478a9188343c1cf41fb3e9269df872dd
diff --git a/metagraph/integration_tests/base.py b/metagraph/integration_tests/base.py
index c655bf7055..84f3f5aafe 100644
--- a/metagraph/integration_tests/base.py
+++ b/metagraph/integration_tests/base.py
@@ -106,7 +106,9 @@ def _clean(graph, output, extra_params=''):
     def _annotate_graph(input, graph_path, output, anno_repr,
                         separate=False, no_fork_opt=False, no_anchor_opt=False):
         target_anno = anno_repr
-        if anno_repr in {'row_sparse'} or anno_repr.endswith('brwt') or anno_repr.startswith('row_diff'):
+        if (anno_repr in {'row_sparse', 'column_coord'} or
+                anno_repr.endswith('brwt') or
+                anno_repr.startswith('row_diff')):
             target_anno = anno_repr
             anno_repr = 'column'
         elif anno_repr in {'flat', 'rbfish'}:
@@ -117,6 +119,9 @@ def _annotate_graph(input, graph_path, output, anno_repr,
                     -i {graph_path} --anno-type {anno_repr} \
                     -o {output} {input}'
 
+        if target_anno.endswith('_coord'):
+            command += ' --coordinates'
+
         with_counts = target_anno.endswith('int_brwt')
         if with_counts:
             command += ' --count-kmers'
diff --git a/metagraph/integration_tests/test_align.py b/metagraph/integration_tests/test_align.py
index fa789f5e70..b7fcb25522 100644
--- a/metagraph/integration_tests/test_align.py
+++ b/metagraph/integration_tests/test_align.py
@@ -60,6 +60,66 @@ def test_simple_align_all_graphs(self, representation):
         self.assertEqual(last_split[1], "AACAGAGAATTGTTTAAATTACAATCTTAGCTATGGGTGCTAAAGGTGGAGTTATAGACTTTTTCACTGATTTGTCGTTGGAAAAAGCTTTTCATCTCGGGTTTACAAGTCTGGTGTATTTGTTTATACTAGAAGGACAGGCGCATTTGA")
         self.assertEqual(last_split[4], "22")
 
+    @parameterized.expand(GRAPH_TYPES)
+    def test_simple_align_map_all_graphs(self, representation):
+
+        self._build_graph(input=TEST_DATA_DIR + '/genome.MT.fa',
+                          output=self.tempdir.name + '/genome.MT',
+                          k=11, repr=representation,
+                          extra_params="--mask-dummy")
+
+        res = self._get_stats(self.tempdir.name + '/genome.MT' + graph_file_extension[representation])
+        params_str = res.stdout.decode().split('\n')[2:]
+        self.assertEqual('k: 11', params_str[0])
+        self.assertEqual('nodes (k): 16438', params_str[1])
+        self.assertEqual('mode: basic', params_str[2])
+
+        stats_command = '{exe} align -i {graph} --map --count-kmers {reads}'.format(
+            exe=METAGRAPH,
+            graph=self.tempdir.name + '/genome.MT' + graph_file_extension[representation],
+            reads=TEST_DATA_DIR + '/genome_MT1.fq',
+        )
+        res = subprocess.run(stats_command.split(), stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        params_str = res.stdout.decode().rstrip().split('\n')
+        self.assertEqual(len(params_str), 6)
+        self.assertEqual(params_str[0], 'MT-10/1\t1/140/1')
+        self.assertEqual(params_str[1], 'MT-8/1\t140/140/140')
+        self.assertEqual(params_str[2], 'MT-6/1\t140/140/140')
+        self.assertEqual(params_str[3], 'MT-4/1\t0/140/0')
+        self.assertEqual(params_str[4], 'MT-2/1\t140/140/140')
+        self.assertEqual(params_str[5], 'MT-11/1\t1/140/1')
+
+    @parameterized.expand(GRAPH_TYPES)
+    def test_simple_align_map_canonical_all_graphs(self, representation):
+
+        self._build_graph(input=TEST_DATA_DIR + '/genome.MT.fa',
+                          output=self.tempdir.name + '/genome.MT',
+                          k=11, repr=representation, mode='canonical',
+                          extra_params="--mask-dummy")
+
+        res = self._get_stats(self.tempdir.name + '/genome.MT' + graph_file_extension[representation])
+        params_str = res.stdout.decode().split('\n')[2:]
+        self.assertEqual('k: 11', params_str[0])
+        self.assertEqual('nodes (k): 32782', params_str[1])
+        self.assertEqual('mode: canonical', params_str[2])
+
+        stats_command = '{exe} align -i {graph} --map --count-kmers {reads}'.format(
+            exe=METAGRAPH,
+            graph=self.tempdir.name + '/genome.MT' + graph_file_extension[representation],
+            reads=TEST_DATA_DIR + '/genome_MT1.fq',
+        )
+        res = subprocess.run(stats_command.split(), stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        params_str = res.stdout.decode().rstrip().split('\n')
+        self.assertEqual(len(params_str), 6)
+        self.assertEqual(params_str[0], 'MT-10/1\t140/140/140')
+        self.assertEqual(params_str[1], 'MT-8/1\t140/140/140')
+        self.assertEqual(params_str[2], 'MT-6/1\t140/140/140')
+        self.assertEqual(params_str[3], 'MT-4/1\t129/140/129')
+        self.assertEqual(params_str[4], 'MT-2/1\t140/140/139')
+        self.assertEqual(params_str[5], 'MT-11/1\t2/140/2')
+
     @parameterized.expand(['succinct'])
     def test_simple_align_json_all_graphs(self, representation):
 
diff --git a/metagraph/integration_tests/test_query.py b/metagraph/integration_tests/test_query.py
index 1780680900..159cab88e3 100644
--- a/metagraph/integration_tests/test_query.py
+++ b/metagraph/integration_tests/test_query.py
@@ -6,6 +6,7 @@
 from tempfile import TemporaryDirectory
 import glob
 import os
+import numpy as np
 from helpers import get_test_class_name
 from base import TestingBase, METAGRAPH, TEST_DATA_DIR, graph_file_extension
 
@@ -16,6 +17,7 @@
 PROTEIN_MODE = os.readlink(METAGRAPH).endswith("_Protein")
 
 anno_file_extension = {'column': '.column.annodbg',
+                       'column_coord': '.column_coord.annodbg',
                        'row': '.row.annodbg',
                        'row_diff': '.row_diff.annodbg',
                        'row_sparse': '.row_sparse.annodbg',
@@ -503,6 +505,257 @@ def test_batch_query_with_tiny_batch(self):
         self.assertEqual(res.returncode, 0)
         self.assertEqual(len(res.stdout), 136959)
 
+    def test_query_coordinates(self):
+        if not self.anno_repr.endswith('_coord'):
+            self.skipTest('annotation does not support coordinates')
+
+        query_command = f'{METAGRAPH} query --query-coords \
+                            -i {self.tempdir.name}/graph{graph_file_extension[self.graph_repr]} \
+                            -a {self.tempdir.name}/annotation{anno_file_extension[self.anno_repr]} \
+                            --discovery-fraction 0.05 {TEST_DATA_DIR}/transcripts_100.fa'
+
+        res = subprocess.run(query_command.split(), stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        self.assertEqual(len(res.stdout), 2155983)
+
+        query_command = f'{METAGRAPH} query --query-coords \
+                            -i {self.tempdir.name}/graph{graph_file_extension[self.graph_repr]} \
+                            -a {self.tempdir.name}/annotation{anno_file_extension[self.anno_repr]} \
+                            --discovery-fraction 0.95 {TEST_DATA_DIR}/transcripts_100.fa'
+
+        res = subprocess.run(query_command.split(), stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        self.assertEqual(len(res.stdout), 687712)
+
+
+@parameterized_class(('graph_repr', 'anno_repr'),
+    input_values=product(
+        [repr for repr in GRAPH_TYPES if not (repr == 'bitmap' and PROTEIN_MODE)],
+        ['int_brwt', 'row_diff_int_brwt']
+    ),
+    class_name_func=get_test_class_name
+)
+class TestQueryCounts(TestingBase):
+    @classmethod
+    def setUpClass(cls):
+        cls.tempdir = TemporaryDirectory()
+
+        cls.kmer_counts_1 = {
+            'AAA': 1,
+            'AAC': 2,
+            'ACC': 3,
+            'CCC': 4,
+            'CCG': 5,
+            'CGG': 6,
+            'GGG': 7,
+            'GGT': 8,
+            'GTT': 9,
+            'TTT': 10,
+            'TTA': 11,
+            'TAA': 12,
+        }
+        cls.kmer_counts_2 = {
+            'AAA': 11,
+            'AAC': 12,
+            'ACC': 13,
+            'CCC': 14,
+            'CCG': 15,
+            'CGG': 16,
+            'GGG': 17,
+            'GGT': 18,
+            'GTT': 19,
+            'TTT': 20,
+        }
+        fasta_file = cls.tempdir.name + '/file.fa'
+        with open(fasta_file, 'w') as f:
+            for kmer, count in cls.kmer_counts_1.items():
+                f.write(f'>L1\n{kmer}\n' * count)
+
+            for kmer, count in cls.kmer_counts_2.items():
+                f.write(f'>L2\n{kmer}\n' * count)
+
+        cls.k = 3
+
+        cls.with_bloom = False
+        if cls.graph_repr == 'succinct_bloom':
+            cls.graph_repr = 'succinct'
+            cls.with_bloom = True
+
+        cls.mask_dummy = False
+        if cls.graph_repr == 'succinct_mask':
+            cls.graph_repr = 'succinct'
+            cls.mask_dummy = True
+
+        construct_command = f"{METAGRAPH} build {'--mask-dummy' if cls.mask_dummy else ''} -p {NUM_THREADS} \
+                                --graph {cls.graph_repr} -k {cls.k} -o {cls.tempdir.name}/graph {fasta_file}"
+
+        res = subprocess.run([construct_command], shell=True)
+        assert(res.returncode == 0)
+
+        stats_command = '{exe} stats {graph}'.format(
+            exe=METAGRAPH,
+            graph=cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+        )
+        res = subprocess.run(stats_command.split(), stdout=PIPE)
+        assert(res.returncode == 0)
+        params_str = res.stdout.decode().split('\n')[2:]
+        assert('k: 3' == params_str[0])
+        if cls.graph_repr != 'succinct' or cls.mask_dummy:
+            assert('nodes (k): 12' == params_str[1])
+        assert('mode: basic' == params_str[2])
+
+        if cls.with_bloom:
+            convert_command = '{exe} transform -o {outfile} --initialize-bloom {bloom_param} {input}'.format(
+                exe=METAGRAPH,
+                outfile=cls.tempdir.name + '/graph',
+                bloom_param='--bloom-fpp 0.1',
+                input=cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+            )
+            res = subprocess.run([convert_command], shell=True)
+            assert(res.returncode == 0)
+
+        def check_suffix(anno_repr, suffix):
+            match = anno_repr.endswith(suffix)
+            if match:
+                anno_repr = anno_repr[:-len(suffix)]
+            return anno_repr, match
+
+        cls.anno_repr, separate = check_suffix(cls.anno_repr, '_separate')
+        cls.anno_repr, no_fork_opt = check_suffix(cls.anno_repr, '_no_fork_opt')
+        cls.anno_repr, no_anchor_opt = check_suffix(cls.anno_repr, '_no_anchor_opt')
+
+        cls._annotate_graph(
+            fasta_file,
+            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+            cls.tempdir.name + '/annotation',
+            cls.anno_repr,
+            separate,
+            no_fork_opt,
+            no_anchor_opt
+        )
+
+        # check annotation
+        anno_stats_command = '{exe} stats -a {annotation}'.format(
+            exe=METAGRAPH,
+            annotation=cls.tempdir.name + '/annotation' + anno_file_extension[cls.anno_repr],
+        )
+        res = subprocess.run(anno_stats_command.split(), stdout=PIPE)
+        assert(res.returncode == 0)
+        params_str = res.stdout.decode().split('\n')[2:]
+        assert('labels:  2' == params_str[0])
+        if cls.graph_repr != 'hashfast' and (cls.graph_repr != 'succinct' or cls.mask_dummy):
+            assert('objects: 12' == params_str[1])
+        assert('representation: ' + cls.anno_repr == params_str[3])
+
+    def test_count_query(self):
+        query_file = self.tempdir.name + '/query.fa'
+        queries = [
+            'AAA',
+            'AAAA',
+            'AAAAAAAAAAAAA',
+            'CCC',
+            'CCCC',
+            'CCCCCCCCCCCCC',
+            'TTT',
+            'AAACCCGGGTTT',
+            'AAACCCGGGTTTTTT',
+            'AAACCCGGGTTTAAA',
+            'TTTAAACCCGGG',
+            'ACACACACACACATTTAAACCCGGG',
+        ]
+        for discovery_rate in np.linspace(0, 1, 5):
+            expected_output = ''
+            with open(query_file, 'w') as f:
+                for i, s in enumerate(queries):
+                    f.write(f'>s{i}\n{s}\n')
+                    expected_output += f'{i}\ts{i}'
+                    def get_count(d, kmer):
+                        try:
+                            return d[kmer]
+                        except:
+                            return 0
+
+                    num_kmers = len(s) - self.k + 1
+
+                    num_matches_1 = sum([get_count(self.kmer_counts_1, s[i:i + self.k]) > 0 for i in range(num_kmers)])
+                    count_1 = sum([get_count(self.kmer_counts_1, s[i:i + self.k]) for i in range(len(s) - self.k + 1)])
+
+                    num_matches_2 = sum([get_count(self.kmer_counts_2, s[i:i + self.k]) > 0 for i in range(num_kmers)])
+                    count_2 = sum([get_count(self.kmer_counts_2, s[i:i + self.k]) for i in range(len(s) - self.k + 1)])
+
+                    for (c, i, n) in [(count_1, 1, num_matches_1), (count_2, 0, num_matches_2)]:
+                        if n >= discovery_rate * num_kmers:
+                            expected_output += f'\t<L{2-i}>:{c}'
+
+                    expected_output += '\n'
+
+            query_command = f'{METAGRAPH} query --fast --count-kmers \
+                            -i {self.tempdir.name}/graph{graph_file_extension[self.graph_repr]} \
+                            -a {self.tempdir.name}/annotation{anno_file_extension[self.anno_repr]} \
+                            --discovery-fraction {discovery_rate} {query_file}'
+
+            res = subprocess.run(query_command.split(), stdout=PIPE)
+            self.assertEqual(res.returncode, 0)
+            self.assertEqual(res.stdout.decode(), expected_output)
+
+    def test_count_quantiles(self):
+        query_file = self.tempdir.name + '/query.fa'
+        queries = [
+            'AAA',
+            'AAAA',
+            'AAAAAAAAAAAAA',
+            'CCC',
+            'CCCC',
+            'CCCCCCCCCCCCC',
+            'TTT',
+            'AAACCCGGGTTT',
+            'AAACCCGGGTTTTTT',
+            'AAACCCGGGTTTAAA',
+            'TTTAAACCCGGG',
+            'ACACACACACACATTTAAACCCGGG',
+        ]
+        quantiles = np.linspace(0, 1, 100)
+        expected_output = ''
+        with open(query_file, 'w') as f:
+            for i, s in enumerate(queries):
+                f.write(f'>s{i}\n{s}\n')
+                expected_output += f'{i}\ts{i}\t<L1>'
+                def get_count(d, kmer):
+                    try:
+                        return d[kmer]
+                    except:
+                        return 0
+                for p in quantiles:
+                    counts = [get_count(self.kmer_counts_1, s[i:i + self.k]) for i in range(len(s) - self.k + 1)]
+                    expected_output += f':{np.quantile(counts, p, interpolation="lower")}'
+                expected_output += f'\t<L2>'
+                for p in quantiles:
+                    counts = [get_count(self.kmer_counts_2, s[i:i + self.k]) for i in range(len(s) - self.k + 1)]
+                    expected_output += f':{np.quantile(counts, p, interpolation="lower")}'
+                expected_output += '\n'
+
+        query_command = f'{METAGRAPH} query --fast --count-quantiles <SET_BELOW> \
+                        -i {self.tempdir.name}/graph{graph_file_extension[self.graph_repr]} \
+                        -a {self.tempdir.name}/annotation{anno_file_extension[self.anno_repr]} \
+                        --discovery-fraction 0.0 {query_file}'
+
+        query_command = query_command.split()
+        query_command[4] = ' '.join([str(p) for p in quantiles])
+        res = subprocess.run(query_command, stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        self.assertEqual(res.stdout.decode(), expected_output)
+
+        query_command = f'{METAGRAPH} query --fast --count-quantiles <SET_BELOW> \
+                        -i {self.tempdir.name}/graph{graph_file_extension[self.graph_repr]} \
+                        -a {self.tempdir.name}/annotation{anno_file_extension[self.anno_repr]} \
+                        --discovery-fraction 1.0 {query_file}'
+
+        query_command = query_command.split()
+        query_command[4] = ' '.join([str(p) for p in quantiles])
+        res = subprocess.run(query_command, stdout=PIPE)
+        self.assertEqual(res.returncode, 0)
+        self.assertEqual(len(res.stdout.decode()), 5230)
+
 
 @parameterized_class(('graph_repr', 'anno_repr'),
     input_values=(product(list(set(GRAPH_TYPES) - {'hashstr'}), ANNO_TYPES) +
diff --git a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp
index 3942571a39..c9174e188d 100644
--- a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp
+++ b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp
@@ -56,6 +56,39 @@ ColumnMajor::get_rows(const std::vector<Row> &row_ids) const {
     return rows;
 }
 
+Vector<std::pair<ColumnMajor::Column, uint64_t>>
+ColumnMajor::get_column_ranks(Row row) const {
+    assert(row < num_rows() || !columns_.size());
+
+    Vector<std::pair<Column, uint64_t>> result;
+    for (size_t i = 0; i < columns_.size(); ++i) {
+        assert(columns_[i]);
+
+        if (uint64_t r = columns_[i]->conditional_rank1(row))
+            result.emplace_back(i, r);
+    }
+    return result;
+}
+
+std::vector<Vector<std::pair<ColumnMajor::Column, uint64_t>>>
+ColumnMajor::get_column_ranks(const std::vector<Row> &row_ids) const {
+    std::vector<Vector<std::pair<Column, uint64_t>>> result(row_ids.size());
+
+    for (size_t j = 0; j < columns_.size(); ++j) {
+        assert(columns_[j]);
+        const bit_vector &col = *columns_[j];
+
+        for (size_t i = 0; i < row_ids.size(); ++i) {
+            assert(row_ids[i] < num_rows());
+
+            if (uint64_t r = col.conditional_rank1(row_ids[i]))
+                result[i].emplace_back(j, r);
+        }
+    }
+
+    return result;
+}
+
 std::vector<ColumnMajor::Column>
 ColumnMajor::slice_rows(const std::vector<Row> &row_ids) const {
     std::vector<Column> slice;
@@ -91,6 +124,7 @@ bool ColumnMajor::load(std::istream &in) {
         for (auto &c : columns_) {
             assert(!c);
 
+            // TODO: switch to bit_vector_smart?
             c = std::make_unique<bit_vector_sd>();
             if (!c->load(in))
                 return false;
diff --git a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp
index 91d7451ec0..7aa95bf639 100644
--- a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp
+++ b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp
@@ -23,6 +23,10 @@ class ColumnMajor : public BinaryMatrix {
     bool get(Row row, Column column) const override;
     SetBitPositions get_row(Row row) const override;
     std::vector<SetBitPositions> get_rows(const std::vector<Row> &rows) const override;
+    // query row and get ranks of each set bit in its column
+    Vector<std::pair<Column, uint64_t>> get_column_ranks(Row row) const;
+    std::vector<Vector<std::pair<Column, uint64_t>>>
+    get_column_ranks(const std::vector<Row> &rows) const;
     std::vector<Row> get_column(Column column) const override;
     // get all selected rows appended with -1 and concatenated
     std::vector<Column> slice_rows(const std::vector<Row> &rows) const override;
diff --git a/metagraph/src/annotation/int_matrix/base/int_matrix.cpp b/metagraph/src/annotation/int_matrix/base/int_matrix.cpp
index 07562f983a..11b8418a2f 100644
--- a/metagraph/src/annotation/int_matrix/base/int_matrix.cpp
+++ b/metagraph/src/annotation/int_matrix/base/int_matrix.cpp
@@ -7,29 +7,31 @@ namespace matrix {
 
 IntMatrix::RowValues
 IntMatrix::sum_row_values(const std::vector<std::pair<Row, size_t>> &index_counts,
-                          size_t min,
-                          size_t cap) const {
-    assert(cap >= min);
-
-    if (!cap)
-        return {};
-
-    min = std::max(min, size_t(1));
+                          size_t min_count) const {
+    min_count = std::max(min_count, size_t(1));
 
+    std::vector<Row> rows;
+    rows.reserve(index_counts.size());
     size_t total_sum = 0;
-    for (const auto &pair : index_counts) {
-        total_sum += pair.second;
+    for (const auto &[i, count] : index_counts) {
+        total_sum += count;
+        rows.push_back(i);
     }
 
-    if (total_sum < min)
+    if (total_sum < min_count)
         return {};
 
     std::vector<size_t> sum_row(num_columns(), 0);
+    std::vector<size_t> counts(num_columns(), 0);
 
-    for (auto [i, count] : index_counts) {
-        for (auto [j, value] : get_row_values(i)) {
+    auto row_values = get_row_values(rows);
+
+    for (size_t t = 0; t < index_counts.size(); ++t) {
+        auto [i, count] = index_counts[t];
+        for (const auto &[j, value] : row_values[t]) {
             assert(j < sum_row.size());
             sum_row[j] += count * value;
+            counts[j] += count;
         }
     }
 
@@ -37,14 +39,47 @@ IntMatrix::sum_row_values(const std::vector<std::pair<Row, size_t>> &index_count
     result.reserve(sum_row.size());
 
     for (size_t j = 0; j < num_columns(); ++j) {
-        if (sum_row[j] >= min) {
-            result.emplace_back(j, std::min(sum_row[j], cap));
+        if (counts[j] >= min_count) {
+            result.emplace_back(j, sum_row[j]);
         }
     }
 
     return result;
 }
 
+
+// return sizes of all non-empty tuples in the row
+MultiIntMatrix::RowValues MultiIntMatrix::get_row_values(Row row) const {
+    RowTuples row_tuples = get_row_tuples(row);
+
+    RowValues row_values(row_tuples.size());
+
+    for (size_t i = 0; i < row_tuples.size(); ++i) {
+        row_values[i].first = row_tuples[i].first;
+        row_values[i].second = row_tuples[i].second.size();
+    }
+
+    return row_values;
+}
+
+// for each row return the sizes of all non-empty tuples
+std::vector<MultiIntMatrix::RowValues>
+MultiIntMatrix::get_row_values(const std::vector<Row> &rows) const {
+    std::vector<RowTuples> row_tuples = get_row_tuples(rows);
+
+    std::vector<RowValues> row_values(row_tuples.size());
+
+    for (size_t i = 0; i < row_tuples.size(); ++i) {
+        row_values[i].resize(row_tuples[i].size());
+        for (size_t j = 0; j < row_tuples[i].size(); ++j) {
+            row_values[i][j].first = row_tuples[i][j].first;
+            row_values[i][j].second = row_tuples[i][j].second.size();
+        }
+    }
+
+    return row_values;
+}
+
 } // namespace matrix
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/annotation/int_matrix/base/int_matrix.hpp b/metagraph/src/annotation/int_matrix/base/int_matrix.hpp
index 1f595b5082..00beb51f6e 100644
--- a/metagraph/src/annotation/int_matrix/base/int_matrix.hpp
+++ b/metagraph/src/annotation/int_matrix/base/int_matrix.hpp
@@ -23,13 +23,38 @@ class IntMatrix : public binmat::BinaryMatrix {
     virtual std::vector<RowValues>
     get_row_values(const std::vector<Row> &rows) const = 0;
 
-    // Get all columns for which the sum of the values in queried rows
-    // is greater than or equal to |min|. Stop counting if the sum is
-    // greater than |cap|.
+    // sum up values for each column with at least |min_count| non-zero values
     virtual RowValues
     sum_row_values(const std::vector<std::pair<Row, size_t>> &index_counts,
-                   size_t min = 1,
-                   size_t cap = std::numeric_limits<size_t>::max()) const;
+                   size_t min_count = 1) const;
+};
+
+
+// Entries are tuples and their aggregated `values` are tuple sizes
+class MultiIntMatrix : public IntMatrix {
+  public:
+    typedef SmallVector<uint64_t> Tuple;
+    typedef Vector<std::pair<Column, Tuple>> RowTuples;
+
+    virtual ~MultiIntMatrix() {}
+
+    // return tuple sizes (if not zero) at each entry
+    virtual RowValues get_row_values(Row row) const;
+
+    virtual std::vector<RowValues>
+    get_row_values(const std::vector<Row> &rows) const;
+
+    // return total number of attributes in all tuples
+    virtual uint64_t num_attributes() const = 0;
+
+    // return entries of the matrix -- where each entry is a set of integers
+    virtual RowTuples get_row_tuples(Row row) const = 0;
+
+    virtual std::vector<RowTuples>
+    get_row_tuples(const std::vector<Row> &rows) const = 0;
+
+    virtual bool load_tuples(std::istream &in) = 0;
+    virtual void serialize_tuples(std::ostream &out) const = 0;
 };
 
 } // namespace matrix
diff --git a/metagraph/src/annotation/int_matrix/rank_extended/tuple_csc_matrix.hpp b/metagraph/src/annotation/int_matrix/rank_extended/tuple_csc_matrix.hpp
new file mode 100644
index 0000000000..08c5a78b11
--- /dev/null
+++ b/metagraph/src/annotation/int_matrix/rank_extended/tuple_csc_matrix.hpp
@@ -0,0 +1,203 @@
+#ifndef __TUPLE_CSC_MATRIX_HPP__
+#define __TUPLE_CSC_MATRIX_HPP__
+
+#include <vector>
+
+#include "annotation/int_matrix/base/int_matrix.hpp"
+#include "common/logger.hpp"
+
+
+namespace mtg {
+namespace annot {
+namespace matrix {
+
+/**
+ * Multi-Value Compressed Sparse Column Matrix (column-rank extended)
+ *
+ * Matrix which stores the non-empty tuples externally and indexes their
+ * positions in a binary matrix. These values are indexed by rank1 called
+ * on binary columns of the indexing matrix.
+ */
+template <class BaseMatrix,
+          class Values = sdsl::int_vector<>,
+          class Delims = bit_vector_smart>
+class TupleCSCMatrix : public MultiIntMatrix {
+  public:
+    TupleCSCMatrix() {}
+
+    TupleCSCMatrix(BaseMatrix&& index_matrix)
+      : binary_matrix_(std::move(index_matrix)) {}
+
+    // return tuple sizes (if not zero) at each entry
+    RowValues get_row_values(Row row) const;
+
+    std::vector<RowValues>
+    get_row_values(const std::vector<Row> &rows) const;
+
+    uint64_t num_attributes() const;
+
+    // return entries of the matrix -- where each entry is a set of integers
+    RowTuples get_row_tuples(Row row) const;
+
+    std::vector<RowTuples>
+    get_row_tuples(const std::vector<Row> &rows) const;
+
+    uint64_t num_columns() const { return binary_matrix_.num_columns(); }
+    uint64_t num_rows() const { return binary_matrix_.num_rows(); }
+    uint64_t num_relations() const { return binary_matrix_.num_relations(); }
+
+    // row is in [0, num_rows), column is in [0, num_columns)
+    bool get(Row row, Column column) const { return binary_matrix_.get(row, column); }
+    SetBitPositions get_row(Row row) const { return binary_matrix_.get_row(row); }
+    std::vector<SetBitPositions> get_rows(const std::vector<Row> &rows) const {
+        return binary_matrix_.get_rows(rows);
+    }
+    std::vector<Row> get_column(Column column) const {
+        return binary_matrix_.get_column(column);
+    }
+    // get all selected rows appended with -1 and concatenated
+    std::vector<Column> slice_rows(const std::vector<Row> &rows) const {
+        return binary_matrix_.slice_rows(rows);
+    }
+
+    bool load(std::istream &in);
+    void serialize(std::ostream &out) const;
+
+    bool load_tuples(std::istream &in);
+    void serialize_tuples(std::ostream &out) const;
+
+    const BaseMatrix& get_binary_matrix() const { return binary_matrix_; }
+
+  private:
+    BaseMatrix binary_matrix_;
+    std::vector<Delims> delimiters_;
+    std::vector<Values> column_values_;
+};
+
+
+template <class BaseMatrix, class Values, class Delims>
+inline typename TupleCSCMatrix<BaseMatrix, Values, Delims>::RowValues
+TupleCSCMatrix<BaseMatrix, Values, Delims>::get_row_values(Row row) const {
+    const auto &column_ranks = binary_matrix_.get_column_ranks(row);
+    RowValues row_values;
+    row_values.reserve(column_ranks.size());
+    for (auto [j, r] : column_ranks) {
+        assert(r >= 1 && "matches can't have zero-rank");
+        size_t tuple_size = delimiters_[j].select1(r + 1) - delimiters_[j].select1(r) - 1;
+        row_values.emplace_back(j, tuple_size);
+    }
+    return row_values;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline std::vector<typename TupleCSCMatrix<BaseMatrix, Values, Delims>::RowValues>
+TupleCSCMatrix<BaseMatrix, Values, Delims>::get_row_values(const std::vector<Row> &rows) const {
+    const auto &column_ranks = binary_matrix_.get_column_ranks(rows);
+    std::vector<RowValues> row_values(rows.size());
+    // TODO: reshape?
+    for (size_t i = 0; i < rows.size(); ++i) {
+        row_values[i].reserve(column_ranks[i].size());
+        for (auto [j, r] : column_ranks[i]) {
+            assert(r >= 1 && "matches can't have zero-rank");
+            size_t tuple_size = delimiters_[j].select1(r + 1) - delimiters_[j].select1(r) - 1;
+            row_values[i].emplace_back(j, tuple_size);
+        }
+    }
+    return row_values;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+uint64_t TupleCSCMatrix<BaseMatrix, Values, Delims>::num_attributes() const {
+    uint64_t num_attributes = 0;
+    for (size_t j = 0; j < column_values_.size(); ++j) {
+        num_attributes += column_values_[j].size();
+    }
+    return num_attributes;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline typename TupleCSCMatrix<BaseMatrix, Values, Delims>::RowTuples
+TupleCSCMatrix<BaseMatrix, Values, Delims>::get_row_tuples(Row row) const {
+    const auto &column_ranks = binary_matrix_.get_column_ranks(row);
+    RowTuples row_tuples;
+    row_tuples.reserve(column_ranks.size());
+    for (auto [j, r] : column_ranks) {
+        assert(r >= 1 && "matches can't have zero-rank");
+        size_t begin = delimiters_[j].select1(r) + 1 - r;
+        size_t end = delimiters_[j].select1(r + 1) - r;
+        Tuple tuple;
+        tuple.reserve(end - begin);
+        for (size_t t = begin; t < end; ++t) {
+            tuple.push_back(column_values_[j][t]);
+        }
+        row_tuples.emplace_back(j, std::move(tuple));
+    }
+    return row_tuples;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline std::vector<typename TupleCSCMatrix<BaseMatrix, Values, Delims>::RowTuples>
+TupleCSCMatrix<BaseMatrix, Values, Delims>::get_row_tuples(const std::vector<Row> &rows) const {
+    const auto &column_ranks = binary_matrix_.get_column_ranks(rows);
+    std::vector<RowTuples> row_tuples(rows.size());
+    // TODO: reshape?
+    for (size_t i = 0; i < rows.size(); ++i) {
+        row_tuples[i].reserve(column_ranks[i].size());
+        for (auto [j, r] : column_ranks[i]) {
+            assert(r >= 1 && "matches can't have zero-rank");
+            size_t begin = delimiters_[j].select1(r) + 1 - r;
+            size_t end = delimiters_[j].select1(r + 1) - r;
+            Tuple tuple;
+            tuple.reserve(end - begin);
+            for (size_t t = begin; t < end; ++t) {
+                tuple.push_back(column_values_[j][t]);
+            }
+            row_tuples[i].emplace_back(j, std::move(tuple));
+        }
+    }
+    return row_tuples;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline bool TupleCSCMatrix<BaseMatrix, Values, Delims>::load(std::istream &in) {
+    return binary_matrix_.load(in) && load_tuples(in);
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline bool TupleCSCMatrix<BaseMatrix, Values, Delims>::load_tuples(std::istream &in) {
+    delimiters_.clear();
+    column_values_.clear();
+
+    delimiters_.resize(num_columns());
+    column_values_.resize(num_columns());
+    for (size_t j = 0; j < column_values_.size(); ++j) {
+        try {
+            delimiters_[j].load(in);
+            column_values_[j].load(in);
+        } catch (...) {
+            common::logger->error("Couldn't load tuple attributes for column {}", j);
+            return false;
+        }
+    }
+    return true;
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline void TupleCSCMatrix<BaseMatrix, Values, Delims>::serialize(std::ostream &out) const {
+    binary_matrix_.serialize(out);
+    serialize_tuples(out);
+}
+
+template <class BaseMatrix, class Values, class Delims>
+inline void TupleCSCMatrix<BaseMatrix, Values, Delims>::serialize_tuples(std::ostream &out) const {
+    for (size_t j = 0; j < column_values_.size(); ++j) {
+        delimiters_[j].serialize(out);
+        column_values_[j].serialize(out);
+    }
+}
+
+} // namespace matrix
+} // namespace annot
+} // namespace mtg
+
+#endif // __TUPLE_CSC_MATRIX_HPP__
diff --git a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
index 211e5722a9..32f06e4909 100644
--- a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
+++ b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
@@ -214,5 +214,7 @@ template class StaticBinRelAnnotator<matrix::IntRowDiff<matrix::CSCMatrix<binmat
 
 template class StaticBinRelAnnotator<matrix::CSRMatrix, std::string>;
 
+template class StaticBinRelAnnotator<matrix::TupleCSCMatrix<binmat::ColumnMajor>, std::string>;
+
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
index 889afa4ca2..d1f53f31ac 100644
--- a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
+++ b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
@@ -16,6 +16,7 @@
 #include "annotation/int_matrix/rank_extended/csc_matrix.hpp"
 #include "annotation/int_matrix/row_diff/int_row_diff.hpp"
 #include "annotation/int_matrix/csr_matrix/csr_matrix.hpp"
+#include "annotation/int_matrix/rank_extended/tuple_csc_matrix.hpp"
 
 
 namespace mtg {
@@ -51,6 +52,8 @@ typedef StaticBinRelAnnotator<matrix::IntRowDiff<matrix::CSCMatrix<binmat::BRWT,
 
 typedef StaticBinRelAnnotator<matrix::CSRMatrix, std::string> IntRowAnnotator;
 
+typedef StaticBinRelAnnotator<matrix::TupleCSCMatrix<binmat::ColumnMajor>, std::string> ColumnCoordAnnotator;
+
 
 template <>
 inline const std::string RowFlatAnnotator::kExtension = ".flat.annodbg";
@@ -80,6 +83,8 @@ template <>
 inline const std::string IntRowDiffBRWTAnnotator::kExtension = ".row_diff_int_brwt.annodbg";
 template <>
 inline const std::string IntRowAnnotator::kExtension = ".int_csr.annodbg";
+template <>
+inline const std::string ColumnCoordAnnotator::kExtension = ".column_coord.annodbg";
 
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
index 9a5f419262..8cdf19f6ce 100644
--- a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
+++ b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
@@ -766,10 +766,11 @@ bitmap_builder& ColumnCompressed<Label>::decompress_builder(size_t j) {
 
     flushed_ = false;
 
-    try {
+    if (auto cached = cached_columns_.TryGet(j)) {
         // check the  the cached bitmap builder
-        return *cached_columns_.Get(j);
-    } catch (...) {
+        return **cached;
+
+    } else {
         assert(j <= bitmatrix_.size());
 
         bitmap_builder *vector;
@@ -800,7 +801,7 @@ bitmap_builder& ColumnCompressed<Label>::decompress_builder(size_t j) {
         }
 
         cached_columns_.Put(j, vector);
-        return *cached_columns_.Get(j);
+        return *vector;
     }
 }
 
diff --git a/metagraph/src/cli/align.cpp b/metagraph/src/cli/align.cpp
index 8e2d429f9b..ef25a15b0d 100644
--- a/metagraph/src/cli/align.cpp
+++ b/metagraph/src/cli/align.cpp
@@ -211,15 +211,14 @@ void map_sequences_in_file(const std::string &file,
 
         if (config.count_kmers) {
             std::sort(graphindices.begin(), graphindices.end());
-            size_t num_unique_matching_kmers = std::inner_product(
-                graphindices.begin() + 1, graphindices.end(),
-                graphindices.begin(),
-                size_t(graphindices.front() != DeBruijnGraph::npos),
-                std::plus<size_t>(),
-                [](DeBruijnGraph::node_index next, DeBruijnGraph::node_index prev) {
-                    return next != DeBruijnGraph::npos && next != prev;
-                }
-            );
+            size_t num_unique_matching_kmers = 0;
+            auto prev = DeBruijnGraph::npos;
+            for (auto i : graphindices) {
+                if (i != DeBruijnGraph::npos && i != prev)
+                    ++num_unique_matching_kmers;
+
+                prev = i;
+            }
             *out << read_stream->name.s << "\t"
                  << num_discovered << "/" << num_kmers << "/"
                  << num_unique_matching_kmers << "\n";
@@ -236,6 +235,7 @@ void map_sequences_in_file(const std::string &file,
             // map input subsequences to multiple nodes
             for (size_t i = 0; i + graph.get_k() <= read_stream->seq.l; ++i) {
                 // TODO: make more efficient
+                // TODO: canonicalization
                 std::string_view subseq(read_stream->seq.s + i, config.alignment_length);
 
                 dbg->call_nodes_with_suffix_matching_longest_prefix(
@@ -463,10 +463,10 @@ int align_to_graph(Config *config) {
                 num_bytes_read += it->seq.l;
             }
 
-            auto process_batch = [&](SeqBatch batch, uint64_t size) {
+            auto process_batch = [&](SeqBatch batch) {
                 auto aln_graph = graph;
-                if (graph->get_mode() == DeBruijnGraph::PRIMARY)
-                    aln_graph = std::make_shared<CanonicalDBG>(aln_graph, size);
+                if (auto *canonical = dynamic_cast<CanonicalDBG*>(graph.get()))
+                    aln_graph = std::make_shared<CanonicalDBG>(*canonical);
 
                 auto aligner = build_aligner(*aln_graph, aligner_config);
 
@@ -500,15 +500,14 @@ int align_to_graph(Config *config) {
                     }
 
                     thread_pool.enqueue(process_batch,
-                                        SeqBatch(last_mv_it, std::make_move_iterator(it)),
-                                        mbatch_size);
+                                        SeqBatch(last_mv_it, std::make_move_iterator(it)));
                     ++num_minibatches;
                 }
 
                 logger->trace("Num minibatches: {}, minibatch size: {} KB",
                               num_minibatches, mbatch_size / 1e3);
             } else {
-                thread_pool.enqueue(process_batch, std::move(seq_batch), batch_size);
+                thread_pool.enqueue(process_batch, std::move(seq_batch));
             }
         };
 
diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp
index a3f45a8717..d5c2536db4 100644
--- a/metagraph/src/cli/config/config.cpp
+++ b/metagraph/src/cli/config/config.cpp
@@ -180,6 +180,10 @@ Config::Config(int argc, char *argv[]) {
             for (const auto &border : utils::split_string(get_value(i++), " ")) {
                 count_slice_quantiles.push_back(std::stod(border));
             }
+        } else if (!strcmp(argv[i], "--count-quantiles")) {
+            for (const auto &p : utils::split_string(get_value(i++), " ")) {
+                count_quantiles.push_back(std::stod(p));
+            }
         } else if (!strcmp(argv[i], "--aggregate-columns")) {
             aggregate_columns = true;
         } else if (!strcmp(argv[i], "--intersected-anno")) {
@@ -196,6 +200,8 @@ Config::Config(int argc, char *argv[]) {
             discovery_fraction = std::stof(get_value(i++));
         } else if (!strcmp(argv[i], "--query-presence")) {
             query_presence = true;
+        } else if (!strcmp(argv[i], "--query-coords")) {
+            query_coords = true;
         } else if (!strcmp(argv[i], "--filter-present")) {
             filter_present = true;
         } else if (!strcmp(argv[i], "--count-labels")) {
@@ -696,6 +702,8 @@ std::string Config::annotype_to_string(AnnotationType state) {
             return "int_brwt";
         case IntRowDiffBRWT:
             return "row_diff_int_brwt";
+        case ColumnCoord:
+            return "column_coord";
     }
     throw std::runtime_error("Never happens");
 }
@@ -729,6 +737,8 @@ Config::AnnotationType Config::string_to_annotype(const std::string &string) {
         return AnnotationType::IntBRWT;
     } else if (string == "row_diff_int_brwt") {
         return AnnotationType::IntRowDiffBRWT;
+    } else if (string == "column_coord") {
+        return AnnotationType::ColumnCoord;
     } else {
         std::cerr << "Error: unknown annotation representation" << std::endl;
         exit(1);
@@ -791,6 +801,7 @@ DeBruijnGraph::Mode Config::string_to_graphmode(const std::string &string) {
 
 void Config::print_usage(const std::string &prog_name, IdentityType identity) {
     const char annotation_list[] = "\t\t( column, brwt, rb_brwt, int_brwt,\n"
+                                   "\t\t  column_coord,\n"
                                    "\t\t  row_diff, row_diff_brwt, row_diff_sparse, row_diff_int_brwt,\n"
                                    "\t\t  row, flat, row_sparse, rbfish, bin_rel_wt, bin_rel_wt_sdsl )";
 
@@ -1168,6 +1179,10 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) {
             fprintf(stderr, "\n");
             fprintf(stderr, "\t   --count-labels \t\tcount labels for k-mers from querying sequences [off]\n");
             fprintf(stderr, "\t   --count-kmers \t\tweight k-mers with their annotated counts (requires count annotation) [off]\n");
+            fprintf(stderr, "\t   --count-quantiles [FLOAT ...] \tk-mer count quantiles to compute for each label [off]\n"
+                            "\t                                 \t\tExample: --count-quantiles '0 0.33 0.5 0.66 1'\n"
+                            "\t                                 \t\t(0 corresponds to MIN, 1 corresponds to MAX)\n");
+            fprintf(stderr, "\t   --query-coords \t\tquery k-mer coordinates (requires coord annotation) [off]\n");
             fprintf(stderr, "\t   --print-signature \t\tprint vectors indicating present/absent k-mers [off]\n");
             fprintf(stderr, "\t   --num-top-labels \t\tmaximum number of frequent labels to print [off]\n");
             fprintf(stderr, "\t   --discovery-fraction [FLOAT] fraction of labeled k-mers required for annotation [0.7]\n");
@@ -1186,7 +1201,7 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) {
             fprintf(stderr, "\t   --align-min-path-score [INT]\t\t\tthe minimum score that a reported path can have [0]\n");
             fprintf(stderr, "\t   --align-edit-distance \t\t\tuse unit costs for scoring matrix [off]\n");
             fprintf(stderr, "\t   --align-max-nodes-per-seq-char [FLOAT]\tmaximum number of nodes to consider per sequence character [12.0]\n");
-            fprintf(stderr, "\t   --align-max-ram [FLOAT]\t\tmaximum amount of RAM used per alignment in MB [200.0]\n");
+            fprintf(stderr, "\t   --align-max-ram [FLOAT]\t\t\tmaximum amount of RAM used per alignment in MB [200.0]\n");
             fprintf(stderr, "\n");
             fprintf(stderr, "\t   --batch-align \t\talign against query graph [off]\n");
             fprintf(stderr, "\t   --max-hull-forks [INT]\tmaximum number of forks to take when expanding query graph [4]\n");
diff --git a/metagraph/src/cli/config/config.hpp b/metagraph/src/cli/config/config.hpp
index 5649a2ef18..4a19a3ab32 100644
--- a/metagraph/src/cli/config/config.hpp
+++ b/metagraph/src/cli/config/config.hpp
@@ -39,6 +39,7 @@ class Config {
     bool count_kmers = false;
     bool print_signature = false;
     bool query_presence = false;
+    bool query_coords = false;
     bool filter_present = false;
     bool dump_text_anno = false;
     bool sparse = false;
@@ -133,6 +134,7 @@ class Config {
     double min_fraction = 0.0;
     double max_fraction = 1.0;
     std::vector<double> count_slice_quantiles;
+    std::vector<double> count_quantiles;
 
     std::vector<std::string> fnames;
     std::vector<std::string> anno_labels;
@@ -198,6 +200,7 @@ class Config {
         RbBRWT,
         IntBRWT,
         IntRowDiffBRWT,
+        ColumnCoord,
     };
 
     enum GraphType {
diff --git a/metagraph/src/cli/load/load_annotated_graph.cpp b/metagraph/src/cli/load/load_annotated_graph.cpp
index 68b1977930..3c9e0ec90c 100644
--- a/metagraph/src/cli/load/load_annotated_graph.cpp
+++ b/metagraph/src/cli/load/load_annotated_graph.cpp
@@ -17,7 +17,6 @@ namespace mtg {
 namespace cli {
 
 using namespace mtg::graph;
-
 using mtg::common::logger;
 
 
@@ -52,6 +51,7 @@ std::unique_ptr<AnnotatedDBG> initialize_annotated_dbg(std::shared_ptr<DeBruijnG
                           config.infbase);
             exit(1);
         }
+
         // row_diff annotation is special, as it must know the graph structure
         using namespace annot::binmat;
         BinaryMatrix &matrix = const_cast<BinaryMatrix &>(annotation_temp->get_matrix());
diff --git a/metagraph/src/cli/load/load_annotation.cpp b/metagraph/src/cli/load/load_annotation.cpp
index 9482b47fd8..d9865ad60a 100644
--- a/metagraph/src/cli/load/load_annotation.cpp
+++ b/metagraph/src/cli/load/load_annotation.cpp
@@ -21,6 +21,9 @@ Config::AnnotationType parse_annotation_type(const std::string &filename) {
     if (utils::ends_with(filename, annot::ColumnCompressed<>::kExtension)) {
         return Config::AnnotationType::ColumnCompressed;
 
+    } else if (utils::ends_with(filename, annot::ColumnCoordAnnotator::kExtension)) {
+        return Config::AnnotationType::ColumnCoord;
+
     } else if (utils::ends_with(filename, annot::RowDiffColumnAnnotator::kExtension)) {
         return Config::AnnotationType::RowDiff;
 
@@ -137,6 +140,10 @@ initialize_annotation(Config::AnnotationType anno_type,
             annotation.reset(new annot::IntRowDiffBRWTAnnotator());
             break;
         }
+        case Config::ColumnCoord: {
+            annotation.reset(new annot::ColumnCoordAnnotator());
+            break;
+        }
     }
 
     return annotation;
diff --git a/metagraph/src/cli/query.cpp b/metagraph/src/cli/query.cpp
index 41d8631d17..13322d0f7d 100644
--- a/metagraph/src/cli/query.cpp
+++ b/metagraph/src/cli/query.cpp
@@ -60,7 +60,9 @@ std::string QueryExecutor::execute_query(const std::string &seq_name,
                                          double discovery_fraction,
                                          std::string anno_labels_delimiter,
                                          const AnnotatedDBG &anno_graph,
-                                         bool with_kmer_counts) {
+                                         bool with_kmer_counts,
+                                         const std::vector<double> &count_quantiles,
+                                         bool query_coords) {
     std::string output;
     output.reserve(1'000);
 
@@ -84,7 +86,43 @@ std::string QueryExecutor::execute_query(const std::string &seq_name,
 
         output += '\n';
 
-    } else if (count_labels) {
+    } else if (query_coords) {
+        auto result = anno_graph.get_kmer_coordinates(sequence,
+                                                      num_top_labels,
+                                                      discovery_fraction);
+
+        if (!result.size() && suppress_unlabeled)
+            return "";
+
+        output += seq_name;
+
+        for (const auto &[label, tuples] : result) {
+            output += "\t<" + label + ">";
+            for (const auto &coords : tuples) {
+                output += fmt::format(":{}", fmt::join(coords, ","));
+            }
+        }
+
+        output += '\n';
+
+    } else if (count_quantiles.size()) {
+        auto result = anno_graph.get_label_count_quantiles(sequence,
+                                                           num_top_labels,
+                                                           discovery_fraction,
+                                                           count_quantiles);
+
+        if (!result.size() && suppress_unlabeled)
+            return "";
+
+        output += seq_name;
+
+        for (const auto &[label, quantiles] : result) {
+            output += fmt::format("\t<{}>:{}", label, fmt::join(quantiles, ":"));
+        }
+
+        output += '\n';
+
+    } else if (count_labels || with_kmer_counts) {
         auto top_labels = anno_graph.get_top_labels(sequence,
                                                     num_top_labels,
                                                     discovery_fraction,
@@ -96,10 +134,7 @@ std::string QueryExecutor::execute_query(const std::string &seq_name,
         output += seq_name;
 
         for (const auto &[label, count] : top_labels) {
-            output += "\t<";
-            output += label;
-            output += ">:";
-            output += fmt::format_int(count).c_str();
+            output += fmt::format("\t<{}>:{}", label, count);
         }
 
         output += '\n';
@@ -893,7 +928,8 @@ std::string query_sequence(size_t id, std::string name, std::string seq,
                                         config.count_labels, config.print_signature,
                                         config.suppress_unlabeled, config.num_top_labels,
                                         config.discovery_fraction, config.anno_labels_delimiter,
-                                        anno_graph, config.count_kmers);
+                                        anno_graph, config.count_kmers, config.count_quantiles,
+                                        config.query_coords);
 }
 
 void QueryExecutor::query_fasta(const string &file,
@@ -903,6 +939,10 @@ void QueryExecutor::query_fasta(const string &file,
     seq_io::FastaParser fasta_parser(file, config_.forward_and_reverse);
 
     if (config_.fast) {
+        if (config_.query_coords) {
+            logger->error("Querying coordinates in batch mode is not supported");
+            exit(1);
+        }
         // Construct a query graph and query against it
         batched_query_fasta(fasta_parser, callback);
         return;
diff --git a/metagraph/src/cli/query.hpp b/metagraph/src/cli/query.hpp
index c007d9bada..f39b16feaa 100644
--- a/metagraph/src/cli/query.hpp
+++ b/metagraph/src/cli/query.hpp
@@ -5,6 +5,7 @@
 #include <functional>
 #include <memory>
 #include <string>
+#include <vector>
 
 class ThreadPool;
 
@@ -63,7 +64,9 @@ class QueryExecutor {
                                      double discovery_fraction,
                                      std::string anno_labels_delimiter,
                                      const graph::AnnotatedDBG &anno_graph,
-                                     bool with_kmer_counts = false);
+                                     bool with_kmer_counts = false,
+                                     const std::vector<double> &count_quantiles = {},
+                                     bool query_coords = false);
 
   private:
     const Config &config_;
diff --git a/metagraph/src/cli/stats.cpp b/metagraph/src/cli/stats.cpp
index 36429282b3..c9961a49b1 100644
--- a/metagraph/src/cli/stats.cpp
+++ b/metagraph/src/cli/stats.cpp
@@ -167,6 +167,12 @@ void print_stats(const Annotator &annotation) {
               << utils::split_string(annotation.file_extension(), ".").at(0) << std::endl;
 
     using namespace annot::binmat;
+    using mtg::annot::matrix::MultiIntMatrix;
+
+    if (const auto *mat_coord = dynamic_cast<const MultiIntMatrix *>(&annotation.get_matrix())) {
+        std::cout << "================== COORDINATES STATS ===================" << std::endl;
+        std::cout << "coordinates: " << mat_coord->num_attributes() << std::endl;
+    }
 
     if (const auto *rbmat = dynamic_cast<const RainbowMatrix *>(&annotation.get_matrix())) {
         std::cout << "================= RAINBOW MATRIX STATS =================" << std::endl;
diff --git a/metagraph/src/cli/transform_annotation.cpp b/metagraph/src/cli/transform_annotation.cpp
index 8eda812d5e..402cd744e3 100644
--- a/metagraph/src/cli/transform_annotation.cpp
+++ b/metagraph/src/cli/transform_annotation.cpp
@@ -652,6 +652,27 @@ int transform_annotation(Config *config) {
                 assert(false);
                 break;
             }
+            case Config::ColumnCoord: {
+                auto label_encoder = annotator->get_label_encoder();
+                auto tuple_matrix = std::make_unique<matrix::TupleCSCMatrix<binmat::ColumnMajor>>(
+                                            annotator->release_matrix());
+                if (files.size() > 1) {
+                    logger->error("Merging coordinates from multiple columns is not supported");
+                    exit(1);
+                }
+                auto coords_fname = utils::remove_suffix(files.at(0),
+                                                         ColumnCompressed<>::kExtension)
+                                        + ColumnCompressed<>::kCoordExtension;
+                std::ifstream in(coords_fname);
+                tuple_matrix->load_tuples(in);
+
+                ColumnCoordAnnotator column_coord(std::move(tuple_matrix), label_encoder);
+
+                logger->trace("Annotation converted in {} sec", timer.elapsed());
+                column_coord.serialize(config->outfbase);
+                logger->trace("Serialized to {}", config->outfbase);
+                break;
+            }
             case Config::RowDiffBRWT: {
                 logger->error("Convert to row_diff first, and then to row_diff_brwt");
                 return 0;
diff --git a/metagraph/src/graph/annotated_dbg.cpp b/metagraph/src/graph/annotated_dbg.cpp
index a9480dc2dc..1461a9baf2 100644
--- a/metagraph/src/graph/annotated_dbg.cpp
+++ b/metagraph/src/graph/annotated_dbg.cpp
@@ -15,12 +15,15 @@
 #include "common/aligned_vector.hpp"
 #include "common/vectors/vector_algorithm.hpp"
 #include "common/vector_map.hpp"
+#include "common/logger.hpp"
 
 
 namespace mtg {
 namespace graph {
 
+using mtg::common::logger;
 using mtg::annot::matrix::IntMatrix;
+using mtg::annot::matrix::MultiIntMatrix;
 
 typedef AnnotatedDBG::Label Label;
 typedef std::pair<Label, size_t> StringCountPair;
@@ -294,6 +297,207 @@ AnnotatedDBG::get_top_labels(std::string_view sequence,
     return top_labels;
 }
 
+std::vector<std::pair<std::string, std::vector<size_t>>>
+AnnotatedDBG::get_label_count_quantiles(std::string_view sequence,
+                                        size_t num_top_labels,
+                                        double presence_ratio,
+                                        const std::vector<double> &count_quantiles) const {
+    assert(presence_ratio >= 0.);
+    assert(presence_ratio <= 1.);
+    assert(check_compatibility());
+    if (!std::is_sorted(count_quantiles.begin(), count_quantiles.end()))
+        throw std::runtime_error("Quantiles must be sorted");
+    if (count_quantiles.at(0) < 0. || count_quantiles.back() > 1.)
+        throw std::runtime_error("Quantiles must be in range [0, 1]");
+
+    if (sequence.size() < dbg_.get_k())
+        return {};
+
+    std::vector<row_index> rows;
+    size_t num_kmers = sequence.size() - dbg_.get_k() + 1;
+    rows.reserve(num_kmers);
+
+    graph_->map_to_nodes(sequence, [&](node_index i) {
+        if (i > 0)
+            rows.push_back(graph_to_anno_index(i));
+    });
+
+    uint64_t min_count = std::max(1.0, std::ceil(presence_ratio * num_kmers));
+    if (rows.size() < min_count)
+        return {};
+
+    std::vector<size_t> q_low(count_quantiles.size());
+    for (size_t i = 0; i < count_quantiles.size(); ++i) {
+        q_low[i] = (num_kmers - 1) * count_quantiles[i];
+    }
+
+    const auto *int_matrix = dynamic_cast<const IntMatrix *>(&annotator_->get_matrix());
+    if (!int_matrix) {
+        logger->error("k-mer counts are not indexed in this annotator");
+        exit(1);
+    }
+
+    VectorMap<size_t, std::vector<uint64_t>> code_to_counts;
+    for (const auto &row_values : int_matrix->get_row_values(rows)) {
+        for (const auto &[column, count] : row_values) {
+            code_to_counts[column].push_back(count);
+        }
+    }
+
+    std::vector<std::pair<size_t, std::vector<uint64_t>>> code_counts;
+    code_counts.reserve(code_to_counts.size());
+    for (auto &[j, counts] : code_to_counts.values_container()) {
+        // filter by the number of matched k-mers
+        if (counts.size() >= min_count)
+            code_counts.emplace_back(j, std::move(counts));
+    }
+    // sort by the number of matched k-mers
+    std::sort(code_counts.begin(), code_counts.end(),
+              [](const auto &x, const auto &y) {
+                  return x.second.size() > y.second.size()
+                      || (x.second.size() == y.second.size() && x.first < y.first);
+              });
+    // keep only the first |num_top_labels| top labels
+    if (code_counts.size() > num_top_labels)
+        code_counts.resize(num_top_labels);
+
+    std::vector<std::pair<Label, std::vector<size_t>>> label_quantiles;
+    label_quantiles.reserve(code_counts.size());
+    // Quantiles are defined as `count[i]` where `i < q * N <= i + 1`
+    for (auto &[j, counts] : code_counts) {
+        std::sort(counts.begin(), counts.end());
+        const size_t num_zeros = num_kmers - counts.size();
+
+        label_quantiles.emplace_back(annotator_->get_label_encoder().decode(j),
+                                     std::vector<size_t>(q_low.size()));
+
+        std::vector<size_t> &quantiles = label_quantiles.back().second;
+        for (size_t q = 0; q < q_low.size(); ++q) {
+            if (q_low[q] < num_zeros) {
+                quantiles[q] = 0;
+            } else {
+                quantiles[q] = counts[q_low[q] - num_zeros];
+            }
+        }
+    }
+
+    return label_quantiles;
+}
+
+std::vector<std::pair<std::string, std::vector<SmallVector<uint64_t>>>>
+AnnotatedDBG::get_kmer_coordinates(std::string_view sequence,
+                                   size_t num_top_labels,
+                                   double presence_ratio) const {
+    assert(presence_ratio >= 0.);
+    assert(presence_ratio <= 1.);
+    assert(check_compatibility());
+
+    if (sequence.size() < dbg_.get_k())
+        return {};
+
+    std::vector<node_index> path;
+    size_t num_kmers = sequence.size() - dbg_.get_k() + 1;
+    path.reserve(num_kmers);
+
+    graph_->map_to_nodes(sequence, [&](node_index i) {
+        path.push_back(i);
+    });
+
+    return get_kmer_coordinates(path, num_top_labels, presence_ratio);
+}
+
+std::vector<std::pair<std::string, std::vector<SmallVector<uint64_t>>>>
+AnnotatedDBG::get_kmer_coordinates(const std::vector<node_index> &path,
+                                   size_t num_top_labels,
+                                   double presence_ratio) const {
+    assert(presence_ratio >= 0.);
+    assert(presence_ratio <= 1.);
+    assert(check_compatibility());
+
+    if (!path.size())
+        return {};
+
+    std::vector<row_index> rows;
+    rows.reserve(path.size());
+
+    std::vector<size_t> ids;
+    ids.reserve(path.size());
+
+    for (node_index i : path) {
+        if (i > 0) {
+            ids.push_back(rows.size());
+            rows.push_back(graph_to_anno_index(i));
+        } else {
+            ids.push_back(-1);
+        }
+    }
+
+    uint64_t min_count = std::max(1.0, std::ceil(presence_ratio * path.size()));
+    if (rows.size() < min_count)
+        return {};
+
+    const auto *tuple_matrix = dynamic_cast<const MultiIntMatrix *>(&annotator_->get_matrix());
+    if (!tuple_matrix) {
+        logger->error("k-mer coordinates are not indexed in this annotator");
+        exit(1);
+    }
+
+    auto rows_tuples = tuple_matrix->get_row_tuples(rows);
+
+    VectorMap<size_t, size_t> code_to_count;
+    for (const auto &row_tuples : rows_tuples) {
+        for (const auto &[column, tuple] : row_tuples) {
+            code_to_count[column] += 1;
+        }
+    }
+
+    auto code_counts = code_to_count.values_container();
+    // sort by the number of matched k-mers
+    std::sort(code_counts.begin(), code_counts.end(),
+              [](const auto &x, const auto &y) {
+                  return x.second > y.second || (x.second == y.second && x.first < y.first);
+              });
+
+    // keep only the first |num_top_labels| top labels
+    if (code_counts.size() > num_top_labels)
+        code_counts.resize(num_top_labels);
+
+    // filter by the number of matched k-mers
+    code_counts.erase(
+        std::upper_bound(code_counts.begin(), code_counts.end(), min_count,
+                         [](uint64_t min, const auto &x) { return x.second < min; }),
+        code_counts.end()
+    );
+
+    code_to_count = VectorMap<size_t, size_t>(code_counts.begin(), code_counts.end());
+
+    std::vector<std::pair<std::string, std::vector<SmallVector<uint64_t>>>> result(code_to_count.size());
+
+    for (size_t j = 0; j < result.size(); ++j) {
+        result[j].first = annotator_->get_label_encoder().decode(code_counts[j].first);
+    }
+
+    for (size_t i : ids) {
+        for (size_t j = 0; j < result.size(); ++j) {
+            // append empty tuple
+            result[j].second.emplace_back();
+        }
+
+        // leave all tuples empty if the k-mer is missing
+        if (i == (size_t)-1)
+            continue;
+
+        // set the non-empty tuples
+        for (auto &[j, tuple] : rows_tuples[i]) {
+            auto it = code_to_count.find(j);
+            if (it != code_to_count.end())
+                result[it - code_to_count.begin()].second.back() = std::move(tuple);
+        }
+    }
+
+    return result;
+}
+
 std::vector<std::pair<Label, sdsl::bit_vector>>
 AnnotatedDBG::get_top_label_signatures(std::string_view sequence,
                                        size_t num_top_labels,
diff --git a/metagraph/src/graph/annotated_dbg.hpp b/metagraph/src/graph/annotated_dbg.hpp
index 2142f74c7b..1ffd58d4d8 100644
--- a/metagraph/src/graph/annotated_dbg.hpp
+++ b/metagraph/src/graph/annotated_dbg.hpp
@@ -9,6 +9,7 @@
 
 #include "representation/base/sequence_graph.hpp"
 #include "annotation/representation/base/annotation.hpp"
+#include "common/vector.hpp"
 
 
 namespace mtg {
@@ -116,6 +117,22 @@ class AnnotatedDBG : public AnnotatedSequenceGraph {
                    size_t min_count = 0,
                    bool with_kmer_counts = false) const;
 
+    std::vector<std::pair<Label, std::vector<size_t>>>
+    get_label_count_quantiles(std::string_view sequence,
+                              size_t num_top_labels,
+                              double presence_ratio,
+                              const std::vector<double> &count_quantiles) const;
+
+    std::vector<std::pair<Label, std::vector<SmallVector<uint64_t>>>>
+    get_kmer_coordinates(std::string_view sequence,
+                         size_t num_top_labels,
+                         double presence_ratio) const;
+
+    std::vector<std::pair<Label, std::vector<SmallVector<uint64_t>>>>
+    get_kmer_coordinates(const std::vector<node_index> &path,
+                         size_t num_top_labels,
+                         double presence_ratio) const;
+
     std::vector<std::pair<Label, sdsl::bit_vector>>
     get_top_label_signatures(std::string_view sequence,
                              size_t num_top_labels,
diff --git a/metagraph/src/graph/representation/canonical_dbg.cpp b/metagraph/src/graph/representation/canonical_dbg.cpp
index de33089f54..e5567bc107 100644
--- a/metagraph/src/graph/representation/canonical_dbg.cpp
+++ b/metagraph/src/graph/representation/canonical_dbg.cpp
@@ -15,18 +15,35 @@ CanonicalDBG::CanonicalDBG(std::shared_ptr<const DeBruijnGraph> graph, size_t ca
       : const_graph_ptr_(graph),
         offset_(graph_.max_index()),
         k_odd_(graph_.get_k() % 2),
+        has_sentinel_(false),
         alphabet_encoder_({ graph_.alphabet().size() }),
-        child_node_cache_(cache_size),
-        parent_node_cache_(cache_size),
-        is_palindrome_cache_(k_odd_ ? 0 : cache_size) {
-    if (graph->get_mode() != DeBruijnGraph::PRIMARY)
-        throw std::runtime_error("Only primary graphs can be wrapped in CanonicalDBG");
+        cache_size_(cache_size),
+        child_node_cache_(cache_size_),
+        parent_node_cache_(cache_size_),
+        is_palindrome_cache_(k_odd_ ? 0 : cache_size_) {
+    if (graph->get_mode() != DeBruijnGraph::PRIMARY) {
+        logger->error("Only primary graphs can be wrapped in CanonicalDBG");
+        exit(1);
+    }
 
     for (size_t i = 0; i < graph_.alphabet().size(); ++i) {
         alphabet_encoder_[graph_.alphabet()[i]] = i;
+        if (graph_.alphabet()[i] == boss::BOSS::kSentinel)
+            has_sentinel_ = true;
     }
 }
 
+CanonicalDBG::CanonicalDBG(const CanonicalDBG &canonical)
+      : const_graph_ptr_(canonical.const_graph_ptr_),
+        offset_(canonical.offset_),
+        k_odd_(canonical.k_odd_),
+        has_sentinel_(canonical.has_sentinel_),
+        alphabet_encoder_(canonical.alphabet_encoder_),
+        cache_size_(canonical.cache_size_),
+        child_node_cache_(cache_size_),
+        parent_node_cache_(cache_size_),
+        is_palindrome_cache_(k_odd_ ? 0 : cache_size_) {}
+
 CanonicalDBG::CanonicalDBG(std::shared_ptr<DeBruijnGraph> graph, size_t cache_size)
       : CanonicalDBG(std::dynamic_pointer_cast<const DeBruijnGraph>(graph), cache_size) {
     graph_ptr_ = graph;
@@ -48,13 +65,13 @@ uint64_t CanonicalDBG::num_nodes() const {
 void CanonicalDBG
 ::add_sequence(std::string_view sequence,
                const std::function<void(node_index)> &on_insertion) {
-    if (!graph_ptr_)
-        throw std::runtime_error("add_sequence only supported for non-const graphs.");
+    assert(graph_ptr_ && "add_sequence only supported for non-const graphs.");
 
     graph_ptr_->add_sequence(sequence, on_insertion);
     offset_ = graph_.max_index();
     child_node_cache_.Clear();
     parent_node_cache_.Clear();
+    is_palindrome_cache_.Clear();
 }
 
 
@@ -182,8 +199,7 @@ void CanonicalDBG::append_next_rc_nodes(node_index node,
 
     // for each n, check for nAGCCA. If found, define and store the index for
     // TGGCTrc(n) as index(nAGCCA) + offset_
-    const auto *dbg_succ = dynamic_cast<const DBGSuccinct*>(&graph_);
-    if (dbg_succ) {
+    if (const auto *dbg_succ = dynamic_cast<const DBGSuccinct*>(&graph_)) {
         const auto &boss = dbg_succ->get_boss();
         dbg_succ->call_nodes_with_suffix_matching_longest_prefix(
             std::string_view(&rev_seq[1], get_k() - 1),
@@ -195,21 +211,35 @@ void CanonicalDBG::append_next_rc_nodes(node_index node,
 
                 c = kmer::KmerExtractorBOSS::complement(c);
 
-                if (children[c] == npos)
+                if (children[c] == npos) {
                     children[c] = next + offset_;
+                    return;
+                }
+
+                if (k_odd_) {
+                    logger->error(
+                        "Primary graph contains both forward and reverse complement: {} {} -> {} {}\t{} {}",
+                        node, graph_.get_node_sequence(node),
+                        children[c], graph_.get_node_sequence(children[c]),
+                        next, graph_.get_node_sequence(next));
+                    exit(1);
+                }
+
+                is_palindrome_cache_.Put(next, true);
             },
             get_k() - 1
         );
 
     } else {
         for (size_t c = 0; c < alphabet.size(); ++c) {
-            if (children[c] != npos)
-                continue;
-
-            rev_seq[0] = complement(alphabet[c]);
-            node_index next = graph_.kmer_to_node(rev_seq);
-            if (next != npos)
-                children[c] = next + offset_;
+            // Do the checks by directly mapping the sequences of the desired k-mers.
+            // For non-DBGSuccinct graphs, this should be fast enough.
+            if (alphabet[c] != boss::BOSS::kSentinel && children[c] == npos) {
+                rev_seq[0] = complement(alphabet[c]);
+                node_index next = graph_.kmer_to_node(rev_seq);
+                if (next != npos)
+                    children[c] = next + offset_;
+            }
         }
     }
 }
@@ -229,22 +259,21 @@ ::call_outgoing_kmers(node_index node, const OutgoingEdgeCallback &callback) con
 
     const auto &alphabet = graph_.alphabet();
 
-    try {
-        auto children = child_node_cache_.Get(node);
+    if (auto fetch = child_node_cache_.TryGet(node)) {
         for (size_t c = 0; c < alphabet.size(); ++c) {
-            if (children[c] != npos)
-                callback(children[c], alphabet[c]);
+            if ((*fetch)[c] != npos)
+                callback((*fetch)[c], alphabet[c]);
         }
 
-    } catch (...) {
-        size_t max_num_edges_left = alphabet.size();
-        std::vector<node_index> children(max_num_edges_left);
+    } else {
+        std::vector<node_index> children(alphabet.size(), npos);
+        size_t max_num_edges_left = children.size() - has_sentinel_;
 
         graph_.call_outgoing_kmers(node, [&](node_index next, char c) {
-            if (c != boss::BOSS::kSentinel)
+            if (c != boss::BOSS::kSentinel) {
                 children[alphabet_encoder_[c]] = next;
-
-            --max_num_edges_left;
+                --max_num_edges_left;
+            }
         });
 
         if (max_num_edges_left)
@@ -301,21 +330,36 @@ void CanonicalDBG::append_prev_rc_nodes(node_index node,
 
                     c = kmer::KmerExtractorBOSS::complement(c);
 
-                    if (parents[c] == npos)
+                    if (parents[c] == npos) {
                         parents[c] = prev + offset_;
+                        return;
+                    }
+
+                    if (k_odd_) {
+                        logger->error(
+                            "Primary graph contains both forward and reverse complement: {} {} -> {} {}\t{} {}",
+                            node, graph_.get_node_sequence(node),
+                            parents[c], graph_.get_node_sequence(parents[c]),
+                            prev, graph_.get_node_sequence(prev)
+                        );
+                        exit(1);
+                    }
+
+                    is_palindrome_cache_.Put(prev, true);
                 }
             });
         }
 
     } else {
         for (size_t c = 0; c < alphabet.size(); ++c) {
-            if (parents[c] != npos)
-                continue;
-
-            rev_seq.back() = complement(alphabet[c]);
-            node_index prev = graph_.kmer_to_node(rev_seq);
-            if (prev != npos)
-                parents[c] = prev + offset_;
+            // Do the checks by directly mapping the sequences of the desired k-mers.
+            // For non-DBGSuccinct graphs, this should be fast enough.
+            if (alphabet[c] != boss::BOSS::kSentinel && parents[c] == npos) {
+                rev_seq.back() = complement(alphabet[c]);
+                node_index prev = graph_.kmer_to_node(rev_seq);
+                if (prev != npos)
+                    parents[c] = prev + offset_;
+            }
         }
     }
 }
@@ -335,22 +379,21 @@ ::call_incoming_kmers(node_index node, const IncomingEdgeCallback &callback) con
 
     const auto &alphabet = graph_.alphabet();
 
-    try {
-        auto parents = parent_node_cache_.Get(node);
+    if (auto fetch = parent_node_cache_.TryGet(node)) {
         for (size_t c = 0; c < alphabet.size(); ++c) {
-            if (parents[c] != npos)
-                callback(parents[c], alphabet[c]);
+            if ((*fetch)[c] != npos)
+                callback((*fetch)[c], alphabet[c]);
         }
 
-    } catch (...) {
-        size_t max_num_edges_left = alphabet.size();
-        std::vector<node_index> parents(max_num_edges_left);
+    } else {
+        std::vector<node_index> parents(alphabet.size(), npos);
+        size_t max_num_edges_left = parents.size() - has_sentinel_;
 
         graph_.call_incoming_kmers(node, [&](node_index prev, char c) {
-            if (c != boss::BOSS::kSentinel)
+            if (c != boss::BOSS::kSentinel) {
                 parents[alphabet_encoder_[c]] = prev;
-
-            --max_num_edges_left;
+                --max_num_edges_left;
+            }
         });
 
         if (max_num_edges_left)
@@ -499,10 +542,10 @@ DeBruijnGraph::node_index CanonicalDBG::reverse_complement(node_index node) cons
 
     }
 
-    try {
-        return is_palindrome_cache_.Get(node) ? node : node + offset_;
+    if (auto fetch = is_palindrome_cache_.TryGet(node)) {
+        return *fetch ? node : node + offset_;
 
-    } catch (...) {
+    } else {
         std::string seq = graph_.get_node_sequence(node);
         std::string rev_seq = seq;
         ::reverse_complement(rev_seq.begin(), rev_seq.end());
@@ -524,7 +567,7 @@ void CanonicalDBG::reverse_complement(std::string &seq,
 
     std::vector<node_index> rev_path(path.size());
     std::transform(path.begin(), path.end(), rev_path.rbegin(), [&](node_index i) {
-        return reverse_complement(i);
+        return i ? reverse_complement(i) : i;
     });
     std::swap(path, rev_path);
 }
diff --git a/metagraph/src/graph/representation/canonical_dbg.hpp b/metagraph/src/graph/representation/canonical_dbg.hpp
index bd576e4c51..10d1d769d4 100644
--- a/metagraph/src/graph/representation/canonical_dbg.hpp
+++ b/metagraph/src/graph/representation/canonical_dbg.hpp
@@ -46,6 +46,17 @@ class CanonicalDBG : public DeBruijnGraph {
      */
     CanonicalDBG(std::shared_ptr<DeBruijnGraph> graph, size_t cache_size = 100'000);
 
+    /**
+     * Copy constructor for CanonicalDBG. This creates a new wrapper with empty caches.
+     * @param canonical the graph to copy
+     */
+    CanonicalDBG(const CanonicalDBG &canonical);
+
+    // caches cannot be resized or moved, so disable these constructors
+    CanonicalDBG& operator=(const CanonicalDBG &canonical) = delete;
+    CanonicalDBG(CanonicalDBG&&) = delete;
+    CanonicalDBG& operator=(CanonicalDBG&&) = delete;
+
     virtual ~CanonicalDBG() {}
 
     virtual void add_sequence(std::string_view sequence,
@@ -144,11 +155,14 @@ class CanonicalDBG : public DeBruijnGraph {
     const DeBruijnGraph &graph_ = *const_graph_ptr_;
     size_t offset_;
     bool k_odd_;
+    bool has_sentinel_;
 
     std::shared_ptr<DeBruijnGraph> graph_ptr_;
 
     std::array<size_t, 256> alphabet_encoder_;
 
+    size_t cache_size_;
+
     // cache the results of call_outgoing_kmers
     mutable caches::fixed_sized_cache<node_index, std::vector<node_index>,
                                       caches::LRUCachePolicy<node_index>> child_node_cache_;
diff --git a/metagraph/tests/graph/all/test_dbg_helpers.hpp b/metagraph/tests/graph/all/test_dbg_helpers.hpp
index 7a008935ea..c43234f986 100644
--- a/metagraph/tests/graph/all/test_dbg_helpers.hpp
+++ b/metagraph/tests/graph/all/test_dbg_helpers.hpp
@@ -111,14 +111,7 @@ typedef ::testing::Types<DBGBitmap,
                          DBGSuccinctBloom<4, 1>,
                          DBGSuccinctBloom<4, 50>> StableGraphTypes;
 
-typedef ::testing::Types<DBGHashFast,
-                         DBGSuccinct,
-#if _PROTEIN_GRAPH
-                         DBGSuccinctIndexed<2>,
-#else
-                         DBGSuccinctIndexed<10>,
-#endif
-                         DBGSuccinctBloomFPR<1, 10>> FewGraphTypes;
+typedef ::testing::Types<DBGHashFast, DBGSuccinct> FewGraphTypes;
 
 } // namespace test
 } // namespace mtg
diff --git a/metagraph/tests/graph/test_aligner.cpp b/metagraph/tests/graph/test_aligner.cpp
index b147452a68..4e02e48f7d 100644
--- a/metagraph/tests/graph/test_aligner.cpp
+++ b/metagraph/tests/graph/test_aligner.cpp
@@ -277,8 +277,7 @@ TYPED_TEST(DBGAlignerTest, align_straight_forward_and_reverse_complement) {
     check_json_dump_load(*graph, path, paths.get_query(), paths.get_query(PICK_REV_COMP));
 
     check_extend(graph, aligner.get_config(), paths, query);
-
-    auto ext_paths = DBGAligner<UniMEMSeeder<>>(*graph, config_fwd_and_rev).align(query);
+    auto ext_paths = get_extend(graph, config_fwd_and_rev, paths, query);
 
     EXPECT_TRUE(std::equal(paths.begin(), paths.end(),
                            ext_paths.begin(), ext_paths.end()));
diff --git a/metagraph/tests/graph/test_aligner_helpers.hpp b/metagraph/tests/graph/test_aligner_helpers.hpp
index ff1c9d17b2..f1be4bb099 100644
--- a/metagraph/tests/graph/test_aligner_helpers.hpp
+++ b/metagraph/tests/graph/test_aligner_helpers.hpp
@@ -86,7 +86,18 @@ inline void check_extend(std::shared_ptr<const DeBruijnGraph> graph,
     ASSERT_EQ(paths.size(), unimem_paths.size());
 
     for (size_t i = 0; i < paths.size(); ++i) {
-        EXPECT_EQ(paths[i], unimem_paths[i]) << paths[i] << "\n" << unimem_paths[i];
+        if (graph->get_mode() == DeBruijnGraph::CANONICAL) {
+            auto alt = unimem_paths[i];
+            alt.reverse_complement(*graph, paths.get_query(!alt.get_orientation()));
+            if (alt.size()) {
+                EXPECT_TRUE(paths[i] == unimem_paths[i] || paths[i] == alt)
+                     << i << "\n" << paths[i] << "\n" << unimem_paths[i];
+            } else {
+                EXPECT_EQ(paths[i], unimem_paths[i]) << i << "\n" << paths[i] << "\n" << unimem_paths[i];
+            }
+        } else {
+            EXPECT_EQ(paths[i], unimem_paths[i]) << i << "\n" << paths[i] << "\n" << unimem_paths[i];
+        }
     }
 }
 
diff --git a/metagraph/tests/graph/test_canonical_dbg.cpp b/metagraph/tests/graph/test_canonical_dbg.cpp
index b2b4b27787..ebfa255c8e 100644
--- a/metagraph/tests/graph/test_canonical_dbg.cpp
+++ b/metagraph/tests/graph/test_canonical_dbg.cpp
@@ -558,8 +558,7 @@ TYPED_TEST(CanonicalDBGTest, CallUnitigs) {
 }
 
 TYPED_TEST(CanonicalDBGTest, WrapCanonicalGraphFail) {
-    EXPECT_THROW(CanonicalDBG(build_graph_batch<TypeParam>(3, {}, DeBruijnGraph::CANONICAL)),
-                 std::runtime_error);
+    ASSERT_DEATH(CanonicalDBG(build_graph_batch<TypeParam>(3, {}, DeBruijnGraph::CANONICAL)), "");
 }
 
 // TODO: A different combination of forward and reverse complement k-mers may be