Single command haplotype sampling

Brewfile

@@ -22,3 +22,5 @@ brew "boost"
 brew "pybind11"
 brew "pandoc"
 brew "openssl"
+tap "brewsci/bio"
+brew "brewsci/bio/kmc"

Dockerfile

@@ -44,7 +44,7 @@ RUN apt-get -qq -y update && apt-get -qq -y upgrade && apt-get -y install \
     samtools curl unzip redland-utils librdf-dev cmake pkg-config wget gtk-doc-tools \
     raptor2-utils rasqal-utils bison flex gawk libgoogle-perftools-dev liblz4-dev liblzma-dev \
     libcairo2-dev libpixman-1-dev libffi-dev libcairo-dev libprotobuf-dev libboost-all-dev \
-    tabix bcftools libzstd-dev pybind11-dev python3-pybind11 pandoc libssl-dev
+    tabix bcftools libzstd-dev pybind11-dev python3-pybind11 pandoc libssl-dev kmc
 ###DEPS_END###
 
 FROM packages AS build

README.md

@@ -91,7 +91,7 @@ On other distros, or if you do not have root access, you will need to perform th
                          automake gettext autopoint libtool jq bsdmainutils bc rs parallel \
                          npm curl unzip redland-utils librdf-dev bison flex gawk lzma-dev \
                          liblzma-dev liblz4-dev libffi-dev libcairo-dev libboost-all-dev \
-                         libzstd-dev pybind11-dev python3-pybind11 libssl-dev
+                         libzstd-dev pybind11-dev python3-pybind11 libssl-dev kmc
 
 At present, you will need GCC version 9 or greater, with support for C++17, to compile vg. (Check your version with `gcc --version`.) GCC up to 11.4.0 is supported.
 

doc/test-docs.sh

@@ -16,6 +16,9 @@ cd "${HERE}/../test"
 echo txm --jobs 1 "${HERE}/../README.md"
 txm --jobs 1 "${HERE}/../README.md"
 
+# Test extra test READMEs
+(cd haplotype-sampling && txm --jobs 1 "README.md")
+
 # Run all the wiki tests
 find "${HERE}/wiki" -name "*.md" | xargs -n 1 -t txm --jobs 1
 

src/index_registry.hpp

@@ -12,6 +12,8 @@
 #include <stdexcept>
 #include <limits>
 
