rna_quality_control.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_RNA_QUALITY_CONTROL_HPP
#define SCRAN_QC_RNA_QUALITY_CONTROL_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <type_traits>
#include <cstddef>
 
#include "tatami/tatami.hpp"
 
#include "find_median_mad.hpp"
#include "per_cell_qc_metrics.hpp"
#include "choose_filter_thresholds.hpp"
#include "utils.hpp"
 
namespace scran_qc {
 
struct ComputeRnaQcMetricsOptions {
    int num_threads = 1;
};
 
template<typename Sum_ = double, typename Detected_ = int, typename Proportion_ = double>
struct ComputeRnaQcMetricsBuffers {
    Sum_* sum = NULL;
 
    Detected_* detected = NULL;
 
    std::vector<Proportion_*> subset_proportion;
};
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_, typename Proportion_>
void compute_rna_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& output,
    const ComputeRnaQcMetricsOptions& options)
{
    const auto NC = mat.ncol();
    const auto nsubsets = subsets.size();
 
    PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_> tmp;
    tmp.sum = output.sum;
    tmp.detected = output.detected;
 
    constexpr bool same_type = std::is_same<Sum_, Proportion_>::value;
    typename std::conditional<same_type, bool, std::vector<std::vector<Sum_> > >::type placeholder_subset;
 
    if (output.subset_proportion.size()) {
        // Providing space for the subset sums if they're not the same type.
        if constexpr(same_type) {
            tmp.subset_sum = output.subset_proportion;
        } else {
            sanisizer::resize(placeholder_subset, nsubsets);
            sanisizer::resize(tmp.subset_sum, nsubsets);
            for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                auto& b = placeholder_subset[s];
                tatami::resize_container_to_Index_size(b, NC);
                tmp.subset_sum[s] = b.data();
            }
        }
    }
 
    PerCellQcMetricsOptions opt;
    opt.num_threads = options.num_threads;
    per_cell_qc_metrics(mat, subsets, tmp, opt);
 
    for (I<decltype(nsubsets)> s = 0 ; s < nsubsets; ++s) {
        const auto dest = output.subset_proportion[s];
        if (dest) {
            const auto src = tmp.subset_sum[s];
            for (Index_ c = 0; c < NC; ++c) {
                dest[c] = static_cast<Proportion_>(src[c]) / static_cast<Proportion_>(tmp.sum[c]);
            }
        }
    }
}
 
template<typename Sum_ = double, typename Detected_ = int, typename Proportion_ = double>
struct ComputeRnaQcMetricsResults {
    std::vector<Sum_> sum;
 
    std::vector<Detected_> detected;
 
    std::vector<std::vector<Proportion_> > subset_proportion;
};
 
