Merge branch 'release/0.6.2-beta'

CMakeLists.txt

@@ -8,7 +8,7 @@ project(octopus)
 
 set(octopus_VERSION_MAJOR 0)
 set(octopus_VERSION_MINOR 6)
-set(octopus_VERSION_PATCH 1)
+set(octopus_VERSION_PATCH 2)
 set(octopus_VERSION_RELEASE beta)
 
 # Get the current working branch

src/core/callers/cancer_caller.cpp

@@ -142,9 +142,11 @@ CancerCaller::infer_latents(const std::vector<Haplotype>& haplotypes,
     if (debug_log_) stream(*debug_log_) << "There are " << result->germline_genotypes_.size() << " candidate germline genotypes";
     evaluate_germline_model(*result, haplotype_likelihoods);
     evaluate_cnv_model(*result, haplotype_likelihoods);
-    fit_somatic_model(*result, haplotype_likelihoods);
-    evaluate_noise_model(*result, haplotype_likelihoods);
-    set_model_posteriors(*result);
+    if (haplotypes.size() > 1) {
+        fit_somatic_model(*result, haplotype_likelihoods);
+        evaluate_noise_model(*result, haplotype_likelihoods);
+        set_model_posteriors(*result);
+    }
     return result;
 }
 
@@ -410,8 +412,9 @@ void filter_with_germline_model(std::vector<CancerGenotype<Haplotype>>& genotype
 
 void CancerCaller::generate_cancer_genotypes(Latents& latents, const HaplotypeLikelihoodArray& haplotype_likelihoods) const
 {
-    const auto& germline_genotypes = latents.germline_genotypes_;
     const auto num_haplotypes = latents.haplotypes_.get().size();
+    if (num_haplotypes == 1) return;
+    const auto& germline_genotypes = latents.germline_genotypes_;
     const auto num_germline_genotypes = germline_genotypes.size();
     const auto max_possible_cancer_genotypes = num_haplotypes * num_germline_genotypes;
     const auto max_allowed_cancer_genotypes = std::max(parameters_.max_genotypes, num_germline_genotypes);
@@ -1547,15 +1550,17 @@ void CancerCaller::Latents::compute_haplotype_posteriors() const
             result.at(haplotype) += model_posteriors_.cnv * p.get<1>();
         }
     }
