Normalization and filtering

Normalization and filtering#

Library-size normalization and expression filtering for count data.

Counts are samples x features (rows are samples, columns are genes): finite, non-negative, integer or not. They are converted to float64, which holds every whole number up to 2**53 (about 9.0e15) exactly; a larger integer count is rounded to the nearest float64, a relative change below 1.2e-16, and never wraps around. R, whose integers have 32 bits, holds counts that large as doubles too. A library size, a total of a feature’s counts or a count per million beyond the float64 range (about 1.8e308) raises OverflowError instead of becoming infinite or NaN. The functions reproduce edgeR and DESeq2:

  • normalization factors as edgeR’s normLibSizes (calcNormFactors) computes them: tmm_factors(), the trimmed mean of M values (Robinson and Oshlack 2010), upper_quartile_factors() (Bullard et al. 2010) and rle_factors() (Anders and Huber 2010). They multiply the library sizes into effective library sizes, library_sizes * norm_factors, and have geometric mean one;

  • DESeq2’s size factors, which take the place of the library sizes: median_ratio_size_factors(), the median-of-ratios method of Anders and Huber (2010);

  • counts per million and log counts per million as edgeR’s cpm computes them: cpm() and log_cpm();

  • edgeR’s filterByExpr (Chen, Lun and Smyth 2016): filter_by_expression().

The routines repeat R’s floating-point operations wherever an operation decides a discrete outcome: the rank-based trimming and the reference sample of TMM, and the comparisons of filterByExpr. Sums and means that R accumulates in extended precision (sum, mean, colSums, rowMeans) are accumulated with error compensation and rounded once. Logarithms and exponentials come from the C math library, as in R. With these, the results agree with edgeR and DESeq2 to the last bit on the validation inputs (numba set up to use Intel’s SVML library may change the last bit). Conventions that the documentation of edgeR, DESeq2 or R leaves open were found by running them as black boxes; comments mark them “found by probing”.

References

