diff --git a/metagraph/src/annotation/binary_matrix/base/binary_matrix.cpp b/metagraph/src/annotation/binary_matrix/base/binary_matrix.cpp
index 8b7ec4ea06..fb32ce927b 100644
--- a/metagraph/src/annotation/binary_matrix/base/binary_matrix.cpp
+++ b/metagraph/src/annotation/binary_matrix/base/binary_matrix.cpp
@@ -1,6 +1,7 @@
 #include "binary_matrix.hpp"
 
 #include <ips4o.hpp>
+#include <progress_bar.hpp>
 
 #include "common/vectors/bitmap.hpp"
 #include "common/serialization.hpp"
@@ -38,11 +39,14 @@ BinaryMatrix::slice_rows(const std::vector<Row> &row_ids) const {
 }
 
 void BinaryMatrix::call_columns(const std::vector<Column> &column_ids,
-                                const std::function<void(size_t, const bitmap&)> &callback,
+                                const std::function<void(size_t, const bitmap_generator&)> &callback,
                                 size_t num_threads) const {
+    ProgressBar progress_bar(column_ids.size(), "Parsing columns",
+                             std::cerr, !common::get_verbose());
     #pragma omp parallel for num_threads(num_threads) schedule(dynamic)
     for (size_t k = 0; k < column_ids.size(); ++k) {
         callback(k, bitmap_generator(get_column(column_ids[k]), num_rows()));
+        ++progress_bar;
     }
 }
 
diff --git a/metagraph/src/annotation/binary_matrix/base/binary_matrix.hpp b/metagraph/src/annotation/binary_matrix/base/binary_matrix.hpp
index b9d4815bc2..f5f8a166ce 100644
--- a/metagraph/src/annotation/binary_matrix/base/binary_matrix.hpp
+++ b/metagraph/src/annotation/binary_matrix/base/binary_matrix.hpp
@@ -9,7 +9,7 @@
 #include "common/vector.hpp"
 
 
