crispr_quality_control.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_CRISPR_QUALITY_CONTROL_HPP
#define SCRAN_QC_CRISPR_QUALITY_CONTROL_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <type_traits>
 
#include "tatami/tatami.hpp"
 
#include "find_median_mad.hpp"
#include "per_cell_qc_metrics.hpp"
#include "choose_filter_thresholds.hpp"
 
namespace scran_qc {
 
struct ComputeCrisprQcMetricsOptions {
    int num_threads = 1;
};
 
template<typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int>
struct ComputeCrisprQcMetricsBuffers {
    Sum_* sum;
 
    Detected_* detected;
 
    Value_* max_value;
 
    Index_* max_index;
};
 
template<typename Value_, typename Index_, typename Sum_, typename Detected_>
void compute_crispr_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& output,
    const ComputeCrisprQcMetricsOptions& options)
{
    PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_> tmp;
    tmp.sum = output.sum;
    tmp.detected = output.detected;
    tmp.max_value = output.max_value;
    tmp.max_index = output.max_index;
 
    PerCellQcMetricsOptions opt;
    opt.num_threads = options.num_threads;
    per_cell_qc_metrics(mat, std::vector<uint8_t*>{}, tmp, opt);
}
 
template<typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int>
struct ComputeCrisprQcMetricsResults {
    std::vector<Sum_> sum;
 
    std::vector<Detected_> detected;
 
    std::vector<Value_> max_value;
 
    std::vector<Index_> max_index;
};
 
template<typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int>
ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_> compute_crispr_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const ComputeCrisprQcMetricsOptions& options)
{
    auto NC = mat.ncol();
    ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_> x;
    ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_> output;
 
    output.sum.resize(NC);
    x.sum = output.sum.data();
 
    output.detected.resize(NC);
    x.detected = output.detected.data();
 
    output.max_value.resize(NC);
    x.max_value = output.max_value.data();
 
    output.max_index.resize(NC);
    x.max_index = output.max_index.data();
 
    compute_crispr_qc_metrics(mat, x, options);
    return output;
}
 
struct ComputeCrisprQcFiltersOptions {
    double max_value_num_mads = 3;
};
 
namespace internal {
 
template<typename Float_, class Host_, typename Sum_, typename Detected_, typename Value_, typename Index_, typename BlockSource_>
void crispr_populate(Host_& host, size_t n, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& res, BlockSource_ block, const ComputeCrisprQcFiltersOptions& options) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    auto buffer = [&]() {
        if constexpr(unblocked) {
            return std::vector<Float_>(n);
        } else {
            return FindMedianMadWorkspace<Float_, size_t>(n, block);
        }
    }();
 
    // Subsetting to the observations in the top 50% of proportions.
    static_assert(std::is_floating_point<Float_>::value);
    std::vector<Float_> maxprop;
    maxprop.reserve(n);
    for (size_t i = 0; i < n; ++i) {
        maxprop.push_back(static_cast<Float_>(res.max_value[i]) / static_cast<Float_>(res.sum[i]));
    }
 
    FindMedianMadOptions fopt;
    fopt.median_only = true;
    auto prop_res = [&]() {
        if constexpr(unblocked) {
            return find_median_mad(n, maxprop.data(), buffer.data(), fopt);
        } else {
            return find_median_mad_blocked(n, maxprop.data(), block, &buffer, fopt);
        }
    }();
 
    for (size_t i = 0; i < n; ++i) {
        auto limit = [&]() {
            if constexpr(unblocked){
                return prop_res.median;
            } else {
                return prop_res[block[i]].median;
            }
        }();
        if (maxprop[i] >= limit) {
            maxprop[i] = res.max_value[i];
        } else {
            maxprop[i] = std::numeric_limits<Float_>::quiet_NaN(); // ignored during threshold calculation.
        }
    }
 
    // Filtering on the max counts.
    ChooseFilterThresholdsOptions copt;
    copt.num_mads = options.max_value_num_mads;
    copt.log = true;
    copt.upper = false;
    host.get_max_value() = [&]() {
        if constexpr(unblocked) {
            return choose_filter_thresholds(n, maxprop.data(), buffer.data(), copt).lower;
        } else {
            return internal::strip_threshold<true>(choose_filter_thresholds_blocked(n, maxprop.data(), block, &buffer, copt));
        }
    }();
}
 
