scran_variances
Model per-gene variance in expression
Loading...
Searching...
No Matches
fit_variance_trend.hpp
Go to the documentation of this file.
1#ifndef SCRAN_VARIANCES_FIT_VARIANCE_TREND_H
2#define SCRAN_VARIANCES_FIT_VARIANCE_TREND_H
3
4#include <algorithm>
5#include <vector>
6#include <array>
7#include <cstddef>
8
10#include "sanisizer/sanisizer.hpp"
11
12#include "utils.hpp"
13
19namespace scran_variances {
20
32 double minimum_mean = 0.1;
33
38 bool mean_filter = true;
39
45 bool transform = true;
46
52 double span = 0.3;
53
62 bool use_minimum_width = false;
63
69 double minimum_width = 1;
70
77
82 int num_threads = 1;
83};
84
90template<typename Float_>
96
97 std::vector<unsigned char> sort_workspace;
98
99 std::vector<Float_> xbuffer, ybuffer;
103};
104
130template<typename Float_>
132 const std::size_t n,
133 const Float_* const mean,
134 const Float_* const variance,
135 Float_* const fitted,
136 Float_* const residuals,
138 const FitVarianceTrendOptions& options)
139{
140 auto& xbuffer = workspace.xbuffer;
141 sanisizer::resize(xbuffer, n);
142 auto& ybuffer = workspace.ybuffer;
143 sanisizer::resize(ybuffer, n);
144
145 const auto quad = [](Float_ x) -> Float_ {
146 return x * x * x * x;
147 };
148
149 std::size_t counter = 0;
150 const Float_ min_mean = options.minimum_mean;
151 for (decltype(I(n)) i = 0; i < n; ++i) {
152 if (!options.mean_filter || mean[i] >= min_mean) {
153 xbuffer[counter] = mean[i];
154 if (options.transform) {
155 ybuffer[counter] = std::pow(variance[i], 0.25); // Using the same quarter-root transform that limma::voom uses.
156 } else {
157 ybuffer[counter] = variance[i];
158 }
159 ++counter;
160 }
161 }
162
163 if (counter < 2) {
164 throw std::runtime_error("not enough observations above the minimum mean");
165 }
166
167 auto& sorter = workspace.sorter;
168 sorter.set(counter, xbuffer.data());
169 auto& work = workspace.sort_workspace;
170 sorter.permute(std::array<Float_*, 2>{ xbuffer.data(), ybuffer.data() }, work);
171
173 if (options.use_minimum_width) {
174 smooth_opt.span = options.minimum_window_count;
175 smooth_opt.span_as_proportion = false;
176 smooth_opt.minimum_width = options.minimum_width;
177 } else {
178 smooth_opt.span = options.span;
179 }
180 smooth_opt.num_threads = options.num_threads;
181
182 // Using the residual array to store the robustness weights as a placeholder;
183 // we'll be overwriting this later.
184 WeightedLowess::compute(counter, xbuffer.data(), ybuffer.data(), fitted, residuals, smooth_opt);
185
186 // Determining the left edge before we unpermute.
187 const Float_ left_x = xbuffer[0];
188 const Float_ left_fitted = (options.transform ? quad(fitted[0]) : fitted[0]);
189
190 sorter.unpermute(fitted, work);
191
192 // Walking backwards to shift the elements back to their original position
193 // (i.e., before filtering on the mean) on the same array. We need to walk
194 // backwards to ensure that writing to the original position on this array
195 // doesn't clobber the first 'counter' positions containing the fitted
196 // values, at least not until each value is shifted to its original place.
197 for (auto i = n; i > 0; --i) {
198 auto j = i - 1;
199 if (!options.mean_filter || mean[j] >= min_mean) {
200 --counter;
201 fitted[j] = (options.transform ? quad(fitted[counter]) : fitted[counter]);
202 } else {
203 fitted[j] = mean[j] / left_x * left_fitted; // draw a y = x line to the origin from the left of the fitted trend.
204 }
205 }
206
207 for (decltype(I(n)) i = 0; i < n; ++i) {
208 residuals[i] = variance[i] - fitted[i];
209 }
210 return;
211}
212
221template<typename Float_>
227
228 FitVarianceTrendResults(const std::size_t n) :
229 fitted(sanisizer::cast<decltype(I(fitted.size()))>(n)
230#ifdef SCRAN_VARIANCES_TEST_INIT
231 , SCRAN_VARIANCES_TEST_INIT
232#endif
233 ),
234 residuals(sanisizer::cast<decltype(I(residuals.size()))>(n)
235#ifdef SCRAN_VARIANCES_TEST_INIT
236 , SCRAN_VARIANCES_TEST_INIT
237#endif
238 )
239 {}
247 std::vector<Float_> fitted;
248
252 std::vector<Float_> residuals;
253};
254
267template<typename Float_>
268FitVarianceTrendResults<Float_> fit_variance_trend(const std::size_t n, const Float_* const mean, const Float_* const variance, const FitVarianceTrendOptions& options) {
271 fit_variance_trend(n, mean, variance, output.fitted.data(), output.residuals.data(), work, options);
272 return output;
273}
274
275}
276
277#endif
void compute(const std::size_t num_points, const Data_ *x, const PrecomputedWindows< Data_ > &windows, const Data_ *y, Data_ *fitted, Data_ *robust_weights, const Options< Data_ > &opt)
Variance modelling for single-cell expression data.
Definition choose_highly_variable_genes.hpp:15
void fit_variance_trend(const std::size_t n, const Float_ *const mean, const Float_ *const variance, Float_ *const fitted, Float_ *const residuals, FitVarianceTrendWorkspace< Float_ > &workspace, const FitVarianceTrendOptions &options)
Definition fit_variance_trend.hpp:131
Options for fit_variance_trend().
Definition fit_variance_trend.hpp:24
int num_threads
Definition fit_variance_trend.hpp:82
double minimum_width
Definition fit_variance_trend.hpp:69
bool use_minimum_width
Definition fit_variance_trend.hpp:62
double span
Definition fit_variance_trend.hpp:52
bool transform
Definition fit_variance_trend.hpp:45
int minimum_window_count
Definition fit_variance_trend.hpp:76
bool mean_filter
Definition fit_variance_trend.hpp:38
double minimum_mean
Definition fit_variance_trend.hpp:32
Results of fit_variance_trend().
Definition fit_variance_trend.hpp:222
std::vector< Float_ > residuals
Definition fit_variance_trend.hpp:252
std::vector< Float_ > fitted
Definition fit_variance_trend.hpp:247
Workspace for fit_variance_trend().
Definition fit_variance_trend.hpp:91