-    const auto credible_frequency = parameters_.get().min_expected_somatic_frequency;
-    const auto conditional_somatic_prob = compute_credible_somatic_mass(somatic_model_inferences_.posteriors.alphas, somatic_ploidy_, credible_frequency);
-    // Contribution from somatic model
-    for (const auto& p : zip(cancer_genotypes_, somatic_model_inferences_.posteriors.genotype_probabilities)) {
-        for (const auto& haplotype : p.get<0>().germline().copy_unique_ref()) {
-            result.at(haplotype) += model_posteriors_.somatic * p.get<1>();
-        }
-        for (const auto& haplotype : p.get<0>().somatic().copy_unique_ref()) {
-            result.at(haplotype) += model_posteriors_.somatic * conditional_somatic_prob * p.get<1>();
+    if (!cancer_genotypes_.empty()) {
+        const auto credible_frequency = parameters_.get().min_expected_somatic_frequency;
+        const auto conditional_somatic_prob = compute_credible_somatic_mass(somatic_model_inferences_.posteriors.alphas, somatic_ploidy_, credible_frequency);
+        // Contribution from somatic model
+        for (const auto& p : zip(cancer_genotypes_, somatic_model_inferences_.posteriors.genotype_probabilities)) {
+            for (const auto& haplotype : p.get<0>().germline().copy_unique_ref()) {
+                result.at(haplotype) += model_posteriors_.somatic * p.get<1>();
+            }
+            for (const auto& haplotype : p.get<0>().somatic().copy_unique_ref()) {
+                result.at(haplotype) += model_posteriors_.somatic * conditional_somatic_prob * p.get<1>();
+            }
         }
     }
     haplotype_posteriors_ = std::make_shared<Latents::HaplotypeProbabilityMap>(std::move(result));

src/core/tools/hapgen/haplotype_generator.cpp

@@ -912,11 +912,14 @@ void HaplotypeGenerator::populate_tree_with_novel_alleles()
     assert(!empty(novel_active_alleles));
     if (overlaps(novel_active_alleles.front(), active_region_)) {
         // Then there are duplicate insertions at the front of the novel region that must be removed
-        while (overlaps(novel_active_alleles.front(), active_region_)) {
+        while (!empty(novel_active_alleles) && overlaps(novel_active_alleles.front(), active_region_)) {
             novel_active_alleles.advance_begin(1);
-            assert(!empty(novel_active_alleles));
         }
-        novel_active_region = encompassing_region(novel_active_alleles);
+        if (!empty(novel_active_alleles)) {
+            novel_active_region = encompassing_region(novel_active_alleles);
+        } else {
+            novel_active_region = tail_region(*next_active_region_);
+        }
     }
     auto last_added_novel_itr = extend_tree_until(novel_active_alleles, tree_, policies_.haplotype_limits.holdout);
     if (last_added_novel_itr != std::cend(novel_active_alleles)) {

src/core/tools/hapgen/haplotype_tree.cpp

@@ -429,14 +429,7 @@ void HaplotypeTree::clear(const GenomicRegion& region)
     if (octopus::contains(region, tree_region)) {
         clear();
     } else if (overlaps(region, tree_region)) {
-        haplotype_leaf_cache_.clear();
-        std::list<Vertex> new_leafs {};
-        for (const Vertex leaf : haplotype_leafs_) {
-            const auto p = clear(leaf, contig_region(region));
-            if (p.second) new_leafs.push_back(p.first);
-        }
-        haplotype_leafs_ = new_leafs;
-        tree_region_ = boost::none;
+        clear_overlapped(region.contig_region());
     }
 }
 
@@ -498,6 +491,11 @@ HaplotypeTree::Vertex HaplotypeTree::get_previous_allele(const Vertex allele) co
     return *p.first;
 }
 
+bool HaplotypeTree::is_leaf(const Vertex v) const
+{
+    return boost::out_degree(v, tree_) == 0;
+}
+
 bool HaplotypeTree::is_bifurcating(const Vertex v) const
 {
     return boost::out_degree(v, tree_) > 1;
@@ -664,6 +662,20 @@ HaplotypeTree::find_equal_haplotype_leaf(const LeafIterator first, const LeafIte
                         });
 }
 
+void HaplotypeTree::clear_overlapped(const ContigRegion& region)
+{
+    haplotype_leaf_cache_.clear();
+    std::list<Vertex> new_leafs {};
+    for (const Vertex leaf : haplotype_leafs_) {
+        const auto p = clear(leaf, region);
+        if (p.second) new_leafs.push_back(p.first);
+    }
+    // As the tree is cleared, a  branch stub could be appended to a previous new leaf node
+    new_leafs.erase(std::remove_if(std::begin(new_leafs), std::end(new_leafs), [this] (Vertex v) { return !is_leaf(v); }), std::end(new_leafs));
+    haplotype_leafs_ = std::move(new_leafs);
+    tree_region_ = boost::none;
+}
+
 std::pair<HaplotypeTree::Vertex, bool>
 HaplotypeTree::clear(const Vertex leaf, const ContigRegion& region)
 {
@@ -677,7 +689,7 @@ HaplotypeTree::clear(const Vertex leaf, const ContigRegion& region)
 std::pair<HaplotypeTree::Vertex, bool>
 HaplotypeTree::clear_external(Vertex leaf, const ContigRegion& region)
 {
-    assert(boost::out_degree(leaf, tree_) == 0);
+    assert(is_leaf(leaf));
     while (leaf != root_) {
         if (boost::out_degree(leaf, tree_) > 0) {
             return std::make_pair(leaf, false);
@@ -694,6 +706,7 @@ HaplotypeTree::clear_external(Vertex leaf, const ContigRegion& region)
 std::pair<HaplotypeTree::Vertex, bool>
 HaplotypeTree::clear_internal(const Vertex leaf, const ContigRegion& region)
 {
+    assert(is_leaf(leaf));
     // TODO: we can optimise this for cases where region overlaps the leftmost alleles in the tree
     if (leaf == root_ || is_after(region, tree_[leaf])) {
         return std::make_pair(leaf, true);
@@ -748,7 +761,7 @@ HaplotypeTree::clear_internal(const Vertex leaf, const ContigRegion& region)
                                      });
         if (it == vertex_range.second) break;
         allele_to_move_to = *it; // i.e. move forward
-        if (boost::out_degree(allele_to_move, tree_) == 0) break;
+        if (is_leaf(allele_to_move)) break;
         // Safe to remove forward as we made this branch earlier via copies
         allele_to_move = remove_forward(allele_to_move);
     }

src/core/tools/hapgen/haplotype_tree.hpp

@@ -110,6 +110,7 @@ class HaplotypeTree
     using LeafIterator  = decltype(haplotype_leafs_)::const_iterator;
     using CacheIterator = decltype(haplotype_leaf_cache_)::iterator;
 
+    bool is_leaf(Vertex v) const;
     bool is_bifurcating(Vertex v) const;
     Vertex remove_forward(Vertex u);
     Vertex remove_backward(Vertex v);
@@ -126,6 +127,7 @@ class HaplotypeTree
                                            const Haplotype& haplotype) const;
     LeafIterator find_equal_haplotype_leaf(LeafIterator first, LeafIterator last,
                                            const Haplotype& haplotype) const;
+    void clear_overlapped(const ContigRegion& region);
     std::pair<Vertex, bool> clear(Vertex leaf, const ContigRegion& region);
     std::pair<Vertex, bool> clear_external(Vertex leaf, const ContigRegion& region);
     std::pair<Vertex, bool> clear_internal(Vertex leaf, const ContigRegion& region);