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 <cstddef>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
#include "per_cell_qc_metrics.hpp"
#include "choose_filter_thresholds.hpp"
#include "utils.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<const unsigned char*>{}, 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)
{
    const auto NC = mat.ncol();
    ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_> x;
    ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_> output;
 
    tatami::resize_container_to_Index_size(output.sum, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.sum = output.sum.data();
 
    tatami::resize_container_to_Index_size(output.detected, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.detected = output.detected.data();
 
    tatami::resize_container_to_Index_size(output.max_value, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.max_value = output.max_value.data();
 
    tatami::resize_container_to_Index_size(output.max_index, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.max_index = output.max_index.data();
 
    compute_crispr_qc_metrics(mat, x, options);
    return output;
}
 
struct ComputeCrisprQcFiltersOptions {
    double max_value_num_mads = 3;
};
 
template<typename Float_, class Host_, typename Sum_, typename Detected_, typename Value_, typename Index_, typename BlockSource_>
void compute_crispr_qc_filters_internal(
    Host_& host,
    const std::size_t num_cells,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& res,
    BlockSource_ block,
    const std::size_t num_blocks,
    const ComputeCrisprQcFiltersOptions& options
) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    auto buffer = [&]{
        if constexpr(unblocked) {
            return sanisizer::create<std::vector<Float_> >(num_cells);
        } else {
            return ChooseFilterThresholdsBlockedWorkspace<Float_>(num_cells, block, num_blocks);
        }
    }();
 
    // Subsetting to the observations in the top 50% of proportions.
    static_assert(std::is_floating_point<Float_>::value);
    std::vector<Float_> maxprop;
    maxprop.reserve(num_cells);
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
        maxprop.push_back(static_cast<Float_>(res.max_value[i]) / static_cast<Float_>(res.sum[i]));
    }
 
    auto prop_res = [&]{
        if constexpr(unblocked) {
            std::copy_n(maxprop.begin(), num_cells, buffer.begin());
            return quickstats::median<Float_>(num_cells, buffer.data());
        } else {
            std::vector<Float_> output;
            output.reserve(num_blocks);
            process_blocks_for_choose_filter_thresholds(
                num_cells,
                maxprop.data(),
                block,
                num_blocks,
                buffer,
                [&](const std::size_t len, Float_* const ptr) -> void {
                    output.push_back(quickstats::median<Float_>(len, ptr));
                }
            );
            return output;
        }
    }();
 
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
        auto limit = [&]{
            if constexpr(unblocked){
                return prop_res;
            } else {
                return prop_res[block[i]];
            }
        }();
        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(num_cells, maxprop.data(), buffer.data(), copt).lower;
        } else {
            return extract_filter_thresholds<true>(choose_filter_thresholds_blocked(num_cells, maxprop.data(), block, num_blocks, buffer, copt));
        }
    }();
}
 
template<class Host_, typename Sum_, typename Detected_, typename Value_, typename Index_, typename BlockSource_, typename Output_>
void apply_crispr_qc_filters_internal(
    const Host_& host,
    const std::size_t n,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics,
    BlockSource_ block,
    Output_* const output
) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    std::fill_n(output, n, 1);
 
    const auto& mv = host.get_max_value();
    for (I<decltype(n)> 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_> crispr_qc_results_to_buffers(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(const std::size_t num_cells, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics, Output_* const output) const {
        apply_crispr_qc_filters_internal(*this, num_cells, metrics, false, output);
    }
 
    template<typename Sum_, typename Detected_, typename Value_, typename Index_, typename Output_>
    void filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics, Output_* const output) const {
        return filter(metrics.max_value.size(), crispr_qc_results_to_buffers(metrics), output);
    }
 
    template<typename Output_ = unsigned char, typename Sum_, typename Detected_, typename Value_, typename Index_>
    std::vector<Output_> filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics) const {
        auto output = sanisizer::create<std::vector<Output_> >(metrics.max_value.size()
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        filter(metrics, output.data());
        return output;
    }
};
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Value_, typename Index_>
CrisprQcFilters<Float_> compute_crispr_qc_filters(
    const std::size_t num_cells,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics,
    const ComputeCrisprQcFiltersOptions& options)
{
    CrisprQcFilters<Float_> output;
    compute_crispr_qc_filters_internal<Float_>(output, num_cells, metrics, false, 0, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Value_, typename Index_>
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(), crispr_qc_results_to_buffers(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(const std::size_t num_cells, const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics, const Block_* const block, Output_* const output) const {
        apply_crispr_qc_filters_internal(*this, num_cells, 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_* const block, Output_* const output) const {
        filter(metrics.max_value.size(), crispr_qc_results_to_buffers(metrics), block, output);
    }
 
    template<typename Output_ = unsigned char, typename Sum_, typename Detected_, typename Value_, typename Index_, typename Block_>
    std::vector<Output_> filter(const ComputeCrisprQcMetricsResults<Sum_, Detected_, Value_, Index_>& metrics, const Block_* const block) const {
        auto output = sanisizer::create<std::vector<Output_> >(metrics.max_value.size()
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        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(
    const std::size_t num_cells,
    const ComputeCrisprQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& metrics,
    const Block_* const block,
    const std::size_t num_blocks,
    const ComputeCrisprQcFiltersOptions& options
) {
    CrisprQcBlockedFilters<Float_> output;
    compute_crispr_qc_filters_internal<Float_>(output, num_cells, metrics, block, num_blocks, 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_* const block,
    const std::size_t num_blocks,
    const ComputeCrisprQcFiltersOptions& options
) {
    return compute_crispr_qc_filters_blocked(metrics.max_value.size(), crispr_qc_results_to_buffers(metrics), block, num_blocks, options);
}
 
}
 
#endif