template<typename Sum_ = double, typename Detected_ = int, typename Proportion_ = double, typename Value_, typename Index_, typename Subset_>
ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_> compute_rna_qc_metrics(const tatami::Matrix<Value_, Index_>& mat, const std::vector<Subset_>& subsets, const ComputeRnaQcMetricsOptions& options) {
    const auto NC = mat.ncol();
    ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_> buffers;
    ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_> output;
 
    tatami::resize_container_to_Index_size(output.sum, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    buffers.sum = output.sum.data();
 
    tatami::resize_container_to_Index_size(output.detected, NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    buffers.detected = output.detected.data();
 
    const auto nsubsets = subsets.size();
    sanisizer::resize(buffers.subset_proportion, nsubsets);
    sanisizer::resize(output.subset_proportion, nsubsets);
    for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
        tatami::resize_container_to_Index_size(output.subset_proportion[s], NC
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffers.subset_proportion[s] = output.subset_proportion[s].data();
    }
 
    compute_rna_qc_metrics(mat, subsets, buffers, options);
    return output;
}
 
struct ComputeRnaQcFiltersOptions {
    double detected_num_mads = 3;
 
    double sum_num_mads = 3;
 
    double subset_proportion_num_mads = 3;
};
 
namespace internal {
 
template<typename Float_, class Host_, typename Sum_, typename Detected_, typename Proportion_, typename BlockSource_>
void rna_populate(Host_& host, const std::size_t n, const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& res, BlockSource_ block, const ComputeRnaQcFiltersOptions& options) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    auto buffer = [&]{
        if constexpr(unblocked) {
            return sanisizer::create<std::vector<Float_> >(n);
        } else {
            return FindMedianMadWorkspace<Float_>(n, block);
        }
    }();
 
    {
        ChooseFilterThresholdsOptions opts;
        opts.num_mads = options.sum_num_mads;
        opts.log = true;
        opts.upper = false;
        host.get_sum() = [&]{
            if constexpr(unblocked) {
                return choose_filter_thresholds(n, res.sum, buffer.data(), opts).lower;
            } else {
                return internal::strip_threshold<true>(choose_filter_thresholds_blocked(n, res.sum, block, &buffer, opts));
            }
        }();
    }
 
    {
        ChooseFilterThresholdsOptions opts;
        opts.num_mads = options.detected_num_mads;
        opts.log = true;
        opts.upper = false;
        host.get_detected() = [&]{
            if constexpr(unblocked) {
                return choose_filter_thresholds(n, res.detected, buffer.data(), opts).lower;
            } else {
                return internal::strip_threshold<true>(choose_filter_thresholds_blocked(n, res.detected, block, &buffer, opts));
            }
        }();
    }
 
    {
        ChooseFilterThresholdsOptions opts;
        opts.num_mads = options.subset_proportion_num_mads;
        opts.lower = false;
 
        const auto nsubsets = res.subset_proportion.size();
        auto& subhost = host.get_subset_proportion();
        sanisizer::resize(subhost, nsubsets);
        for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
            const auto sub = res.subset_proportion[s];
            subhost[s] = [&]{
                if constexpr(unblocked) {
                    return choose_filter_thresholds(n, sub, buffer.data(), opts).upper;
                } else {
                    return internal::strip_threshold<false>(choose_filter_thresholds_blocked(n, sub, block, &buffer, opts));
                }
            }();
        }
    }
}
 
template<class Host_, typename Sum_, typename Detected_, typename Proportion_, typename BlockSource_, typename Output_>
void rna_filter(const Host_& host, const std::size_t n, const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& metrics, BlockSource_ block, Output_* const output) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    std::fill_n(output, n, 1);
 
    const auto& sum = host.get_sum();
    for (I<decltype(n)> i = 0; i < n; ++i) {
        auto thresh = [&]{
            if constexpr(unblocked) {
                return sum;
            } else {
                return sum[block[i]];
            }
        }();
        output[i] = output[i] && (metrics.sum[i] >= thresh);
    }
 
    const auto& detected = host.get_detected();
    for (I<decltype(n)> i = 0; i < n; ++i) {
        auto thresh = [&]{
            if constexpr(unblocked) {
                return detected;
            } else {
                return detected[block[i]];
            }
        }();
        output[i] = output[i] && (metrics.detected[i] >= thresh);
    }
 
    const auto nsubsets = metrics.subset_proportion.size();
    for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
        const auto sub = metrics.subset_proportion[s];
        const auto& sthresh = host.get_subset_proportion()[s];
        for (I<decltype(n)> i = 0; i < n; ++i) {
            const auto thresh = [&]{
                if constexpr(unblocked) {
                    return sthresh;
                } else {
                    return sthresh[block[i]];
                }
            }();
            output[i] = output[i] && (sub[i] <= thresh);
        }
    }
}
 
template<typename Sum_, typename Detected_, typename Proportion_>
ComputeRnaQcMetricsBuffers<const Sum_, const Detected_, const Proportion_> to_buffer(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics) {
    ComputeRnaQcMetricsBuffers<const Sum_, const Detected_, const Proportion_> buffer;
    buffer.sum = metrics.sum.data();
    buffer.detected = metrics.detected.data();
    buffer.subset_proportion.reserve(metrics.subset_proportion.size());
    for (const auto& s : metrics.subset_proportion) {
        buffer.subset_proportion.push_back(s.data());
    }
    return buffer;
}
 
}
template<typename Float_ = double>
class RnaQcFilters {
public:
    Float_ get_sum() const {
        return my_sum;
    }
 
    Float_ get_detected() const {
        return my_detected;
    }
 