Robinson MD, Oshlack A (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11, R25. doi:10.1186/gb-2010-11-3-r25

Bullard JH, Purdom E, Hansen KD, Dudoit S (2010). Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics 11, 94. doi:10.1186/1471-2105-11-94

Anders S, Huber W (2010). Differential expression analysis for sequence count data. Genome Biology 11, R106. doi:10.1186/gb-2010-11-10-r106

Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5, 1438. doi:10.12688/f1000research.8987.2

cytepool.normalization.cpm(counts, *, library_sizes=None, norm_factors=None)#

Return counts per million, as edgeR’s cpm computes them.

The value is 1e6 * y / (N f) for library size N and normalization factor f.

Parameters:

  • counts – Counts, samples x features.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • norm_factors – Normalization factor of each sample, multiplied into the library sizes; defaults to ones. Must be positive.

Returns:

Float64 array of shape (n_samples, n_features).

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If an argument is out of range, or a sample has library size zero.

  • OverflowError – If a library size, or a count per million as edgeR forms it, exceeds the float64 range.

Explained in: Count models

cytepool.normalization.filter_by_expression(counts, *, group=None, design=None, library_sizes=None, norm_factors=None, min_count=10, min_total_count=15, large_n=10, min_prop=0.7)#

Return which features to keep, by edgeR’s filterByExpr rule.

A feature is kept when its CPM reaches min_count / median(N) * 1e6 in at least n samples and its total count over all samples is at least min_total_count, where N are the effective library sizes (Chen, Lun and Smyth 2016). n is the smallest group size or, with a design, the smallest inverse leverage 1 / max(h), where h is the diagonal of the hat matrix of the design with an intercept column added, as R’s hat() adds it (found by probing); the two agree for a one-way design. When n exceeds large_n it is relaxed to large_n + (n - large_n) * min_prop. Without group and design all samples form one group, as in edgeR, with a warning.

Parameters:

  • counts – Counts, samples x features.

  • group – Group label of each sample; labels must not be missing.

  • design – Design matrix, samples x coefficients, or any object with a matrix attribute holding one, such as a Design. With the intercept column added, its rank is decided as R’s qr() decides it. Only the fixed effects set n: the random effects of a Design are ignored, as edgeR defines n for the model matrix of the fixed effects. Leaving them out can only raise n, and so keep fewer features. When a random effect separates groups whose expression may differ, such as the tissues of a variance partition, pass those groups as group instead.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • norm_factors – Normalization factor of each sample, multiplied into the library sizes; defaults to ones. Must be positive.

  • min_count – Count needed, at the median library size, in at least n samples. At least 0.

  • min_total_count – Total count needed over all samples. At least 0.

  • large_n – Number of samples per group considered large. At least 0.

  • min_prop – Fraction of the samples beyond large_n in which a feature must reach the cutoff, in [0, 1].

Returns:

Boolean array of shape (n_features,), True for features to keep.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If an argument is out of range, both group and design are given (edgeR silently ignores the design then), or a sample has library size zero.

  • OverflowError – If a library size, or the total count of a feature, exceeds the float64 range.

Used in: Quick start, Bulk RNA-seq

Explained in: Choosing a method, Count models

cytepool.normalization.library_sizes(counts)#

Return the library size of each sample: its total count over all features.

The sums are accumulated with error compensation and rounded once, like R’s colSums, so they are exact for whole-number counts whose total is at most 2**53.

Parameters:

counts – Counts, samples x features.

Returns:

Float64 array of shape (n_samples,).

Raises:

  • TypeError – If counts are not numeric, or are a numpy masked array.

  • ValueError – If counts are not a non-empty 2-D array of finite, non-negative values.

  • OverflowError – If a library size exceeds the float64 range.

Used in: Quick start

cytepool.normalization.log_cpm(counts, *, library_sizes=None, norm_factors=None, prior_count=2.0)#

Return log2 counts per million, as edgeR’s cpm(log=TRUE) computes them.

The prior count is scaled by library size, c_j = prior_count * N_j / mean(N) with N the effective library sizes. It is added to every count of sample j and twice to its library size (edgeR’s addPriorCount), giving log2((y + c_j) / (N_j + 2 c_j) * 1e6).

Parameters:

  • counts – Counts, samples x features.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • norm_factors – Normalization factor of each sample, multiplied into the library sizes; defaults to ones. Must be positive.

  • prior_count – Average count added to each observation, at least 0. With 0, zero counts give minus infinity.

Returns:

Float64 array of shape (n_samples, n_features).

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If an argument is out of range, or a sample has library size zero.

  • OverflowError – If a library size, or the sum of the library sizes, exceeds the float64 range.

Explained in: Count models

cytepool.normalization.median_ratio_size_factors(counts, *, method='ratio')#

Return DESeq2’s median-of-ratios size factors.

With method="ratio" (Anders and Huber 2010, equation 5), the size factor of a sample is the median, over the features with positive counts in every sample, of its counts divided by the features’ geometric means across samples. As in DESeq2, the median is taken of the log ratios (found by probing), and the factors are not rescaled.

With method="poscounts", a feature’s geometric mean is the n-th root of the product of its positive counts, n being the number of samples, so features with zeros take part. A sample’s median runs over the features where its count is positive. The factors are then scaled to geometric mean one, as DESeq2’s estimateSizeFactors(), and so DESeq(), returns them (found by probing). The lower-level estimateSizeFactorsForMatrix() returns them before this scaling.

Parameters:

  • counts – Counts, samples x features.

  • method – "ratio" or "poscounts".

Returns:

Float64 array of shape (n_samples,).

Raises:

  • TypeError – If counts are not numeric, or are a numpy masked array.

  • ValueError – If method is unknown; with "ratio", if no feature has positive counts in every sample; with "poscounts", if a sample has no positive count.

Explained in: Count models

cytepool.normalization.rle_factors(counts, *, library_sizes=None)#

Return relative log expression (RLE) factors, as edgeR computes them.

This is the median-ratio method of Anders and Huber (2010) in edgeR’s form. Over the features with positive counts in every sample, a sample’s counts are divided by the features’ geometric means across samples; the factor is the median of these ratios divided by the library size. edgeR takes the median of the ratios, where DESeq2 takes the median of their logarithms (found by probing); the two differ when the number of features is even. The factors are scaled to geometric mean one.

Parameters:

  • counts – Counts, samples x features.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

Returns:

Float64 array of shape (n_samples,) with geometric mean one.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If a sample has library size zero, or no feature has positive counts in every sample (edgeR then returns NaN factors).

  • OverflowError – If a library size exceeds the float64 range.

Explained in: Count models

cytepool.normalization.tmm_factors(counts, *, library_sizes=None, reference=None, log_ratio_trim=0.3, sum_trim=0.05, weighted=True, a_cutoff=-1e10)#

Return TMM normalization factors, as edgeR’s calcNormFactors computes them.

The trimmed mean of M values (Robinson and Oshlack 2010) compares each sample k with a reference sample r over the features with positive counts in both. With proportions p = y / N, where N is the library size, a feature has log ratio M = log2(p_k / p_r) and average log proportion A = (log2 p_k + log2 p_r) / 2. Features with A <= a_cutoff are dropped. Of the n features left, those whose M lies in the lowest or highest log_ratio_trim fraction, or whose A lies in the lowest or highest sum_trim fraction, are trimmed: a feature is kept when its average rank lies within floor(n * trim) + 1 and n - floor(n * trim) for both M and A, so tied values are kept or trimmed together. The factor of sample k is 2 to the mean of the remaining M values, weighted by 1 / v with the delta-method variance v = (N_k - y_k) / (N_k y_k) + (N_r - y_r) / (N_r y_r). A sample with no feature left, or whose log ratios are all smaller than 1e-6 in absolute value, gets factor 1. Finally the factors are scaled to geometric mean one.

Features with zero counts in every sample are left out, as edgeR removes them; this matters only for the upper quartiles that pick the default reference.

Parameters:

  • counts – Counts, samples x features.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • reference – Index of the reference sample. The default is edgeR’s: the sample whose upper quartile of counts divided by its library size is closest to the mean of that ratio over samples (the first one when several are equally close). When the median of the ratios is below 1e-20, as with very sparse counts, it is the sample with the largest sum of square-root counts instead (found by probing).

  • log_ratio_trim – Fraction of features trimmed from each end of the M values, in [0, 0.5).

  • sum_trim – Fraction of features trimmed from each end of the A values, in [0, 0.5).

  • weighted – Whether to weight the M values by their inverse variances.

  • a_cutoff – Features with A at or below this value are not used.

Returns:

Float64 array of shape (n_samples,) with geometric mean one.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If an argument is out of range, or a sample has library size zero.

  • OverflowError – If a library size exceeds the float64 range.

Used in: Quick start, Bulk RNA-seq

Explained in: Count models

cytepool.normalization.upper_quartile_factors(counts, *, library_sizes=None, p=0.75)#

Return upper-quartile normalization factors, as edgeR computes them.

The factor of a sample is the p quantile of its counts (R’s default quantile, type 7) over the features that are not zero in every sample, divided by its library size (Bullard et al. 2010). The factors are scaled to geometric mean one.

Parameters:

  • counts – Counts, samples x features.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • p – Quantile to use, in [0, 1].

Returns:

Float64 array of shape (n_samples,) with geometric mean one.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – If an argument is out of range, a sample has library size zero, or the quantile is zero in some sample (edgeR then returns NaN factors).

  • OverflowError – If a library size exceeds the float64 range.

Explained in: Count models