Consider ambiguous support for AlleleDepth and AlleleFrequency

Author: Daniel Cooke

src/core/csr/measures/allele_depth.cpp

@@ -46,40 +46,34 @@ bool is_evaluable(const VcfRecord& call, const VcfRecord::SampleName& sample)
     return has_called_alt_allele(call, sample);
 }
 
+template <typename T>
+void pop_front(std::vector<T>& v)
+{
+    v.erase(std::cbegin(v));
+}
+
+auto min_support_count(const AlleleSupportMap& support)
+{
+    const static auto support_less = [] (const auto& lhs, const auto& rhs) { return lhs.second.size() < rhs.second.size(); };
+    return std::min_element(std::cbegin(support), std::cend(support), support_less)->second.size();
+}
+
 } // namespace
 
 Measure::ResultType AlleleDepth::do_evaluate(const VcfRecord& call, const FacetMap& facets) const
 {
     const auto& samples = get_value<Samples>(facets.at("Samples"));
-    const auto& assignments = get_value<ReadAssignments>(facets.at("ReadAssignments")).support;
+    const auto& assignments = get_value<ReadAssignments>(facets.at("ReadAssignments"));
     std::vector<boost::optional<int>> result {};
     result.reserve(samples.size());
     for (const auto& sample : samples) {
         boost::optional<int> sample_result {};
         if (is_evaluable(call, sample)) {
-            const auto& sample_assignments = assignments.at(sample);
             std::vector<Allele> alleles; bool has_ref;
-            std::tie(alleles, has_ref) = get_called_alleles(call, sample, true);
-            assert(!alleles.empty());
-            std::vector<int> allele_counts(alleles.size());
-            for (const auto& p : sample_assignments) {
-                const auto& haplotype = p.first;
-                const auto& reads = p.second;
-                const auto haplotype_support_depth = count_overlapped(reads, call);
-                if (haplotype_support_depth > 0) {
-                    std::transform(std::cbegin(alleles), std::cend(alleles), std::cbegin(allele_counts), std::begin(allele_counts),
-                                   [&] (const auto& allele, auto count) {
-                                       if (haplotype.includes(allele)) {
-                                           count += haplotype_support_depth;
-                                       }
-                                       return count;
-                                   });
-                }
-            }
-            auto first_called_count_itr = std::cbegin(allele_counts);
-            if (has_ref) ++first_called_count_itr;
-            assert(first_called_count_itr != std::cend(allele_counts));
-            sample_result = *std::min_element(first_called_count_itr, std::cend(allele_counts));
+            std::tie(alleles, has_ref) = get_called_alleles(call, sample);
+            if (has_ref) pop_front(alleles); // ref always first
+            const auto allele_support = compute_allele_support(alleles, assignments, sample);
+            sample_result = min_support_count(allele_support);
         }
         result.push_back(sample_result);
     }

src/core/csr/measures/allele_frequency.cpp

@@ -8,13 +8,11 @@
 
 #include <boost/variant.hpp>
 
-#include "core/tools/read_assigner.hpp"
 #include "core/types/allele.hpp"
 #include "io/variant/vcf_record.hpp"
-#include "io/variant/vcf_spec.hpp"
-#include "utils/genotype_reader.hpp"
-#include "../facets/samples.hpp"
-#include "../facets/read_assignments.hpp"
+#include "utils/concat.hpp"
+#include "allele_depth.hpp"
+#include "depth.hpp"
 
 namespace octopus { namespace csr {
 
@@ -25,69 +23,19 @@ std::unique_ptr<Measure> AlleleFrequency::do_clone() const
     return std::make_unique<AlleleFrequency>(*this);
 }
 
-namespace {
-
-bool is_canonical(const VcfRecord::NucleotideSequence& allele) noexcept
-{
-    const static VcfRecord::NucleotideSequence deleted_allele {vcfspec::deletedBase};
-    return !(allele == vcfspec::missingValue || allele == deleted_allele);
-}
-
-bool has_called_alt_allele(const VcfRecord& call, const VcfRecord::SampleName& sample)
-{
-    if (!call.has_genotypes()) return true;
-    const auto& genotype = get_genotype(call, sample);
-    return std::any_of(std::cbegin(genotype), std::cend(genotype),
-                       [&] (const auto& allele) { return allele != call.ref() && is_canonical(allele); });
-}
-
-bool is_evaluable(const VcfRecord& call, const VcfRecord::SampleName& sample)
-{
-    return has_called_alt_allele(call, sample);
-}
-
-} // namespace
-
 Measure::ResultType AlleleFrequency::do_evaluate(const VcfRecord& call, const FacetMap& facets) const
 {
-    const auto& samples = get_value<Samples>(facets.at("Samples"));
-    const auto& assignments = get_value<ReadAssignments>(facets.at("ReadAssignments")).support;
-    std::vector<boost::optional<double>> result {};
-    result.reserve(samples.size());
-    for (const auto& sample : samples) {
-        boost::optional<double> sample_result {};
-        if (is_evaluable(call, sample)) {
-            const auto& sample_assignments = assignments.at(sample);
-            std::vector<Allele> alleles; bool has_ref;
-            std::tie(alleles, has_ref) = get_called_alleles(call, sample, true);
-            assert(!alleles.empty());
-            std::size_t read_count {0};
-            std::vector<unsigned> allele_counts(alleles.size());
-            for (const auto& p : sample_assignments) {
-                const auto& haplotype = p.first;
-                const auto& reads = p.second;
-                const auto haplotype_support_depth = count_overlapped(reads, call);
-                if (haplotype_support_depth > 0) {
-                    std::transform(std::cbegin(alleles), std::cend(alleles), std::cbegin(allele_counts), std::begin(allele_counts),
-                                   [&] (const auto& allele, auto count) {
-                                       if (haplotype.includes(allele)) {
-                                           count += haplotype_support_depth;
-                                       }
-                                       return count;
-                                   });
-                    read_count += haplotype_support_depth;
-                }
-            }
-            if (read_count > 0) {
-                auto first_called_count_itr = std::cbegin(allele_counts);
-                if (has_ref) ++first_called_count_itr;
-                assert(first_called_count_itr != std::cend(allele_counts));
-                const auto min_count_itr = std::min_element(first_called_count_itr, std::cend(allele_counts));
-                sample_result = static_cast<double>(*min_count_itr) / read_count;
-            }
+    const auto allele_depths = boost::get<std::vector<boost::optional<int>>>(AlleleDepth{}.evaluate(call, facets));
+    const auto depths = boost::get<std::vector<std::size_t>>(Depth{true, false}.evaluate(call, facets));
+    std::vector<boost::optional<double>> result(allele_depths.size());
+    std::transform(std::cbegin(allele_depths), std::cend(allele_depths), std::cbegin(depths), std::begin(result),
+                   [] (const auto allele_depth, const auto depth) -> boost::optional<double> {
+        if (allele_depth && depth > 0) {
+            return static_cast<double>(*allele_depth) / depth;
+        } else {
+            return boost::none;
         }
-        result.push_back(sample_result);
-    }
+    });
     return result;
 }
 
@@ -108,7 +56,7 @@ std::string AlleleFrequency::do_describe() const
 
 std::vector<std::string> AlleleFrequency::do_requirements() const
 {
-    return {"Samples", "ReadAssignments"};
+    return concat(AlleleDepth{}.requirements(), Depth{true, false}.requirements());
 }
 
 } // namespace csr