+#include <gbwtgraph/utils.h>
+
 namespace vg {
 
 using namespace std;
@@ -36,6 +38,9 @@ using IndexGroup = set<IndexName>;
 
 /**
  * Names a recipe in the collection of registered recipes.
+ *
+ * A max-value number means the indexes are provided and not generated by any
+ * recipe.
  */
 using RecipeName = pair<IndexGroup, size_t>;
 
@@ -128,6 +133,16 @@ struct IndexingParameters {
     static double thread_chunk_inflation_factor;
     // whether indexing algorithms will log progress (if available) [Basic]
     static Verbosity verbosity;
+    // Set of samples to make references and retain during haplotype sampling [{}]
+    static gbwtgraph::sample_name_set haplotype_sampling_reference_samples;
+    // Number of haplotypes to sample during haplotype sampling [Recombinator::NUM_CANDIDATES]
+    static size_t haplotype_sampling_num_haplotypes;
+    // Whether to do diploid sampling during haplotype sampling [true]
+    static bool haplotype_sampling_diploid;
+    // length of k-mer used in haplotype index and kmer counting [29]
+    static int haplotype_sampling_minimizer_k;
+    // length of window used in haplotype index [11]
+    static int haplotype_sampling_minimizer_w;
 };
 
 /**
@@ -167,6 +182,13 @@ class IndexingPlan {
     string output_filepath(const IndexName& identifier) const;
 
     /// Get the suffix with which to save the given index's files.
+    ///
+    /// The path will substitute scopes into any {wildcards} in the index's
+    /// suffix, and scopes not used in the suffix will appear in the
+    /// filename in key order.
+    ///
+    /// The suffix used will be the first one for which all wildcards can be
+    /// substituted. 
     string output_filepath(const IndexName& identifier, size_t chunk, size_t num_chunks) const;
 
     /// Ge the steps of the plan
@@ -175,24 +197,43 @@ class IndexingPlan {
     /// Returns true if the given index is to be intermediate under the given
     /// plan, and false if it is to be preserved.
     bool is_intermediate(const IndexName& identifier) const;
+
+    /// Get the scopes that the given index will be qualified with when
+    /// generated. Only works for indexes that will be generated by the plan,
+    /// not for inputs.
+    const map<string, string> get_scopes(const IndexName& identifier) const;
+
+    /// Get the scopes that the given index will be qualified with when
+    /// generated. Only works for indexes that will be generated by the plan,
+    /// not for inputs.
+    const map<string, string>& get_scopes(const IndexName& identifier);
 
     // TODO: is this where this function wants to live?
     /// The memory limit, with a little slosh for prediction inaccuracy
     int64_t target_memory_usage() const;
-    /// The mmeory limit with no slosh
+    /// The memory limit with no slosh
     int64_t literal_target_memory_usage() const;
 
     /// Returns the recipes in the plan that depend on this index, including the one in which
-    /// it was created (if any)
+    /// it was created (if any).
     set<RecipeName> dependents(const IndexName& identifier) const;
 
 protected:
 
+    /// Add a scope to the scopes the given index will be qualified with when
+    /// generated, if not present already.
+    void add_scope(const IndexName& identifier, const string& key, const string& value);
+
     /// The steps to be invoked in the plan. May be empty before the plan is
     /// actually planned.
     vector<RecipeName> steps;
     /// The indexes to create as outputs.
     set<IndexName> targets;
+    /// The scopes qualifying each index in the plan (for example, if an index
+    /// applies only to a particular sample, it will be scoped to the sample
+    /// name). Stored using an ordered map to ensure a consistent iteration
+    /// order.
+    map<IndexName, map<string, string>> scopes;
 
     /// The registry that the plan is using.
     /// The registry must not move while the plan is in use.
@@ -244,6 +285,17 @@ class IndexRegistry {
 
     /// Register an index containing the given identifier
     void register_index(const IndexName& identifier, const string& suffix);
+
+    /// Register an index containing the given identifier, with multiple possible suffixes.
+    /// The first suffix where all {wildcards} can be substituted with scopes will be used.
+    void register_index(const IndexName& identifier, const vector<string>& suffixes);
+
+    /// Get the names of all brace-enclosed {wildcards} in the given pattern
+    static set<string> get_wildcards(const string& pattern);
+
+    /// Substitute wildcards into the given pattern. All wildcards must have
+    /// values assigned. Extra values not used are allowed.
+    static string substitute_wildcards(const string& pattern, const map<string, string> values);
 
     /// Register a recipe to produce an index using other indexes
     /// or input files. Recipes registered earlier will have higher priority.
@@ -254,12 +306,22 @@ class IndexRegistry {
     /// Indicate one recipe is a broadened version of another. The indexes consumed and produced
     /// by the generalization must be semantically identical to those of the generalizee
     void register_generalization(const RecipeName& generalizer, const RecipeName& generalizee);
-
-    /// Indicate a serialized file that contains some identified index
-    void provide(const IndexName& identifier, const string& filename);
-
-    /// Indicate a list of serialized files that contains some identified index
-    void provide(const IndexName& identifier, const vector<string>& filenames);
+
+    /// Indicate a serialized file that contains some identified index,
+    /// optionally with scopes that propagates to descendant files.
+    ///
+    /// TODO: If scopes contain ".", we can run into problems with combinations
+    /// of different scopes producing the same final string. Right now we only
+    /// use one kind of scope, which avoids this.
+    void provide(const IndexName& identifier, const string& filename, const map<string, string>& scopes = {});
+
+    /// Indicate a list of serialized files that contains some identified index,
+    /// optionally with scopes that propagates to descendant files.
+    ///
+    /// TODO: If scopes contain ".", we can run into problems with combinations
+    /// of different scopes producing the same final string. Right now we only
+    /// use one kind of scope, which avoids this.
+    void provide(const IndexName& identifier, const vector<string>& filenames, const map<string, string>& scopes = {});
 
     /// Remove a provided index
     void reset(const IndexName& identifier);
@@ -269,6 +331,7 @@ class IndexRegistry {
     bool available(const IndexName& identifier) const;
 
     /// Get the possible filename(s) associated with the given index with the given prefix.
+    /// TODO: Get this to account for sample-scoped indexes.
     vector<string> get_possible_filenames(const IndexName& identifier) const;
 
     /// Get the filename(s) associated with the given index. Aborts if the
@@ -317,6 +380,14 @@ class IndexRegistry {
     /// generate a plan to create the indexes
     IndexingPlan make_plan(const IndexGroup& end_products) const;
 
+    /// Check if a recipe identifier correesponds to a recipe.
+    ///
+    /// Recipe identifiers not corresponding to actual recipes are used during
+    /// planning to represent provided inputs.
+    ///
+    /// TODO: Refactor that with some kind of tagged union or optional.
+    bool has_recipe(const RecipeName& recipe_name) const;
+
     /// use a recipe identifier to get the recipe
     const IndexRecipe& get_recipe(const RecipeName& recipe_name) const;
 
@@ -371,12 +442,15 @@ class IndexFile {
 
     /// Create a new IndexFile with a unique identifier
     IndexFile(const IndexName& identifier, const string& suffix);
+
+    /// Create a new IndexFile with a unique identifier, which may use one of several suffixes.
+    IndexFile(const IndexName& identifier, const vector<string>& suffixes);
 
     /// Get the globally unique identifier for this index
     const IndexName& get_identifier() const;
 
-    /// Returns the suffix to be used for this index
-    const string& get_suffix() const;
+    /// Returns the suffixes that can used for this index
+    const vector<string>& get_suffixes() const;
 
     /// Get the filename(s) that contain this index
     const vector<string>& get_filenames() const;
@@ -386,7 +460,18 @@ class IndexFile {
 
     /// Assign constructed filenames to this index
     void assign_constructed(const vector<string>& filenames);
-
+
+    /// Add a new scope to qualify the files in the index, if the scope is
+    /// not already on the index.
+    ///
+    /// Should only be used on files that are actually filled in already;
+    /// scopes for files that don't exist yet are the responsibility of the
+    /// IndexingPlan.
+    void add_scope(const string& key, const string& value);
+
+    /// Get all scopes qualifying the index (such as a sample name).
+    const map<string, string>& get_scopes() const;
+
     /// Returns true if the index has already been built or provided
     bool is_finished() const;
 
@@ -401,11 +486,14 @@ class IndexFile {
     // the global identifier for the
     IndexName identifier;
 
-    // the suffix it adds to output files
-    const string suffix;
+    // the suffixes it can add to output files
+    const vector<string> suffixes;
 
     // the filename(s) associated with the index
     vector<string> filenames;
+
+    // The scopes qualifying the index.
+    map<string, string> scopes;
 
     // keep track of whether the index was provided directly
     bool provided_directly = false;

src/subcommand/autoindex_main.cpp

@@ -390,7 +390,16 @@ int main_autoindex(int argc, char** argv) {
             if (!registry.available(target)) {
                 vector<string> inferred_file_names = registry.get_possible_filenames(target);
                 for (const string& filename : inferred_file_names) {
-                    if (ifstream(filename).is_open()) {
+                    if (target == "Giraffe GBZ" && !ends_with(filename, ".giraffe.gbz")) {
+                        // TODO: Giraffe GBZ indexes can be saved as .<sample>.gbz or .giraffe.gbz.
+                        // But we can't pick up .<sample>.gbz automatically without possibly mistaking a plain GBZ for a Giraffe one.
+                        // And we don't handle haplotype samplign in autoindex yet anyway.
+                        // So only find Giraffe GBZs.
+                        // TODO: Allow the siffixes to have Snakemake-style
+                        // wildcards that populate scopes on the files.
+                        continue;
+                    }
+                    if (file_exists(filename)) {
                         logger.info() << "Guessing that " << filename << " is " << target << endl;
                         registry.provide(target, filename);
                         break;