template<class Host_, typename Sum_, typename Detected_, typename Value_, typename Index_, typename BlockSource_, typename Output_>
void crispr_filter(const Host_& host, size_t n, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics, BlockSource_ block, Output_* output) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    std::fill_n(output, n, 1);
 
    const auto& mv = host.get_max_value();
    for (size_t i = 0; i < n; ++i) {
        auto thresh = [&]() {
            if constexpr(unblocked) {
                return mv;
            } else {
                return mv[block[i]];
            }
        }();
        output[i] = output[i] && (metrics.max_value[i] >= thresh);
    }
}
 
template<typename Sum_, typename Detected_, typename Value_, typename Index_>
ComputeCrisprQcMetricsBuffers<const Sum_, const Detected_, const Value_, const Index_> to_buffer(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics) {
    ComputeCrisprQcMetricsBuffers<const Sum_, const Detected_, const Value_, const Index_> buffer;
    buffer.sum = metrics.sum.data();
    buffer.detected = metrics.detected.data();
    buffer.max_value = metrics.max_value.data();
    buffer.max_index = metrics.max_index.data();
    return buffer;
}
 
}
template<typename Float_ = double>
class CrisprQcFilters {
public:
    Float_ get_max_value() const {
        return my_max_value;
    }
 
    Float_& get_max_value() {
        return my_max_value;
    }
 
private:
    Float_ my_max_value = 0;
 
public:
    template<typename Sum_, typename Detected_, typename Value_, typename Index_, typename Output_>
    void filter(size_t num, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics, Output_* output) const {
        internal::crispr_filter(*this, num, metrics, false, output);
    }
 
    template<typename Sum_, typename Detected_, typename Value_, typename Index_, typename Output_>
    void filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics, Output_* output) const {
        return filter(metrics.max_value.size(), internal::to_buffer(metrics), output);
    }
 
    template<typename Output_ = uint8_t, typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int>
    std::vector<Output_> filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics) const {
        std::vector<Output_> output(metrics.max_value.size());
        filter(metrics, output.data());
        return output;
    }
};
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Value_, typename Index_>
CrisprQcFilters<Float_> compute_crispr_qc_filters(
    size_t num,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics,
    const ComputeCrisprQcFiltersOptions& options)
{
    CrisprQcFilters<Float_> output;
    internal::crispr_populate<Float_>(output, num, metrics, false, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int>
CrisprQcFilters<Float_> compute_crispr_qc_filters(
    const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics,
    const ComputeCrisprQcFiltersOptions& options)
{
    return compute_crispr_qc_filters(metrics.max_value.size(), internal::to_buffer(metrics), options);
}
 
template<typename Float_ = double>
class CrisprQcBlockedFilters {
public:
    const std::vector<Float_>& get_max_value() const {
        return my_max_value;
    }
 
    std::vector<Float_>& get_max_value() {
        return my_max_value;
    }
 
private:
    std::vector<Float_> my_sum;
    std::vector<Float_> my_max_value;
 
public:
    template<typename Sum_, typename Detected_, typename Value_, typename Index_, typename Block_, typename Output_>
    void filter(size_t num, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics, const Block_* block, Output_* output) const {
        internal::crispr_filter(*this, num, metrics, block, output);
    }
 
    template<typename Sum_, typename Detected_, typename Value_, typename Index_, typename Block_, typename Output_>
    void filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics, const Block_* block, Output_* output) const {
        filter(metrics.max_value.size(), internal::to_buffer(metrics), block, output);
    }
 
    template<typename Output_ = uint8_t, typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int, typename Block_ = int>
    std::vector<Output_> filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics, const Block_* block) const {
        std::vector<Output_> output(metrics.max_value.size());
        filter(metrics, block, output.data());
        return output;
    }
};
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Value_, typename Index_, typename Block_>
CrisprQcBlockedFilters<Float_> compute_crispr_qc_filters_blocked(
    size_t num,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics,
    const Block_* block,
    const ComputeCrisprQcFiltersOptions& options)
{
    CrisprQcBlockedFilters<Float_> output;
    internal::crispr_populate<Float_>(output, num, metrics, block, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Value_, typename Index_, typename Block_>
CrisprQcBlockedFilters<Float_> compute_crispr_qc_filters_blocked(
    const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics,
    const Block_* block,
    const ComputeCrisprQcFiltersOptions& options)
{
    return compute_crispr_qc_filters_blocked(metrics.max_value.size(), internal::to_buffer(metrics), block, options);
}
 
}
 
#endif