-class bitmap;
+class bitmap_generator;
 
 namespace mtg {
 namespace annot {
@@ -41,7 +41,7 @@ class BinaryMatrix {
     // For each column id in columns, run callback on its respective index in columns
     // and a bitmap represnting the column
     virtual void call_columns(const std::vector<Column> &columns,
-                              const std::function<void(size_t, const bitmap&)> &callback,
+                              const std::function<void(size_t, const bitmap_generator&)> &callback,
                               size_t num_threads = 1) const;
 
     virtual bool load(std::istream &in) = 0;
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 acf930a484..30f4bf4a13 100644
--- a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp
+++ b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.cpp
@@ -104,7 +104,7 @@ ColumnMajor::slice_rows(const std::vector<Row> &row_ids) const {
 }
 
 void ColumnMajor::call_columns(const std::vector<Column> &columns,
-                               const std::function<void(size_t, const bitmap&)> &callback,
+                               const std::function<void(size_t, const bitmap_generator&)> &callback,
                                size_t num_threads) const {
     #pragma omp parallel for num_threads(num_threads) schedule(dynamic)
     for (size_t j = 0; j < columns.size(); ++j) {
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 a473e6f1d0..4d29932c2b 100644
--- a/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp
+++ b/metagraph/src/annotation/binary_matrix/column_sparse/column_major.hpp
@@ -32,7 +32,7 @@ class ColumnMajor : public BinaryMatrix {
     std::vector<Column> slice_rows(const std::vector<Row> &rows) const override;
 
     void call_columns(const std::vector<Column> &columns,
-                      const std::function<void(size_t, const bitmap&)> &callback,
+                      const std::function<void(size_t, const bitmap_generator&)> &callback,
                       size_t num_threads = 1) const override;
 
     bool load(std::istream &in) override;
diff --git a/metagraph/src/cli/assemble.cpp b/metagraph/src/cli/assemble.cpp
index d2e60bfd8d..6cc7deefc5 100644
--- a/metagraph/src/cli/assemble.cpp
+++ b/metagraph/src/cli/assemble.cpp
@@ -11,6 +11,7 @@
 #include "seq_io/sequence_io.hpp"
 #include "config/config.hpp"
 #include "load/load_graph.hpp"
+#include "load/load_annotation.hpp"
 #include "load/load_annotated_graph.hpp"
 #include "graph/annotated_graph_algorithm.hpp"
 #include "graph/representation/masked_graph.hpp"
@@ -41,10 +42,8 @@ void check_labels(const tsl::hopscotch_set<std::string> &label_set,
         exit(1);
 }
 
-DifferentialAssemblyConfig diff_assembly_config(const Json::Value &experiment,
-                                                const DeBruijnGraph &graph) {
+DifferentialAssemblyConfig diff_assembly_config(const Json::Value &experiment) {
     DifferentialAssemblyConfig diff_config;
-    diff_config.add_complement = graph.get_mode() == DeBruijnGraph::CANONICAL;
     diff_config.label_mask_in_kmer_fraction = experiment.get("in_min_fraction", 1.0).asDouble();
     diff_config.label_mask_in_unitig_fraction = experiment.get("unitig_in_min_fraction", 0.0).asDouble();
     diff_config.label_mask_out_kmer_fraction = experiment.get("out_max_fraction", 0.0).asDouble();
@@ -81,16 +80,11 @@ void call_masked_graphs(const AnnotatedDBG &anno_graph,
     Json::Value diff_json;
     fin >> diff_json;
 
-    tsl::hopscotch_set<std::string> foreground_labels;
-    tsl::hopscotch_set<std::string> background_labels;
-    tsl::hopscotch_set<std::string> shared_foreground_labels;
-    tsl::hopscotch_set<std::string> shared_background_labels;
-
     for (const Json::Value &group : diff_json["groups"]) {
-        if (group["shared_labels"]) {
-            shared_foreground_labels.clear();
-            shared_background_labels.clear();
+        tsl::hopscotch_set<std::string> shared_foreground_labels;
+        tsl::hopscotch_set<std::string> shared_background_labels;
 
+        if (group["shared_labels"]) {
             for (const Json::Value &in_label : group["shared_labels"]["in"]) {
                 shared_foreground_labels.emplace(in_label.asString());
             }
@@ -107,12 +101,13 @@ void call_masked_graphs(const AnnotatedDBG &anno_graph,
             throw std::runtime_error("Missing experiments in group");
 
         for (const Json::Value &experiment : group["experiments"]) {
-            DifferentialAssemblyConfig diff_config = diff_assembly_config(
-                experiment, anno_graph.get_graph()
-            );
+            tsl::hopscotch_set<std::string> foreground_labels;
+            tsl::hopscotch_set<std::string> background_labels;
 
-            foreground_labels.clear();
-            background_labels.clear();
+            DifferentialAssemblyConfig diff_config = diff_assembly_config(experiment);
+
+            std::string exp_name = experiment["name"].asString()
+                                    + (config->enumerate_out_sequences ? "." : "");
 
             for (const Json::Value &in_label : experiment["in"]) {
                 foreground_labels.emplace(in_label.asString());
@@ -125,9 +120,6 @@ void call_masked_graphs(const AnnotatedDBG &anno_graph,
             check_labels(foreground_labels, anno_graph);
             check_labels(background_labels, anno_graph);
 
-            std::string exp_name = experiment["name"].asString()
-                                    + (config->enumerate_out_sequences ? "." : "");
-
             logger->trace("Running assembly: {}", exp_name);
 
             callback(*mask_nodes_by_label(anno_graph,
@@ -135,7 +127,65 @@ void call_masked_graphs(const AnnotatedDBG &anno_graph,
                                           background_labels,
                                           shared_foreground_labels,
                                           shared_background_labels,
-                                          diff_config, num_threads), exp_name);
+                                          diff_config, num_threads,
+                                          config->parallel_nodes), exp_name);
+        }
+    }
+}
+
+void call_masked_graphs(std::shared_ptr<const DeBruijnGraph> graph_ptr,
+                        Config *config,
+                        const CallMaskedGraphHeader &callback) {
+    assert(!config->assembly_config_file.empty());
+
+    std::ifstream fin(config->assembly_config_file);
+    if (!fin.good())
+        throw std::iostream::failure("Failed to read assembly config JSON from " + config->assembly_config_file);
+
+    for (const auto &name : config->fnames) {
+        if (parse_annotation_type(name) != Config::ColumnCompressed) {
+            throw std::runtime_error("Multiple annotation files must be ColumnCompressed");
+        }
+    }
+
+    size_t num_threads = std::max(1u, get_num_threads());
+
+    Json::Value diff_json;
+    fin >> diff_json;
+
+    for (const Json::Value &group : diff_json["groups"]) {
+        if (group["shared_labels"]) {
+            throw std::runtime_error("Shared labels not supported with multiple annotation files");
+        }
+
+        if (!group["experiments"])
+            throw std::runtime_error("Missing experiments in group");
+
+        for (const Json::Value &experiment : group["experiments"]) {
+            tsl::hopscotch_set<std::string> foreground_labels;
+            tsl::hopscotch_set<std::string> background_labels;
+
+            DifferentialAssemblyConfig diff_config = diff_assembly_config(experiment);
+
+            std::string exp_name = experiment["name"].asString()
+                                    + (config->enumerate_out_sequences ? "." : "");
+
+            for (const Json::Value &in_label : experiment["in"]) {
+                foreground_labels.emplace(in_label.asString());
+            }
+
+            for (const Json::Value &out_label : experiment["out"]) {
+                background_labels.emplace(out_label.asString());
+            }
+
+            logger->trace("Running assembly: {}", exp_name);
+
+            callback(*mask_nodes_by_label(graph_ptr,
+                                          config->fnames,
+                                          foreground_labels,
+                                          background_labels,
+                                          diff_config, num_threads,
+                                          config->parallel_nodes), exp_name);
         }
     }
 }
@@ -143,9 +193,17 @@ void call_masked_graphs(const AnnotatedDBG &anno_graph,
 int assemble(Config *config) {
     assert(config);
 
+    if (config->infbase_annotators.size()) {
+        assert(config->assembly_config_file.size());
+        if (!std::filesystem::exists(config->assembly_config_file)) {
+            logger->error("Differential assembly config does not exist\n{}",
+                          config->assembly_config_file);
+            exit(1);
+        }
+    }
+
     const auto &files = config->fnames;
 
-    assert(files.size() == 1);
     assert(config->outfbase.size());
 
     Timer timer;
@@ -156,48 +214,49 @@ int assemble(Config *config) {
     logger->trace("Graph loaded in {} sec", timer.elapsed());
 
     if (config->infbase_annotators.size()) {
-        config->infbase = files.at(0);
-
-        assert(config->assembly_config_file.size());
-        auto anno_graph = initialize_annotated_dbg(graph, *config);
-
         logger->trace("Generating masked graphs...");
+        std::mutex write_mutex;
+        size_t num_threads = std::max(1u, get_num_threads());
 
         std::filesystem::remove(
             utils::remove_suffix(config->outfbase, ".gz", ".fasta") + ".fasta.gz"
         );
 
-        std::mutex write_mutex;
-
-        size_t num_threads = std::max(1u, get_num_threads());
-
-        call_masked_graphs(*anno_graph, config,
-            [&](const graph::MaskedDeBruijnGraph &graph, const std::string &header) {
-                seq_io::FastaWriter writer(config->outfbase, header,
-                                           config->enumerate_out_sequences,
-                                           num_threads > 1, /* async write */
-                                           "a" /* append mode */);
-
-                if (config->unitigs || config->min_tip_size > 1) {
-                    graph.call_unitigs([&](const std::string &unitig, auto&&) {
-                                           std::lock_guard<std::mutex> lock(write_mutex);
-                                           writer.write(unitig);
-                                       },
-                                       num_threads, config->min_tip_size,
-                                       config->kmers_in_single_form);
-                } else {
-                    graph.call_sequences([&](const std::string &seq, auto&&) {
-                                             std::lock_guard<std::mutex> lock(write_mutex);
-                                             writer.write(seq);
-                                         },
-                                         num_threads, config->kmers_in_single_form);
-                }
+        auto graph_callback = [&](const graph::MaskedDeBruijnGraph &graph, const std::string &header) {
+            seq_io::FastaWriter writer(config->outfbase, header,
+                                       config->enumerate_out_sequences,
+                                       num_threads > 1, /* async write */
+                                       "a" /* append mode */);
+
+            if (config->unitigs || config->min_tip_size > 1) {
+                graph.call_unitigs([&](const std::string &unitig, auto&&) {
+                                       std::lock_guard<std::mutex> lock(write_mutex);
+                                       writer.write(unitig);
+                                   },
+                                   num_threads, config->min_tip_size,
+                                   config->kmers_in_single_form);
+            } else {
+                graph.call_sequences([&](const std::string &seq, auto&&) {
+                                         std::lock_guard<std::mutex> lock(write_mutex);
+                                         writer.write(seq);
+                                     },
+                                     num_threads, config->kmers_in_single_form);
             }
-        );
-
+        };
+
+        if (files.size() > 1) {
+            config->fnames.erase(config->fnames.begin(), config->fnames.begin() + 1);
+            call_masked_graphs(graph, config, graph_callback);
+        } else {
+            config->infbase = files.at(0);
+            auto anno_graph = initialize_annotated_dbg(graph, *config);
+            call_masked_graphs(*anno_graph, config, graph_callback);
+        }
         return 0;
     }
 
+    assert(files.size() == 1);
+
     logger->trace("Extracting sequences from graph...");
 
     timer.reset();
diff --git a/metagraph/src/common/vectors/bitmap.cpp b/metagraph/src/common/vectors/bitmap.cpp
index f4888735e4..6d31c5f7c8 100644
--- a/metagraph/src/common/vectors/bitmap.cpp
+++ b/metagraph/src/common/vectors/bitmap.cpp
@@ -356,6 +356,13 @@ uint64_t bitmap_lazy::num_set_bits() const {
     return num_set_bits_;
 }
 
+std::optional<uint64_t> bitmap_lazy::try_get_num_set_bits() const {
+    if (num_set_bits_ == static_cast<size_t>(-1))
+        return std::nullopt;
+
+    return num_set_bits_;
+}
+
 void bitmap_lazy::add_to(sdsl::bit_vector *other) const {
     assert(other);
     assert(other->size() == size());
diff --git a/metagraph/src/common/vectors/bitmap.hpp b/metagraph/src/common/vectors/bitmap.hpp
index b5b855a5eb..9f8d68438a 100644
--- a/metagraph/src/common/vectors/bitmap.hpp
+++ b/metagraph/src/common/vectors/bitmap.hpp
@@ -176,6 +176,7 @@ class bitmap_lazy : public bitmap {
 
     uint64_t size() const { return size_; }
     uint64_t num_set_bits() const;
+    std::optional<uint64_t> try_get_num_set_bits() const;
 
     void add_to(sdsl::bit_vector *other) const;
     void call_ones_in_range(uint64_t begin, uint64_t end,
@@ -197,12 +198,18 @@ class bitmap_generator : public bitmap {
                            std::for_each(s.begin(), s.end(), callback);
                        }) {}
 
+    bitmap_generator(const bitmap &bits) noexcept
+          : size_(bits.size()), num_set_bits_(bits.num_set_bits()),
+            generator_([&bits](const VoidCall<uint64_t> &callback) {
+                           bits.call_ones(callback);
+                       }) {}
+
     bitmap_generator(std::function<void(const VoidCall<uint64_t>&)>&& generator,
                      size_t size,
                      size_t num_set_bits = -1) noexcept;
 
-    bool operator[](uint64_t) const { throw std::runtime_error("Not implemented"); }
-    uint64_t get_int(uint64_t, uint32_t) const { throw std::runtime_error("Not implemented"); }
+    bool operator[](uint64_t) const { throw std::runtime_error("operator[] not implemented for bitmap_generator"); }
+    uint64_t get_int(uint64_t, uint32_t) const { throw std::runtime_error("get_int not implemented for bitmap_generator"); }
 
     uint64_t size() const { return size_; }
     uint64_t num_set_bits() const { return num_set_bits_; }
diff --git a/metagraph/src/graph/annotated_graph_algorithm.cpp b/metagraph/src/graph/annotated_graph_algorithm.cpp
index 60dc814919..0a5ee8d1f7 100644
--- a/metagraph/src/graph/annotated_graph_algorithm.cpp
+++ b/metagraph/src/graph/annotated_graph_algorithm.cpp
@@ -4,6 +4,7 @@
 #include "common/vectors/vector_algorithm.hpp"
 #include "common/vectors/bitmap.hpp"
 #include "graph/representation/masked_graph.hpp"
+#include "annotation/representation/column_compressed/annotate_column_compressed.hpp"
 
 
 namespace mtg {
@@ -12,7 +13,9 @@ namespace graph {
 using mtg::common::logger;
 
 typedef AnnotatedDBG::node_index node_index;
+typedef AnnotatedDBG::Annotator Annotator;
 typedef AnnotatedDBG::Annotator::Label Label;
+using Column = annot::binmat::BinaryMatrix::Column;
 
 typedef std::function<size_t()> LabelCountCallback;
 
@@ -20,21 +23,36 @@ constexpr std::memory_order MO_RELAXED = std::memory_order_relaxed;
 
 
 /**
- * Return an int_vector<>, bit_vector pair, of lengths anno_graph.get_graph().max_index() * 2
- * and anno_graph.get_graph().max_index(), respectively.
+ * Return an int_vector<>, bit_vector of lengths anno_graph.get_graph().max_index() * 2
+ * and anno_graph.get_graph().max_index(), respectively, and a bit_vector of length
+ * equal to the total number of labels.
  * For an index i, the int_vector at indices 2*i and 2*i + 1 represent the
  * number of labels in labels_in and labels_out which the k-mer of index i is
  * annotated with, respectively. The width of the int_vector<> is computed to be
  * wide enough to contain counts up to num_labels.
  * The returned bit_vector is a k-mer mask indicating those k-mers annotated
- * with at least one in-label or out-label.
+ * with at least one in-label. If add_out_labels_to_mask is true, then it also indicates
+ * which those k-mers with at least one out-label.
+ * The second bit vector indicates which of the labels correspond to an "other"
+ * label (i.e., not in- or out-).
  */
-std::pair<sdsl::int_vector<>, sdsl::bit_vector>
-construct_diff_label_count_vector(const AnnotatedDBG &anno_graph,
+std::tuple<sdsl::int_vector<>, sdsl::bit_vector, sdsl::bit_vector>
+construct_diff_label_count_vector(const Annotator &annotator,
                                   const tsl::hopscotch_set<Label> &labels_in,
                                   const tsl::hopscotch_set<Label> &labels_out,
+                                  size_t max_index,
                                   size_t num_labels,
-                                  size_t num_threads);
+                                  size_t num_threads,
+                                  bool add_out_labels_to_mask);
+
+std::tuple<sdsl::int_vector<>, sdsl::bit_vector, sdsl::bit_vector, std::vector<size_t>>
+construct_diff_label_count_vector(const std::vector<std::string> &files,
+                                  const tsl::hopscotch_set<Label> &labels_in,
+                                  const tsl::hopscotch_set<Label> &labels_out,
+                                  size_t max_index,
+                                  size_t num_labels,
+                                  size_t num_threads,
+                                  bool add_out_labels_to_mask);
 
 // Regions of a graph mask which should be kept (i.e., masked in)
 typedef std::vector<std::pair<size_t, size_t>> Intervals;
@@ -57,6 +75,20 @@ make_initial_masked_graph(std::shared_ptr<const DeBruijnGraph> graph_ptr,
                           bool add_complement,
                           size_t num_threads);
 
+// Helper function for calls in interface
+std::shared_ptr<MaskedDeBruijnGraph>
+mask_nodes_by_label(std::shared_ptr<const DeBruijnGraph> graph_ptr,
+                    const AnnotatedDBG *anno_graph,
+                    sdsl::int_vector<>&& counts,
+                    sdsl::bit_vector&& mask,
+                    sdsl::bit_vector&& other_mask,
+                    size_t num_labels,
+                    const tsl::hopscotch_set<Label> &labels_in,
+                    const tsl::hopscotch_set<Label> &labels_out,
+                    const tsl::hopscotch_set<Label> &labels_in_round2,
+                    const tsl::hopscotch_set<Label> &labels_out_round2,
+                    const DifferentialAssemblyConfig &config,
+                    size_t num_threads);
 
 std::shared_ptr<MaskedDeBruijnGraph>
 mask_nodes_by_label(const AnnotatedDBG &anno_graph,
@@ -65,52 +97,152 @@ mask_nodes_by_label(const AnnotatedDBG &anno_graph,
                     const tsl::hopscotch_set<Label> &labels_in_round2,
                     const tsl::hopscotch_set<Label> &labels_out_round2,
                     const DifferentialAssemblyConfig &config,
-                    size_t num_threads) {
-    bool parallel = num_threads > 1;
-    size_t num_labels = anno_graph.get_annotator().num_labels();
+                    size_t num_threads,
+                    size_t num_parallel_files) {
+    num_parallel_files = std::min(num_threads, num_parallel_files);
     size_t num_in_labels = labels_in.size() + labels_in_round2.size();
     size_t num_out_labels = labels_out.size() + labels_out_round2.size();
-    auto graph_ptr = std::dynamic_pointer_cast<const DeBruijnGraph>(
+    auto graph_ptr = std::static_pointer_cast<const DeBruijnGraph>(
         anno_graph.get_graph_ptr()
     );
 
+    bool check_other = config.label_mask_other_unitig_fraction != 1.0;
+    bool unitig_mode = check_other || labels_in_round2.size() || labels_out_round2.size()
+            || config.label_mask_in_unitig_fraction != 0.0
+            || config.label_mask_out_unitig_fraction != 1.0
+            || config.label_mask_other_unitig_fraction != 1.0;
+
+    bool add_complement = graph_ptr->get_mode() == DeBruijnGraph::CANONICAL
+        && (config.add_complement || unitig_mode);
+
     logger->trace("Generating initial mask");
 
     // Construct initial masked graph from union of labels in labels_in
     auto count_vector = construct_diff_label_count_vector(
-        anno_graph, labels_in, labels_out,
-        std::max(num_in_labels, num_out_labels),
-        num_threads
+        anno_graph.get_annotator(), labels_in, labels_out, graph_ptr->max_index(),
+        std::max(num_in_labels, num_out_labels), num_parallel_files,
+        add_complement
     );
-    auto &[counts, init_mask] = count_vector;
+    auto &[counts, init_mask, other_labels] = count_vector;
+    sdsl::bit_vector other_mask(init_mask.size() * check_other, false);
+    return mask_nodes_by_label(graph_ptr, &anno_graph,
+                               std::move(counts), std::move(init_mask),
+                               std::move(other_mask),
+                               anno_graph.get_annotator().num_labels(),
+                               labels_in, labels_out,
+                               labels_in_round2, labels_out_round2,
+                               config, num_threads);
+}
+
+std::shared_ptr<MaskedDeBruijnGraph>
+mask_nodes_by_label(std::shared_ptr<const DeBruijnGraph> graph_ptr,
+                    const std::vector<std::string> &files,
+                    const tsl::hopscotch_set<Label> &labels_in,
+                    const tsl::hopscotch_set<Label> &labels_out,
+                    const DifferentialAssemblyConfig &config,
+                    size_t num_threads,
+                    size_t num_parallel_files) {
+    bool check_other = config.label_mask_other_unitig_fraction != 1.0;
+    bool unitig_mode = check_other
+            || config.label_mask_in_unitig_fraction != 0.0
+            || config.label_mask_out_unitig_fraction != 1.0
+            || config.label_mask_other_unitig_fraction != 1.0;
+
+    bool add_complement = graph_ptr->get_mode() == DeBruijnGraph::CANONICAL
+        && (config.add_complement || unitig_mode);
 
+    logger->trace("Generating initial mask");
+
+    // Construct initial masked graph from union of labels in labels_in
+    auto count_vector = construct_diff_label_count_vector(
+        files, labels_in, labels_out, graph_ptr->max_index(),
+        std::max(labels_in.size(), labels_out.size()), num_parallel_files,
+        add_complement
+    );
+    auto &[counts, init_mask, other_labels, num_labels_per_file] = count_vector;
+    sdsl::bit_vector other_mask(init_mask.size() * check_other, false);
+    if (check_other && sdsl::util::cnt_one_bits(other_labels)) {
+        bool parallel = num_parallel_files > 1;
+        size_t j = 0;
+        std::atomic_thread_fence(std::memory_order_release);
+        for (size_t i = 0; i < files.size(); ++i) {
+            if (count_ones(other_labels, j, j + num_labels_per_file[i])) {
+                annot::ColumnCompressed<>::merge_load({ files[i] },
+                    [&](size_t offset, const auto &, auto&& column) {
+                        auto &[counts, init_mask, other_labels, num_labels_per_file]
+                            = count_vector;
+                        if (other_labels[j + offset]) {
+                            call_ones(init_mask, [&](size_t i) {
+                                if ((*column)[AnnotatedDBG::graph_to_anno_index(i)])
+                                    set_bit(other_mask.data(), i, parallel, MO_RELAXED);
+                            });
+                        }
+                    },
+                    num_parallel_files
+                );
+            }
+
+            j += num_labels_per_file[i];
+        }
+        std::atomic_thread_fence(std::memory_order_acquire);
+    }
+
+    return mask_nodes_by_label(graph_ptr, nullptr,
+                               std::move(counts), std::move(init_mask),
+                               std::move(other_mask),
+                               files.size(),
+                               labels_in, labels_out, {}, {},
+                               config, num_threads);
+}
+
+std::shared_ptr<MaskedDeBruijnGraph>
+mask_nodes_by_label(std::shared_ptr<const DeBruijnGraph> graph_ptr,
+                    const AnnotatedDBG *anno_graph,
+                    sdsl::int_vector<>&& counts,
+                    sdsl::bit_vector&& init_mask,
+                    sdsl::bit_vector&& other_mask,
+                    size_t num_labels,
+                    const tsl::hopscotch_set<Label> &labels_in,
+                    const tsl::hopscotch_set<Label> &labels_out,
+                    const tsl::hopscotch_set<Label> &labels_in_round2,
+                    const tsl::hopscotch_set<Label> &labels_out_round2,
+                    const DifferentialAssemblyConfig &config,
+                    size_t num_threads) {
     // in and out counts are stored interleaved in the counts vector
     assert(counts.size() == init_mask.size() * 2);
 
+    bool parallel = num_threads > 1;
     bool check_other = config.label_mask_other_unitig_fraction != 1.0;
-    sdsl::bit_vector union_mask;
+    bool unitig_mode = check_other || labels_in_round2.size() || labels_out_round2.size()
+            || config.label_mask_in_unitig_fraction != 0.0
+            || config.label_mask_out_unitig_fraction != 1.0
+            || config.label_mask_other_unitig_fraction != 1.0;
 
-    sdsl::bit_vector other_mask(init_mask.size() * check_other, false);
-    auto masked_graph = make_initial_masked_graph(graph_ptr, counts, std::move(init_mask),
-                                                  config.add_complement, num_threads);
+    size_t num_in_labels = labels_in.size() + labels_in_round2.size();
+    size_t num_out_labels = labels_out.size() + labels_out_round2.size();
 
-    if (check_other || labels_in_round2.size() || labels_out_round2.size())
-        union_mask = dynamic_cast<const bitmap_vector&>(masked_graph->get_mask()).data();
+    bool add_complement = graph_ptr->get_mode() == DeBruijnGraph::CANONICAL
+        && (config.add_complement || unitig_mode);
 
-    // check all other labels and post labels
-    if (check_other || labels_in_round2.size() || labels_out_round2.size()) {
+    auto masked_graph = make_initial_masked_graph(graph_ptr, counts, std::move(init_mask),
+                                                  add_complement, num_threads);
+
+    auto mask_or = [&](sdsl::bit_vector &a,
+                       const sdsl::bit_vector &b,
+                       const std::vector<node_index> &id_map) {
+        call_ones(b, [&](size_t i) {
+            if (id_map[i])
+                set_bit(a.data(), id_map[i], parallel, MO_RELAXED);
+        });
+    };
+
+    // check all other labels and round 2 labels
+    if (anno_graph && (check_other || labels_in_round2.size() || labels_out_round2.size())) {
+        sdsl::bit_vector union_mask
+            = static_cast<const bitmap_vector&>(masked_graph->get_mask()).data();
         std::mutex vector_backup_mutex;
         std::atomic_thread_fence(std::memory_order_release);
 
-        auto mask_or = [&](sdsl::bit_vector &a,
-                           const sdsl::bit_vector &b,
-                           const std::vector<node_index> &id_map) {
-            call_ones(b, [&](size_t i) {
-                if (id_map[i])
-                    set_bit(a.data(), id_map[i], parallel, MO_RELAXED);
-            });
-        };
-
         auto count_merge = [&](sdsl::bit_vector &a,
                                const sdsl::bit_vector &b,
                                const std::vector<node_index> &id_map,
@@ -118,15 +250,16 @@ mask_nodes_by_label(const AnnotatedDBG &anno_graph,
             call_ones(b, [&](size_t i) {
                 if (id_map[i]) {
                     set_bit(a.data(), id_map[i], parallel, MO_RELAXED);
-                    atomic_fetch_and_add(count_vector.first, id_map[i] * 2 + offset, 1,
+                    atomic_fetch_and_add(counts, id_map[i] * 2 + offset, 1,
                                          vector_backup_mutex, MO_RELAXED);
                 }
             });
         };
 
         logger->trace("Checking shared and other labels");
-        masked_graph->call_sequences([&](const std::string &contig, const std::vector<node_index> &path) {
-            for (const auto &[label, sig] : anno_graph.get_top_label_signatures(contig, num_labels)) {
+        masked_graph->call_sequences([&](const std::string &contig,
+                                         const std::vector<node_index> &path) {
+            for (const auto &[label, sig] : anno_graph->get_top_label_signatures(contig, num_labels)) {
                 bool found_in = labels_in.count(label);
                 bool found_out = labels_out.count(label);
                 bool found_in_round2 = labels_in_round2.count(label);
@@ -152,22 +285,25 @@ mask_nodes_by_label(const AnnotatedDBG &anno_graph,
     // Filter unitigs from masked graph based on filtration criteria
     logger->trace("Filtering out background");
 
-    size_t min_label_in_count = std::ceil(config.label_mask_in_kmer_fraction * num_in_labels);
-    size_t max_label_out_count = std::floor(config.label_mask_out_kmer_fraction * num_out_labels);
+    size_t min_label_in_count = std::ceil(config.label_mask_in_kmer_fraction
+                                    * num_in_labels);
+    size_t max_label_out_count = std::floor(config.label_mask_out_kmer_fraction
+                                    * num_out_labels);
 
     if (config.label_mask_in_unitig_fraction == 0.0
             && config.label_mask_out_unitig_fraction == 1.0
             && config.label_mask_other_unitig_fraction == 1.0) {
         logger->trace("Filtering by node");
         size_t total_nodes = masked_graph->num_nodes();
-        const auto &in_mask = dynamic_cast<const bitmap_vector&>(masked_graph->get_mask()).data();
+        const auto &in_mask
+            = static_cast<const bitmap_vector&>(masked_graph->get_mask()).data();
         sdsl::bit_vector mask = in_mask;
 
         // TODO: make this part multithreaded
         size_t kept_nodes = 0;
         call_ones(in_mask, [&](node_index node) {
-            uint64_t in_count = count_vector.first[node * 2];
-            uint64_t out_count = count_vector.first[node * 2 + 1];
+            uint64_t in_count = counts[node * 2];
+            uint64_t out_count = counts[node * 2 + 1];
 
             if (in_count >= min_label_in_count && out_count <= max_label_out_count) {
                 ++kept_nodes;
@@ -188,8 +324,6 @@ mask_nodes_by_label(const AnnotatedDBG &anno_graph,
     update_masked_graph_by_unitig(*masked_graph,
         [&](const std::string &, const std::vector<node_index> &path) -> Intervals {
             // return a set of intervals to keep in the graph
-            const auto &counts = count_vector.first;
-
             size_t in_kmer_count = 0;
 
             size_t begin = path.size();
@@ -254,7 +388,6 @@ make_initial_masked_graph(std::shared_ptr<const DeBruijnGraph> graph_ptr,
     // counts interleaved
     assert(counts.size() == mask.size() * 2);
 
-    add_complement |= graph_ptr->get_mode() == DeBruijnGraph::CANONICAL;
     auto masked_graph = std::make_shared<MaskedDeBruijnGraph>(
         graph_ptr,
         add_complement
@@ -298,57 +431,154 @@ make_initial_masked_graph(std::shared_ptr<const DeBruijnGraph> graph_ptr,
 }
 
 std::pair<sdsl::int_vector<>, sdsl::bit_vector>
-construct_diff_label_count_vector(const AnnotatedDBG &anno_graph,
+construct_diff_label_count_vector(const std::function<void(const std::function<void(const Label &, const bitmap&)>)> &generate_columns,
                                   const tsl::hopscotch_set<Label> &labels_in,
                                   const tsl::hopscotch_set<Label> &labels_out,
+                                  size_t max_index,
                                   size_t num_labels,
-                                  size_t num_threads) {
-    size_t width = sdsl::bits::hi(num_labels) + 1;
-    sdsl::bit_vector indicator(anno_graph.get_graph().max_index() + 1, false);
+                                  bool parallel,
+                                  bool add_out_labels_to_mask) {
+    logger->trace("Allocating mask vector");
+    size_t width = (sdsl::bits::hi(num_labels) + 1) * (1 + add_out_labels_to_mask);
+    sdsl::bit_vector indicator(max_index + 1, false);
 
+    logger->trace("Allocating count vector");
     // the in and out counts are stored interleaved
     sdsl::int_vector<> counts = aligned_int_vector(indicator.size() * 2, 0, width, 16);
+    logger->trace("done");
+
+    std::mutex vector_backup_mutex;
+
+    logger->trace("Populating count vector");
+    std::atomic_thread_fence(std::memory_order_release);
+    generate_columns([&](const Label &label, const bitmap &column) {
+        uint8_t col_indicator = static_cast<bool>(labels_in.count(label));
+        if (labels_out.count(label))
+            col_indicator |= 2;
+
+        auto add_in = [&indicator,&counts,&vector_backup_mutex,parallel](Column i) {
+            i = AnnotatedDBG::anno_to_graph_index(i);
+            assert(i != DeBruijnGraph::npos);
+            set_bit(indicator.data(), i, parallel, MO_RELAXED);
+            atomic_fetch_and_add(counts, i * 2, 1, vector_backup_mutex, MO_RELAXED);
+        };
+
+        auto add_out = [&indicator,&counts,&vector_backup_mutex,parallel,add_out_labels_to_mask](Column i) {
+            i = AnnotatedDBG::anno_to_graph_index(i);
+            assert(i != DeBruijnGraph::npos);
+            if (add_out_labels_to_mask)
+                set_bit(indicator.data(), i, parallel, MO_RELAXED);
+
+            atomic_fetch_and_add(counts, i * 2 + 1, 1, vector_backup_mutex, MO_RELAXED);
+        };
+
+        auto add_both = [&indicator,&counts,&vector_backup_mutex,parallel](Column i) {
+            i = AnnotatedDBG::anno_to_graph_index(i);
+            assert(i != DeBruijnGraph::npos);
+            set_bit(indicator.data(), i, parallel, MO_RELAXED);
+            atomic_fetch_and_add(counts, i * 2, 1, vector_backup_mutex, MO_RELAXED);
+            atomic_fetch_and_add(counts, i * 2 + 1, 1, vector_backup_mutex, MO_RELAXED);
+        };
+
+        switch (col_indicator) {
+            case 1: { column.call_ones(add_in); } break;
+            case 2: { column.call_ones(add_out); } break;
+            case 3: { column.call_ones(add_both); } break;
+            default: {}
+        }
+    });
+
+    std::atomic_thread_fence(std::memory_order_acquire);
+    logger->trace("done");
 
-    const auto &label_encoder = anno_graph.get_annotator().get_label_encoder();
-    const auto &binmat = anno_graph.get_annotator().get_matrix();
+    return std::make_pair(std::move(counts), std::move(indicator));
+}
 
-    tsl::hopscotch_map<uint64_t, uint8_t> code_to_indicator;
-    for (const std::string &label_in : labels_in) {
-        code_to_indicator[label_encoder.encode(label_in)] = 1;
+
+std::tuple<sdsl::int_vector<>, sdsl::bit_vector, sdsl::bit_vector>
+construct_diff_label_count_vector(const Annotator &annotator,
+                                  const tsl::hopscotch_set<Label> &labels_in,
+                                  const tsl::hopscotch_set<Label> &labels_out,
+                                  size_t max_index,
+                                  size_t num_labels,
+                                  size_t num_threads,
+                                  bool add_out_labels_to_mask) {
+    std::vector<Column> columns;
+    const auto &encoder = annotator.get_label_encoder();
+    AnnotatedDBG::Annotator::VLabels labels;
+    labels.reserve(labels_in.size() + labels_out.size());
+    columns.reserve(labels_in.size() + labels_out.size());
+    for (const auto &label : labels_in) {
+        labels.emplace_back(label);
+        columns.emplace_back(encoder.encode(label));
     }
-    for (const std::string &label_out : labels_out) {
-        code_to_indicator[label_encoder.encode(label_out)] |= 2;
+    for (const auto &label : labels_out) {
+        if (!labels_in.count(label)) {
+            labels.emplace_back(label);
+            columns.emplace_back(encoder.encode(label));
+        }
     }
 
-    std::vector<uint64_t> label_codes;
-    label_codes.reserve(code_to_indicator.size());
-    for (const auto &[code, indicator] : code_to_indicator) {
-        label_codes.push_back(code);
+    sdsl::bit_vector other_labels(encoder.size(), true);
+    for (Column column : columns) {
+        other_labels[column] = false;
     }
 
-    std::mutex vector_backup_mutex;
-    std::atomic_thread_fence(std::memory_order_release);
-    bool parallel = num_threads > 1;
+    auto [counts, init_mask] = construct_diff_label_count_vector([&](const auto &callback) {
+        annotator.get_matrix().call_columns(columns,
+            [&](size_t j, const bitmap_generator &column) {
+                callback(labels[j], column);
+            }, num_threads
+        );
+    }, labels_in, labels_out, max_index, num_labels, num_threads > 1, add_out_labels_to_mask);
 
-    binmat.call_columns(label_codes,
-        [&](auto col_idx, const bitmap &rows) {
-            uint8_t col_indicator = code_to_indicator[label_codes[col_idx]];
-            rows.call_ones([&](auto r) {
-                node_index i = AnnotatedDBG::anno_to_graph_index(r);
-                set_bit(indicator.data(), i, parallel, MO_RELAXED);
-                if (col_indicator & 1)
-                    atomic_fetch_and_add(counts, i * 2, 1, vector_backup_mutex, MO_RELAXED);
+    return std::make_tuple(std::move(counts), std::move(init_mask), std::move(other_labels));
+}
 
-                if (col_indicator & 2)
-                    atomic_fetch_and_add(counts, i * 2 + 1, 1, vector_backup_mutex, MO_RELAXED);
-            });
-        },
-        num_threads
-    );
+std::tuple<sdsl::int_vector<>, sdsl::bit_vector, sdsl::bit_vector, std::vector<size_t>>
+construct_diff_label_count_vector(const std::vector<std::string> &files,
+                                  const tsl::hopscotch_set<Label> &labels_in,
+                                  const tsl::hopscotch_set<Label> &labels_out,
+                                  size_t max_index,
+                                  size_t num_labels,
+                                  size_t num_threads,
+                                  bool add_out_labels_to_mask) {
+    std::vector<bool> other_labels_vec;
+    other_labels_vec.reserve(files.size());
+
+    std::vector<size_t> num_labels_per_file;
+    num_labels_per_file.reserve(files.size());
+
+    for (size_t i = 0; i < files.size(); ++i) {
+        auto label_encoder = annot::ColumnCompressed<>::read_label_encoder(files[i]);
+        num_labels_per_file.emplace_back(label_encoder.size());
+        for (const auto &label : label_encoder.get_labels()) {
+            other_labels_vec.emplace_back(!labels_in.count(label) && !labels_out.count(label));
+        }
+    }
 
-    std::atomic_thread_fence(std::memory_order_acquire);
+    auto other_labels = to_sdsl(std::move(other_labels_vec));
+
+    auto [counts, init_mask] = construct_diff_label_count_vector([&](const auto &callback) {
+        size_t j = 0;
+        for (size_t i = 0; i < files.size(); ++i) {
+            if (count_ones(other_labels, j, j + num_labels_per_file[i])
+                    != num_labels_per_file[i]) {
+                annot::ColumnCompressed<>::merge_load({ files[i] },
+                    [&](size_t offset, const auto &label, auto&& column) {
+                        if (!other_labels[j + offset])
+                            callback(label, *column);
+                    },
+                    num_threads
+                );
+            }
 
-    return std::make_pair(std::move(counts), std::move(indicator));
+            j += num_labels_per_file[i];
+        }
+    }, labels_in, labels_out, max_index, num_labels, num_threads > 1, add_out_labels_to_mask);
+
+    return std::make_tuple(std::move(counts), std::move(init_mask),
+                           std::move(other_labels), std::move(num_labels_per_file));
 }
 
 
@@ -360,11 +590,12 @@ void update_masked_graph_by_unitig(MaskedDeBruijnGraph &masked_graph,
     std::atomic<uint64_t> num_kept_nodes(0);
     bool parallel = num_threads > 1;
 
-    sdsl::bit_vector mask = dynamic_cast<const bitmap_vector&>(masked_graph.get_mask()).data();
+    sdsl::bit_vector mask = static_cast<const bitmap_vector&>(masked_graph.get_mask()).data();
 
     std::atomic_thread_fence(std::memory_order_release);
 
-    masked_graph.call_unitigs([&](const std::string &unitig, const std::vector<node_index> &path) {
+    masked_graph.call_unitigs([&](const std::string &unitig,
+                                  const std::vector<node_index> &path) {
         total_unitigs.fetch_add(1, MO_RELAXED);
 
         size_t last = 0;
diff --git a/metagraph/src/graph/annotated_graph_algorithm.hpp b/metagraph/src/graph/annotated_graph_algorithm.hpp
index 681caa864a..52cad9ec03 100644
--- a/metagraph/src/graph/annotated_graph_algorithm.hpp
+++ b/metagraph/src/graph/annotated_graph_algorithm.hpp
@@ -1,8 +1,6 @@
 #ifndef __ANNOTATED_GRAPH_ALGORITHM_HPP__
 #define __ANNOTATED_GRAPH_ALGORITHM_HPP__
 
-#include <vector>
-
 #include <tsl/hopscotch_set.h>
 
 #include "graph/annotated_dbg.hpp"
@@ -57,7 +55,17 @@ mask_nodes_by_label(const AnnotatedDBG &anno_graph,
                     const tsl::hopscotch_set<typename AnnotatedDBG::Annotator::Label> &labels_in_round2,
                     const tsl::hopscotch_set<typename AnnotatedDBG::Annotator::Label> &labels_out_round2,
                     const DifferentialAssemblyConfig &config,
-                    size_t num_threads = 1);
+                    size_t num_threads = 1,
+                    size_t num_parallel_files = std::numeric_limits<size_t>::max());
+
+std::shared_ptr<MaskedDeBruijnGraph>
+mask_nodes_by_label(std::shared_ptr<const DeBruijnGraph> graph_ptr,
+                    const std::vector<std::string> &files,
+                    const tsl::hopscotch_set<typename AnnotatedDBG::Annotator::Label> &labels_in,
+                    const tsl::hopscotch_set<typename AnnotatedDBG::Annotator::Label> &labels_out,
+                    const DifferentialAssemblyConfig &config,
+                    size_t num_threads = 1,
+                    size_t num_parallel_files = std::numeric_limits<size_t>::max());
 
 } // namespace graph
 } // namespace mtg