Why Compare Libraries Gene by Gene?
RNA-seq library sizes can differ, but dividing every count by its library size does not always make samples comparable. If a small number of genes occupy a large share of one library, the remaining genes can appear to have lower proportions even when their underlying expression has not changed. TMM, or the trimmed mean of M-values, estimates a relative scaling adjustment by comparing gene-wise proportions between a sample and a reference sample.
The aim is not to declare that most genes are unchanged with certainty. Rather, TMM uses the working assumption that most genes, or a sufficiently balanced set of genes, provide a stable basis for estimating relative library composition. Trimming reduces the influence of genes with extreme relative changes and genes with very low or very high average abundance.
An \(M\)-value of \(1\) means the target gene’s proportion is twice the reference gene’s proportion; an \(M\)-value of \(-1\) means it is half. The \(A\)-value describes the abundance scale, not the direction of the difference. TMM trims on both statistics: genes at the extremes of the \(M\)-value distribution and genes at the extremes of the \(A\)-value distribution are excluded. The final contributing set is the intersection of the two retained sets.
Assumptions and the TMM Estimate
TMM requires count data and library sizes defined consistently with those counts. For a target-reference comparison, genes with a zero in either library do not have a finite log-ratio and are omitted from that comparison. Its central biological and statistical assumption is that the retained genes give a representative estimate of relative composition: broadly, most genes are not changing in one direction, or opposing changes are sufficiently balanced. If nearly all genes shift together, the data alone generally cannot distinguish that biological shift from a change in overall composition.
After trimming, TMM takes a precision-weighted mean of the retained \(M\)-values. A commonly used approximation for the sampling variance of \(M_g^{(j,r)}\), up to a constant common to all genes, is $$ v_g^{(j,r)}=\frac{N_j-Y_{gj}}{N_jY_{gj}}+\frac{N_r-Y_{gr}}{N_rY_{gr}}. $$ The corresponding weight is \(w_g^{(j,r)}=1/v_g^{(j,r)}\). The target’s uncentered TMM factor is $$ f_j^{(r)}=2^{\widehat{m}_{j,r}}, \qquad \widehat{m}_{j,r}= \frac{\sum_{g\in R_{j,r}}w_g^{(j,r)}M_g^{(j,r)}} {\sum_{g\in R_{j,r}}w_g^{(j,r)}}, $$ where \(R_{j,r}\) is the set of genes retained after both trimming steps. A positive weighted mean gives a factor above 1; a negative mean gives a factor below 1.
Why the Weights Approximate Precision
The inverse-variance weights follow from a simple sampling argument. If a gene has count \(Y\) out of a library of size \(N\), its estimated proportion is \(Y/N\). Under a binomial approximation, the variance of that proportion is approximately \(p(1-p)/N\). Applying the delta method to its logarithm gives an approximate variance of \((1-p)/(Np)\), which, on substituting \(p\approx Y/N\), is \(1/Y-1/N\). For two independently sampled libraries, the variances of their log proportions add. Changing from natural logs to base-2 logs multiplies the variance by the same constant \(1/(\ln 2)^2\) for every gene, so that constant cancels when forming normalized weights.
This is an approximation: RNA-seq counts can be more variable than a binomial sampling model allows, and the formula does not model biological variation between replicates. It is used here to weight the gene-wise ratios for the normalization calculation, not as a complete model for differential expression.
Worked Example: Identify the Genes Removed by Trimming
Suppose a hypothetical target and reference each have library size \(N=10{,}000\). The following eight genes are selected from those libraries. Use the stated counts to calculate their \(M\)- and \(A\)-values, then apply a simplified trim of 12.5% from each tail of each statistic. For ties at the upper \(M\)-tail, remove gene 1.
| Gene | Target count | Reference count | \(M\) | \(A\), approximately |
|---|---|---|---|---|
| 1 | 20 | 10 | 1 | -9.466 |
| 2 | 40 | 20 | 1 | -8.466 |
| 3 | 80 | 80 | 0 | -6.966 |
| 4 | 160 | 160 | 0 | -5.966 |
| 5 | 320 | 320 | 0 | -4.966 |
| 6 | 640 | 640 | 0 | -3.966 |
| 7 | 1,280 | 1,280 | 0 | -2.966 |
| 8 | 4,000 | 8,000 | -1 | -0.822 |
For gene 1, the proportions are \(20/10{,}000=0.002\) and \(10/10{,}000=0.001\), so \(M_1=\log_2(0.002/0.001)=1\) and \(A_1=\tfrac12\log_2(0.002\times0.001)\approx-9.466\). For gene 8, the proportions are \(0.4\) and \(0.8\), so \(M_8=\log_2(0.4/0.8)=-1\) and \(A_8=\tfrac12\log_2(0.4\times0.8)\approx-0.822\). The other table entries follow in the same way.
With eight genes, 12.5% of the genes is one gene. The \(M\)-trim removes gene 8 at the low end and, by the stated tie rule, gene 1 at the high end. The \(A\)-trim removes gene 1 at the low-abundance end and gene 8 at the high-abundance end. Taking the intersection leaves genes 2 through 7. This illustrates that the two trims act on different features: gene 2 is extreme in \(M\) but remains because only one upper-tail gene is trimmed, while genes 1 and 8 are removed by the abundance trim as well.
Calculating a Weighted TMM Factor
In practice, trimming proportions are chosen by the implementation and analysis settings. For example, commonly used edgeR defaults trim 30% from each tail of the \(M\)-values and 5% from each tail of the \(A\)-values. These are not mathematical requirements; a small teaching example can use a pre-specified retained set to focus on the weighted calculation.
Worked Example: Compute a TMM Factor From Retained Genes
Consider a hypothetical target-reference comparison with \(N_j=N_r=10{,}000\). Suppose trimming has retained three genes with counts \((Y_{gj},Y_{gr})\) equal to \((100,50)\), \((200,100)\), and \((300,300)\). Assume all three have positive counts and are retained; other genes in the full libraries are not shown.
The first two genes have \(M=\log_2(2)=1\), and the third has \(M=\log_2(1)=0\). For the first gene, the weight is $$ w_1=\left(\frac{10{,}000-100}{10{,}000(100)} +\frac{10{,}000-50}{10{,}000(50)}\right)^{-1} =(0.0099+0.0199)^{-1}\approx33.56. $$ For the second gene, $$ w_2=(0.0049+0.0099)^{-1}\approx67.57, $$ and for the third, $$ w_3=(0.003233+0.003233)^{-1}\approx154.6. $$ The third gene receives the largest weight because its two counts are larger and its estimated ratio is more precise under the sampling approximation.
The weighted mean is $$ \widehat{m}_{j,r} =\frac{33.56(1)+67.57(1)+154.6(0)} {33.56+67.57+154.6} \approx0.3954. $$ Therefore the uncentered factor is $$ f_j^{(r)}=2^{0.3954}\approx1.315. $$ The target’s relative composition factor is about 1.315 before across-sample centering. This calculation does not mean that every target count is multiplied by 1.315; it is used to adjust the target’s library-size offset.
Centering Factors Across Samples
TMM comparisons are relative, so the numerical factors depend on how the factors are centered. For \(J\) samples, if the initial factors are \(f_1,\ldots,f_J\), a common convention divides each by their geometric mean: $$ \widetilde f_j=\frac{f_j}{\left(\prod_{k=1}^{J}f_k\right)^{1/J}}. $$ Then \(\prod_{j=1}^{J}\widetilde f_j=1\). The rescaling preserves relative factors while avoiding an arbitrary overall shift in all effective library sizes.
Worked Example: Center Factors and Find Effective Library Sizes
Suppose three hypothetical samples have initial factors \(1.20\), \(0.80\), and \(1.05\). Their geometric mean is $$ (1.20\times0.80\times1.05)^{1/3} =(1.008)^{1/3}\approx1.00266. $$ The centered factors are therefore approximately \(1.197\), \(0.7979\), and \(1.047\), respectively. Their product is approximately 1, allowing for rounding.
If their library sizes are 20 million, 30 million, and 25 million reads, the effective library sizes are approximately \(20\times1.197=23.94\) million, \(30\times0.7979=23.94\) million, and \(25\times1.047=26.18\) million reads. These effective sizes incorporate the estimated relative composition adjustments. The original observed counts remain counts; the normalization is applied through the effective library size in downstream modeling.
Diagnostics, Alternatives, and Limitations
A TMM estimate can be unreliable if its assumptions are poorly supported or if the retained genes provide little information. Useful checks include plotting \(M\) against \(A\) for each sample-reference comparison, confirming that the retained ratios are not dominated by a narrow group of genes, and checking that factors are not driven by a few observations. Many extreme ratios at low abundance are a reason to check zero handling and filtering choices; a large shift across most genes in one direction may indicate a composition change that trimming cannot resolve.
Median-ratio size-factor estimation, covered in the earlier material on median-ratio size-factor estimation, is a different robust summary of between-sample gene-wise comparisons. TMM uses log-ratios and average abundance, trims both distributions, and weights retained ratios by approximate precision. Neither approach creates information about absolute RNA output when only relative sequencing counts are available. The appropriate normalization method depends on the data and assumptions; the factor calculation alone cannot verify those assumptions.
A frequent mistake is to trim only genes with large \(M\)-values and to describe the result as TMM. That ignores the separate abundance trim, which can remove genes with extreme \(A\)-values even when their \(M\)-values are moderate. Another mistake is to treat a factor above 1 as evidence that the target has more total RNA per cell. The factor describes a relative adjustment between libraries under the normalization assumptions, not a per-cell measurement.
TMM also has practical limits. It cannot correct sample swaps, poor-quality libraries, or an unsuitable experimental design. It may be sensitive when the number of informative shared genes is small, when counts are mostly zero, or when many genes change in a coordinated direction. Inspecting the data and reporting the normalization method and key settings are part of a defensible analysis.
Check Your Understanding
Try each problem, then compare your reasoning with the solutions below.
- If the target and reference proportions for a gene are 0.006 and 0.003, respectively, calculate its \(M\)-value and interpret it.
- What determines whether a gene is removed by the \(A\)-trim, and how does that differ from the \(M\)-trim?
- Suppose three retained genes have \(M\)-values 1, 0, and \(-1\), with weights 2, 5, and 1. Calculate the weighted mean and the uncentered TMM factor.
- Why are factors usually centered to have geometric mean 1 across samples? Does this make the factors estimates of absolute RNA per cell?
Solutions. (1) \(M=\log_2(0.006/0.003)=\log_2(2)=1\); the target proportion is twice the reference proportion for this gene. (2) The \(A\)-trim removes genes in the extreme tails of average log-abundance, whereas the \(M\)-trim removes genes in the extreme tails of target-to-reference log-ratio. A gene can be removed by either trim, and only genes retained by both contribute. (3) The weighted mean is \([2(1)+5(0)+1(-1)]/(2+5+1)=1/8=0.125\). The factor is \(2^{0.125}\approx1.091\). (4) Centering sets a common relative scale for the factors without changing their ratios; it does not measure absolute RNA per cell, which cannot generally be identified from relative sequencing counts alone.