    const std::vector<Float_>& get_subset_proportion() const {
        return my_subset_proportion;
    }
 
    Float_& get_sum() {
        return my_sum;
    }
 
    Float_& get_detected() {
        return my_detected;
    }
 
    std::vector<Float_>& get_subset_proportion() {
        return my_subset_proportion;
    }
 
private:
    Float_ my_sum = 0;
    Float_ my_detected = 0;
    std::vector<Float_> my_subset_proportion;
 
public:
    template<typename Sum_, typename Detected_, typename Proportion_, typename Output_>
    void filter(const std::size_t num, const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& metrics, Output_* const output) const {
        internal::rna_filter(*this, num, metrics, false, output);
    }
 
    template<typename Sum_, typename Detected_, typename Proportion_, typename Output_>
    void filter(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics, Output_* const output) const {
        return filter(metrics.sum.size(), internal::to_buffer(metrics), output);
    }
 
    template<typename Output_ = unsigned char, typename Sum_, typename Detected_, typename Proportion_>
    std::vector<Output_> filter(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics) const {
        auto output = sanisizer::create<std::vector<Output_> >(metrics.sum.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 Proportion_>
RnaQcFilters<Float_> compute_rna_qc_filters(const std::size_t num, const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& metrics, const ComputeRnaQcFiltersOptions& options) {
    RnaQcFilters<Float_> output;
    internal::rna_populate<Float_>(output, num, metrics, false, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Proportion_>
RnaQcFilters<Float_> compute_rna_qc_filters(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics, const ComputeRnaQcFiltersOptions& options) {
    return compute_rna_qc_filters(metrics.sum.size(), internal::to_buffer(metrics), options);
}
 
template<typename Float_ = double>
class RnaQcBlockedFilters {
public:
    const std::vector<Float_>& get_sum() const {
        return my_sum;
    }
 
    const std::vector<Float_>& get_detected() const {
        return my_detected;
    }
 
    const std::vector<std::vector<Float_> >& get_subset_proportion() const {
        return my_subset_proportion;
    }
 
    std::vector<Float_>& get_sum() {
        return my_sum;
    }
 
    std::vector<Float_>& get_detected() {
        return my_detected;
    }
 
    std::vector<std::vector<Float_> >& get_subset_proportion() {
        return my_subset_proportion;
    }
 
private:
    std::vector<Float_> my_sum;
    std::vector<Float_> my_detected;
    std::vector<std::vector<Float_> > my_subset_proportion;
 
public:
    template<typename Sum_, typename Detected_, typename Proportion_, typename Block_, typename Output_>
    void filter(const std::size_t num, const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& metrics, const Block_* const block, Output_* const output) const {
        internal::rna_filter(*this, num, metrics, block, output);
    }
 
    template<typename Sum_, typename Detected_, typename Proportion_, typename Block_, typename Output_>
    void filter(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics, const Block_* const block, Output_* const output) const {
        return filter(metrics.sum.size(), internal::to_buffer(metrics), block, output);
    }
 
    template<typename Output_ = unsigned char, typename Sum_, typename Detected_, typename Proportion_, typename Block_>
    std::vector<Output_> filter(const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics, const Block_* const block) const {
        auto output = sanisizer::create<std::vector<Output_> >(metrics.sum.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 Proportion_, typename Block_>
RnaQcBlockedFilters<Float_> compute_rna_qc_filters_blocked(
    const std::size_t num,
    const ComputeRnaQcMetricsBuffers<Sum_, Detected_, Proportion_>& metrics,
    const Block_* const block,
    const ComputeRnaQcFiltersOptions& options) 
{
    RnaQcBlockedFilters<Float_> output;
    internal::rna_populate<Float_>(output, num, metrics, block, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Proportion_, typename Block_>
RnaQcBlockedFilters<Float_> compute_rna_qc_filters_blocked(
    const ComputeRnaQcMetricsResults<Sum_, Detected_, Proportion_>& metrics,
    const Block_* const block,
    const ComputeRnaQcFiltersOptions& options)
{
    return compute_rna_qc_filters_blocked(metrics.sum.size(), internal::to_buffer(metrics), block, options);
}
 
}
 
#endif