SuggestRnaQcFilters.hpp

Go to the documentation of this file.

#ifndef SCRAN_SUGGEST_RNA_QC_FILTERS_H
#define SCRAN_SUGGEST_RNA_QC_FILTERS_H
 
#include "../utils/macros.hpp"
 
#include <vector>
#include <cstdint>
 
#include "PerCellRnaQcMetrics.hpp"
#include "ComputeMedianMad.hpp"
#include "ChooseOutlierFilters.hpp"
 
namespace scran {
 
class SuggestRnaQcFilters {
public:
    struct Defaults {
        static constexpr double num_mads = 3;
    };
 
private:
    double detected_num_mads = Defaults::num_mads;
    double sums_num_mads = Defaults::num_mads;
    double subset_num_mads = Defaults::num_mads;
 
public:
    SuggestRnaQcFilters& set_detected_num_mads(double n = Defaults::num_mads) {
        detected_num_mads = n;
        return *this;
    }
 
    SuggestRnaQcFilters& set_sums_num_mads(double n = Defaults::num_mads) {
        sums_num_mads = n;
        return *this;
    }
 
    SuggestRnaQcFilters& set_subset_num_mads(double n = Defaults::num_mads) {
        subset_num_mads = n;
        return *this;
    }
 
    SuggestRnaQcFilters& set_num_mads(double n = Defaults::num_mads) {
        detected_num_mads= n;
        subset_num_mads = n;
        sums_num_mads = n;
        return *this;
    }
 
public:
    struct Thresholds {
        std::vector<double> sums;
 
        std::vector<double> detected;
 
        std::vector<std::vector<double> > subset_proportions;
 
    public:
        template<bool overwrite = true, typename Float, typename Integer, typename Output>
        void filter(size_t n, const PerCellRnaQcMetrics::Buffers<Float, Integer>& buffers, Output* output) const {
            if (detected.size() != 1) {
                throw std::runtime_error("should use filter_blocked() for multiple batches");
            }
            filter_<overwrite>(n, buffers, output, [](size_t i) -> size_t { return 0; });
        }
 
        template<typename Output = uint8_t>
        std::vector<Output> filter(const PerCellRnaQcMetrics::Results& metrics) const {
            std::vector<Output> output(metrics.detected.size());
            filter(output.size(), metrics.buffers(), output.data());
            return output;
        }
 
    public:
        template<bool overwrite = true, typename Block, typename Float, typename Integer, typename Output>
        void filter_blocked(size_t n, const Block* block, const PerCellRnaQcMetrics::Buffers<Float, Integer>& buffers, Output* output) const {
            if (block) {
                filter_<overwrite>(n, buffers, output, [&](size_t i) -> Block { return block[i]; });
            } else {
                filter<overwrite>(n, buffers, output);
            }
        }
 
        template<typename Output = uint8_t, typename Block>
        std::vector<Output> filter_blocked(const PerCellRnaQcMetrics::Results& metrics, const Block* block) const {
            std::vector<Output> output(metrics.detected.size());
            filter_blocked(output.size(), block, metrics.buffers(), output.data());
            return output;
        }
 
    private:
        template<bool overwrite, typename Float, typename Integer, typename Output, typename Function>
        void filter_(size_t n, const PerCellRnaQcMetrics::Buffers<Float, Integer>& buffers, Output* output, Function indexer) const {
            size_t nsubsets = subset_proportions.size();
 
            for (size_t i = 0; i < n; ++i) {
                auto b = indexer(i);
                if (quality_control::is_less_than(buffers.sums[i], sums[b])) {
                    output[i] = true;
                    continue;
                }
 
                if (quality_control::is_less_than<double>(buffers.detected[i], detected[b])) {
                    output[i] = true;
                    continue;
                }
 
                bool fail = false;
                for (size_t s = 0; s < nsubsets; ++s) {
                    if (quality_control::is_greater_than(buffers.subset_proportions[s][i], subset_proportions[s][b])) {
                        fail = true;
                        break;
                    }
                }
                if (fail) {
                    output[i] = true;
                    continue;
                }
 
                if constexpr(overwrite) {
                    output[i] = false;
                }
            }
 
            return;
        }
    };
 
public:
    template<typename Float, typename Integer>
    Thresholds run(size_t n, const PerCellRnaQcMetrics::Buffers<Float, Integer>& buffers) const {
        return run_blocked(n, static_cast<int*>(NULL), buffers);
    }
 
    Thresholds run(const PerCellRnaQcMetrics::Results& metrics) const {
        return run(metrics.detected.size(), metrics.buffers());
    }
 
public:
    template<typename Block, typename Float, typename Integer>
    Thresholds run_blocked(size_t n, const Block* block, const PerCellRnaQcMetrics::Buffers<Float, Integer>& buffers) const {
        Thresholds output;
        std::vector<double> workspace(n);
        std::vector<int> starts;
        if (block) {
            starts = ComputeMedianMad::compute_block_starts(n, block);
        }
 
        // Filtering on the total counts.
        {
            ComputeMedianMad meddler;
            meddler.set_log(true);
            auto sums_res = meddler.run_blocked(n, block, starts, buffers.sums, workspace.data());
 
            ChooseOutlierFilters filter;
            filter.set_num_mads(sums_num_mads);
            filter.set_upper(false);
            auto sums_filt = filter.run(std::move(sums_res));
 
            output.sums = std::move(sums_filt.lower);
        }
 
        // Filtering on the detected features.
        {
            ComputeMedianMad meddler;
            meddler.set_log(true);
            auto detected_res = meddler.run_blocked(n, block, starts, buffers.detected, workspace.data());
 
            ChooseOutlierFilters filter;
            filter.set_num_mads(sums_num_mads);
            filter.set_upper(false);
            auto detected_filt = filter.run(std::move(detected_res));
 
            output.detected = std::move(detected_filt.lower);
        }
 
        // Filtering on the subset proportions (no log). 
        {
            ComputeMedianMad meddler;
 
            ChooseOutlierFilters filter;
            filter.set_num_mads(subset_num_mads);
            filter.set_lower(false);
 
            size_t nsubsets = buffers.subset_proportions.size();
            for (size_t s = 0; s < nsubsets; ++s) {
                auto subset_res = meddler.run_blocked(n, block, starts, buffers.subset_proportions[s], workspace.data());
                auto subset_filt = filter.run(std::move(subset_res));
                output.subset_proportions.push_back(std::move(subset_filt.upper));
            }
        }
 
        return output;
    }
 
    template<typename Block>
    Thresholds run_blocked(const PerCellRnaQcMetrics::Results& metrics, const Block* block) const {
        return run_blocked(metrics.detected.size(), block, metrics.buffers());
    }
};
 
}
 
#endif