Matthew E. RitchiepersName, Belinda PhipsonpersName, Di WupersName, Yifang HupersName, Charity W. LawpersName, Wei ShipersName, Gordon K. SmythpersName (2015date). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Researchtitle·j.
Object of study — limmaterm's statistical philosophy, expression-analysis workflow, RNA-seqterm extensions, gene-set methods, object model, and interfaces
Research question — how to analyse small, complex gene-expression experiments flexibly and reliably across microarray and RNA-seqterm technologies
Contribution — a review and synthesis of limmaterm's historical and new capabilities from import and preprocessing through linear modelling, empirical Bayes inference, and gene-set analysis
Cross-disciplinary analysis
content-tags & facets:
phenomenon · claim-type · cross-domain-transfer · mechanism · claim-relation · replication-status · limitation — grounded in the full text| phenomenon | differential gene expression, differential splicing, and higher-order expression signatures in small complex experiments — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | |
| claim-type | existence — “This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described. The article outlines limma 's functionality at each of the main steps in a gene expression analysis, from data import, pre-processing, quality assessment and normalization, through to linear modelling, DE and gene signature analyses.” | |
| method — evidence | — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | |
| object of study — evidence | — “This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described. The article outlines limma 's functionality at each of the main steps in a gene expression analysis, from data import, pre-processing, quality assessment and normalization, through to linear modelling, DE and gene signature analyses.” | |
| research question — evidence | — “Despite the complexity, gene expression studies often involve only a small number of biological replicates. The small but complex nature of gene expression studies poses challenging statistical problems and motivates the use of a number of specialized statistical techniques in order to get the most out of each data set. We have developed the limma software over the past decade to provide a framework for analysing gene expression experiments from beginning to end in a flexible and statistically rigorous way.” | |
| contribution — evidence | — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | |
| mechanism | whole-experiment modelling One integrated model shares information and represents sample correlation across the full design. — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | mechanisms[0] |
| mechanism | global parameters Parameters estimated from all genes link otherwise separate gene-wise models. — “A separate model is fitted for each gene, but the gene-wise models can be linked by global parameters or global hyper-parameters. The use of global parameters is a simple means of sharing information between genes that can be used even for the smallest experiments, because the global parameters can be estimated from the entire data set involving all the genes at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same gene, or correlations between related RNA samples, or variations in quality between the RNA samples.” | mechanisms[1] |
| mechanism | empirical Bayes moderation Each residual variance is shrunk toward pooled global variability to increase effective degrees of freedom. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | mechanisms[2] |
| mechanism | robust empirical Bayes Hyper-variable genes receive separate weighting instead of uniform shrinkage. — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | mechanisms[3] |
| mechanism | quantitative weights Less reliable genes, samples, or observations are down-weighted without being discarded. — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | mechanisms[4] |
| mechanism | voom The empirical count mean-variance trend becomes precision weights for logCPM linear modelling. — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | mechanisms[5] |
| mechanism | normalization Platform-specific procedures align samples and remove technical effects unrelated to biology. — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | mechanisms[6] |
| mechanism | multidimensional scaling Leading-fold-change distances reveal sample clusters and unexpected batch patterns. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | mechanisms[7] |
| mechanism | consensus correlation DuplicateCorrelation estimates one intrablock correlation shared across genes and inserts it into all DE tests. — “limma includes a unique strategy for incorporating the fact that observations or samples may be correlated. The strategy is similar to fitting a random effects model, with the difference that all genes are constrained to share the same intrablock correlation. The duplicateCorrelation function is used to estimate the consensus correlation. The correlation structure is then incorporated into the linear model fit and hence into all tests for DE. Originally the idea was used to estimate the correlation between replicate copies of the same probe on a microarray ( 44 ). The correlation strategy preserves more information than simply averaging the replicate probe copies. More generally, the same idea is also used to model the correlation between related RNA samples, for example repeated measures on the same individual or RNA samples collected at the same time.” | mechanisms[8] |
| mechanism | robust variance shrinkage Gene-specific factors protect outlying large variances while shrinking stable genes more strongly. — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | mechanisms[9] |
| mechanism | camera Inter-gene correlation produces a variance-inflation adjustment for competitive gene-set statistics. — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | mechanisms[10] |
| mechanism | rotation Residuals rotate while non-tested model coefficients remain fixed. — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | mechanisms[11] |
| builds-on | Bioconductor — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | claim_relations[0] |
| extends | microarray analysis — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | claim_relations[1] |
| builds-on | gene-wise linear models — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | claim_relations[2] |
| complements | edgeR — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | claim_relations[3] |
| sense | limma as an integrated R/Bioconductor expression-analysis framework — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | senses[0] |
| sense | linear modelling as one whole-experiment model rather than pairwise fragments — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | senses[1] |
| sense | empirical Bayes as compromise variance estimation across gene-wise and global information — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | senses[2] |
| sense | weights as graduated reliability rather than binary filtering — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | senses[3] |
| sense | voom as empirical mean-variance precision modelling for log-counts — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | senses[4] |
| sense | normalization as technical-effect removal rather than biological-signal removal — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | senses[5] |
| sense | leading fold change as RMS separation among the genes best distinguishing samples — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | senses[6] |
| sense | biological correlation as the error-free association between contrast effects — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | senses[7] |
| sense | camera as competitive testing corrected for inter-gene dependence — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | senses[8] |
| sense | rotation as residual-space simulation that preserves non-tested coefficients — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | senses[9] |
| sense | limma objects as list-based, matrix-analogous containers for each workflow stage — “A simple but appropriate object-oriented paradigm provides users with a consistent analysis interface that is very easy from a user point of view. Many limma functions are generic or operate appropriately on objects of different classes. limma defines a number of classes that have been tailored to handle both microarray and RNA-seq data. The philosophy has been to define simple list-based data objects that can be easily explored and manipulated by users, in the same style as familiar, long-standing core functions in R such as lm and glm.” | senses[10] |
| sense | limma as a review subject and existing infrastructure, not a newly introduced method in this article — “This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described. The article outlines limma 's functionality at each of the main steps in a gene expression analysis, from data import, pre-processing, quality assessment and normalization, through to linear modelling, DE and gene signature analyses.” | senses[11] |
| control | control probes with quantitative weights — held — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” — rules out normalization domination by unreliable probes — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | controls[0] |
| control | robust empirical Bayes shrinkage — held — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” — rules out false significance caused by over-shrinking hyper-variable genes — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | controls[1] |
| control | non-tested linear-model coefficients during rotation — held — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” — rules out variation in nuisance coefficients during a gene-set test — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | controls[2] |
| note | Bioconductor — limma is a core component of Bioconductor. — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | claim_relations[0] |
| note | microarray analysis — Capabilities formerly limited to microarrays now apply to RNA-seq. — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | claim_relations[1] |
| weight | microarray analysis — 0.95 — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | claim_relations[1] |
| note | gene-wise linear models — The package fits one linear model per feature and links them by shared information. — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | claim_relations[2] |
| note | edgeR — limma accepts DGEList objects and edgeR normalization factors. — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | claim_relations[3] |
| weight | edgeR — 0.95 — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | claim_relations[3] |
| note | limma —is-a→ R/Bioconductor software package — Explicit categorical definition. — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concept_edges[0] |
| relation | limma —is-a→ R/Bioconductor software package — limma is an R/Bioconductor software package. — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concept_edges[0] |
| note | Bioconductor —is-a→ R-based open-source software development project — Explicit categorical definition. — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | concept_edges[1] |
| relation | Bioconductor —is-a→ R-based open-source software development project — Bioconductor is an R-based open-source software development project. — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | concept_edges[1] |
| note | linear models —contrasted-with→ piece-meal comparisons between pairs of treatments — Explicit rather-than comparison. — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | concept_edges[2] |
| relation | linear models —contrasted-with→ piece-meal comparisons between pairs of treatments — Whole-experiment linear models are used rather than piece-meal pairwise comparisons. — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | concept_edges[2] |
| weight | linear models —contrasted-with→ piece-meal comparisons between pairs of treatments — 0.95 — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | concept_edges[2] |
| effect direction | empirical Bayes procedures —causes→ statistical power and accuracy — increases — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concept_edges[3] |
| note | empirical Bayes procedures —causes→ statistical power and accuracy — Flat positive improvement claim. — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concept_edges[3] |
| relation | empirical Bayes procedures —causes→ statistical power and accuracy — Flexible empirical Bayes enhancements improve statistical power and accuracy. — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concept_edges[3] |
| weight | empirical Bayes procedures —causes→ statistical power and accuracy — 0.95 — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concept_edges[3] |
| effect direction | quantitative weights —causes→ power to detect differentially expressed genes — increases — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concept_edges[4] |
| note | quantitative weights —causes→ power to detect differentially expressed genes — Explicit positive increase. — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concept_edges[4] |
| relation | quantitative weights —causes→ power to detect differentially expressed genes — The use of quantitative weights increases power to detect differentially expressed genes. — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concept_edges[4] |
| weight | quantitative weights —causes→ power to detect differentially expressed genes — 0.95 — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concept_edges[4] |
| note | Removing this non-specific signal —same-as→ background correction — Explicit naming equivalence. — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | concept_edges[5] |
| relation | Removing this non-specific signal —same-as→ background correction — Removing the non-specific probe signal is called background correction. — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | concept_edges[5] |
| note | leading fold change —operationalised-as→ root-mean-square average of the log-fold-changes — Exact operational definition. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concept_edges[6] |
| relation | leading fold change —operationalised-as→ root-mean-square average of the log-fold-changes — Leading fold change is defined as the root-mean-square average of selected log-fold-changes. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concept_edges[6] |
| note | limma —distinct-from→ permutation or re-sampling-based inference — Explicit negative distinction. — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | concept_edges[7] |
| relation | limma —distinct-from→ permutation or re-sampling-based inference — limma does not do permutation or re-sampling-based inference. — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | concept_edges[7] |
| weight | limma —distinct-from→ permutation or re-sampling-based inference — 0.95 — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | concept_edges[7] |
| note | greater than a threshold —contrasted-with→ different to zero — Explicit rather-than contrast. — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concept_edges[8] |
| relation | greater than a threshold —contrasted-with→ different to zero — Treat tests whether fold change exceeds a threshold rather than merely differs from zero. — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concept_edges[8] |
| weight | greater than a threshold —contrasted-with→ different to zero — 0.95 — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concept_edges[8] |
| note | Gene set enrichment analysis —is-a→ hybrid approach — Explicit categorical classification. — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | concept_edges[9] |
| relation | Gene set enrichment analysis —is-a→ hybrid approach — GSEA is a hybrid competitive and permutation approach. — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | concept_edges[9] |
| outcome | control — held — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | controls[0] |
| outcome | control — held — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | controls[1] |
| outcome | control — held — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | controls[2] |
| relation | limma → information borrowing — limma uses information borrowing to overcome small sample sizes. — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | unhoused_relations[0] |
| weight | limma → information borrowing — 0.95 — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | unhoused_relations[0] |
| relation | parametric empirical Bayes → class of statistical methods — A class of statistical methods is called parametric empirical Bayes. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | unhoused_relations[1] |
| relation | global parameters → gene-wise models — Global parameters link separate gene-wise models. — “A separate model is fitted for each gene, but the gene-wise models can be linked by global parameters or global hyper-parameters. The use of global parameters is a simple means of sharing information between genes that can be used even for the smallest experiments, because the global parameters can be estimated from the entire data set involving all the genes at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same gene, or correlations between related RNA samples, or variations in quality between the RNA samples.” | unhoused_relations[2] |
| relation | parametric empirical Bayes → residual variances — Empirical Bayes borrows strength across genes to moderate residual variances. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | unhoused_relations[3] |
| weight | parametric empirical Bayes → residual variances — 0.95 — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | unhoused_relations[3] |
| relation | voom function → precision weights — Voom converts the mean-variance trend into precision weights. — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | unhoused_relations[4] |
| relation | normalization → systematic effects — Normalization removes systematic technical effects unrelated to biology. — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | unhoused_relations[5] |
| relation | pre-processing → missing values or inflated variances — Pre-processing is designed to avoid missing values and inflated variances. — “Microarray expression data are measured as intensities, which need to be background corrected and normalized before any statistical analysis can be conducted. limma includes a range of background correction and normalization procedures suitable for different types of DNA microarrays or protein arrays. Notable are the maximum likelihood implementation of the normal-exponential convolution model for background correction ( 22 ) and the implementation of loess curves and normalization using quantitative weights. The guiding principle in the pre-processing steps is to preserve information, avoiding missing values or inflated variances ( 23 ). Normalized intensities are offset from zero before transforming to the log-scale to avoid missing values or large variances. Offsets in a range of moderate values have been shown to achieve an effective compromise between noise and bias ( 24 ).” | unhoused_relations[6] |
| weight | pre-processing → missing values or inflated variances — 0.95 — “Microarray expression data are measured as intensities, which need to be background corrected and normalized before any statistical analysis can be conducted. limma includes a range of background correction and normalization procedures suitable for different types of DNA microarrays or protein arrays. Notable are the maximum likelihood implementation of the normal-exponential convolution model for background correction ( 22 ) and the implementation of loess curves and normalization using quantitative weights. The guiding principle in the pre-processing steps is to preserve information, avoiding missing values or inflated variances ( 23 ). Normalized intensities are offset from zero before transforming to the log-scale to avoid missing values or large variances. Offsets in a range of moderate values have been shown to achieve an effective compromise between noise and bias ( 24 ).” | unhoused_relations[6] |
| relation | multi-dimensional scaling → batch effects — MDS helps identify unexpected batch effects before modelling. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | unhoused_relations[7] |
| weight | multi-dimensional scaling → batch effects — 0.95 — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | unhoused_relations[7] |
| relation | robustified shrinkage strategy → false positives — Robust shrinkage reduces false positives from spuriously small variances. — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | unhoused_relations[8] |
| weight | robustified shrinkage strategy → false positives — 0.95 — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | unhoused_relations[8] |
| relation | genas function → biological correlation — Genas estimates the true biological correlation between contrasts. — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | unhoused_relations[9] |
| relation | rank-based tests → optimistic P -values — Conditional-independence assumptions in rank tests yield optimistic P-values. — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | unhoused_relations[10] |
| weight | rank-based tests → optimistic P -values — 0.95 — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | unhoused_relations[10] |
| relation | camera function → optimistic P -values — Camera avoids optimistic P-values by adjusting for inter-gene dependence. — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | unhoused_relations[11] |
| weight | camera function → optimistic P -values — 0.95 — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | unhoused_relations[11] |
| relation | rotation tests → smoothed version of permutation — Rotation can be viewed as a smoothed version of permutation suitable for linear models. — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | unhoused_relations[12] |
| weight | rotation tests → smoothed version of permutation — 0.7 — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | unhoused_relations[12] |
| relation | matrix → information — Coercion of limma classes to matrix entails loss of information. — “Furthermore all of these classes may be coerced to be of class matrix using as.matrix, although this entails loss of information. Fitted model objects of class ‘MArrayLM’ can be coerced to class data.frame using as.data.frame in R. The first five classes belong to the virtual class ‘LargeDataObject’ for which a show method is defined to display the leading rows of a large vector, matrix or data.frame.” | unhoused_relations[13] |
| relation | C code → computationally intensive steps — C code speeds up computationally intensive steps in limma. — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | unhoused_relations[14] |
| weight | C code → computationally intensive steps — 0.95 — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | unhoused_relations[14] |
| relation | voom methodology → RNA-seq data — Voom unlocks most limma analysis methods for RNA-seq data. — “Although originally developed with microarray data in mind, the development of the voom methodology unlocks the majority of analysis methods for use on RNA-seq data, such as random effects modelling and gene set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seq) and experimental design.” | unhoused_relations[15] |
| weight | voom methodology → RNA-seq data — 0.95 — “Although originally developed with microarray data in mind, the development of the voom methodology unlocks the majority of analysis methods for use on RNA-seq data, such as random effects modelling and gene set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seq) and experimental design.” | unhoused_relations[15] |
| no type for this | limma → information borrowing — limma uses information borrowing to overcome small sample sizes. — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | unhoused_relations[0] |
| no type for this | parametric empirical Bayes → class of statistical methods — A class of statistical methods is called parametric empirical Bayes. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | unhoused_relations[1] |
| no type for this | global parameters → gene-wise models — Global parameters link separate gene-wise models. — “A separate model is fitted for each gene, but the gene-wise models can be linked by global parameters or global hyper-parameters. The use of global parameters is a simple means of sharing information between genes that can be used even for the smallest experiments, because the global parameters can be estimated from the entire data set involving all the genes at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same gene, or correlations between related RNA samples, or variations in quality between the RNA samples.” | unhoused_relations[2] |
| no type for this | parametric empirical Bayes → residual variances — Empirical Bayes borrows strength across genes to moderate residual variances. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | unhoused_relations[3] |
| no type for this | voom function → precision weights — Voom converts the mean-variance trend into precision weights. — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | unhoused_relations[4] |
| no type for this | normalization → systematic effects — Normalization removes systematic technical effects unrelated to biology. — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | unhoused_relations[5] |
| no type for this | pre-processing → missing values or inflated variances — Pre-processing is designed to avoid missing values and inflated variances. — “Microarray expression data are measured as intensities, which need to be background corrected and normalized before any statistical analysis can be conducted. limma includes a range of background correction and normalization procedures suitable for different types of DNA microarrays or protein arrays. Notable are the maximum likelihood implementation of the normal-exponential convolution model for background correction ( 22 ) and the implementation of loess curves and normalization using quantitative weights. The guiding principle in the pre-processing steps is to preserve information, avoiding missing values or inflated variances ( 23 ). Normalized intensities are offset from zero before transforming to the log-scale to avoid missing values or large variances. Offsets in a range of moderate values have been shown to achieve an effective compromise between noise and bias ( 24 ).” | unhoused_relations[6] |
| no type for this | multi-dimensional scaling → batch effects — MDS helps identify unexpected batch effects before modelling. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | unhoused_relations[7] |
| no type for this | robustified shrinkage strategy → false positives — Robust shrinkage reduces false positives from spuriously small variances. — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | unhoused_relations[8] |
| no type for this | genas function → biological correlation — Genas estimates the true biological correlation between contrasts. — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | unhoused_relations[9] |
| no type for this | rank-based tests → optimistic P -values — Conditional-independence assumptions in rank tests yield optimistic P-values. — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | unhoused_relations[10] |
| no type for this | camera function → optimistic P -values — Camera avoids optimistic P-values by adjusting for inter-gene dependence. — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | unhoused_relations[11] |
| no type for this | rotation tests → smoothed version of permutation — Rotation can be viewed as a smoothed version of permutation suitable for linear models. — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | unhoused_relations[12] |
| no type for this | matrix → information — Coercion of limma classes to matrix entails loss of information. — “Furthermore all of these classes may be coerced to be of class matrix using as.matrix, although this entails loss of information. Fitted model objects of class ‘MArrayLM’ can be coerced to class data.frame using as.data.frame in R. The first five classes belong to the virtual class ‘LargeDataObject’ for which a show method is defined to display the leading rows of a large vector, matrix or data.frame.” | unhoused_relations[13] |
| no type for this | C code → computationally intensive steps — C code speeds up computationally intensive steps in limma. — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | unhoused_relations[14] |
| no type for this | voom methodology → RNA-seq data — Voom unlocks most limma analysis methods for RNA-seq data. — “Although originally developed with microarray data in mind, the development of the voom methodology unlocks the majority of analysis methods for use on RNA-seq data, such as random effects modelling and gene set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seq) and experimental design.” | unhoused_relations[15] |
| limitation | Gene-expression studies commonly combine complex designs with few biological replicates. — “Despite the complexity, gene expression studies often involve only a small number of biological replicates. The small but complex nature of gene expression studies poses challenging statistical problems and motivates the use of a number of specialized statistical techniques in order to get the most out of each data set. We have developed the limma software over the past decade to provide a framework for analysing gene expression experiments from beginning to end in a flexible and statistically rigorous way.” | limitations[0] |
| limitation | Parametric expression models can never be perfectly correct. — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | limitations[1] |
| limitation | Permutation has low power with few replicates and cannot represent complex null hypotheses, correlation, blocking, or unequal precision. — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | limitations[2] |
| limitation | RPKM input prevents limma from estimating the appropriate mean-variance relationship; raw counts are required before voom. — “Raw read counts are assembled outside limma using tools such as featureCounts ( 29 ), HTSeq-counts ( 30 ) or RSEM ( 31 ). The authors of this article find the Subread ( 32 ) and featureCounts pipeline particularly convenient because it is fast, accurate ( 8 ) and can be run from the R prompt using the Rsubread package. The data input to limma should be counts, rather than popular expression summaries such as reads-per-kilobase-per-million (RPKM), so that limma can estimate the appropriate mean-variance relationship. The voom output can be converted to RPKM values for convenience of interpretation, by subtracting log-gene-lengths, but this should be done after running voom rather than before.” | limitations[3] |
| limitation | Simple foreground-minus-background subtraction is too heavy-handed. — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | limitations[4] |
| limitation | Rank-based gene-set tests assume conditional gene independence and therefore produce optimistic P-values. — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | limitations[5] |
| limitation | Camera's dependence correction reduces optimistic P-values but also lowers statistical power. — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | limitations[6] |
| limitation | Gene-set workflows must reconcile symbols and annotations across studies because symbols change over time. — “The use of gene sets require that gene symbols and annotation be matched between different databases and studies. Gene symbols change over time, so limma includes the alias2Symbol and alias2SymbolTable functions to map gene symbol aliases to current official gene symbols. The id2indices function matches gene identifiers in a collection of gene sets to those in the expression matrix, in a format suitable for input into camera, mroast or romer.” | limitations[7] |
| limitation | Coercing limma data classes to matrices entails information loss. — “Furthermore all of these classes may be coerced to be of class matrix using as.matrix, although this entails loss of information. Fitted model objects of class ‘MArrayLM’ can be coerced to class data.frame using as.data.frame in R. The first five classes belong to the virtual class ‘LargeDataObject’ for which a show method is defined to display the leading rows of a large vector, matrix or data.frame.” | limitations[8] |
| limitation | The appropriate method combination remains dependent on biological question, platform, and design. — “Although originally developed with microarray data in mind, the development of the voom methodology unlocks the majority of analysis methods for use on RNA-seq data, such as random effects modelling and gene set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seq) and experimental design.” | limitations[9] |
Content tags
content-observations:
finding · quantity · variable · definition · key-term — grounded in the full text| finding | limma provides an integrated gene-expression analysis solution across complex designs and small samples. — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | content_tags.findings[0] |
| finding | One downstream pipeline can analyse both RNA-seq and microarray data after preprocessing. — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | content_tags.findings[1] |
| finding | Whole-experiment linear models share information and accommodate sample correlations. — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | content_tags.findings[2] |
| finding | Empirical Bayes moderation stabilizes gene-wise variance inference when replicates are few. — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | content_tags.findings[3] |
| finding | Flexible robust empirical Bayes modelling improves statistical power and accuracy. — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | content_tags.findings[4] |
| finding | Quantitative weights increase power while avoiding ad hoc observation filtering. — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | content_tags.findings[5] |
| finding | Voom converts the mean-variance trend into precision weights for log-transformed RNA-seq counts. — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | content_tags.findings[6] |
| finding | The voom pipeline performs comparably to leading negative-binomial packages with greater large-data speed and reliability. — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | content_tags.findings[7] |
| finding | Moderate preprocessing offsets balance noise and bias while avoiding missing values and inflated variances. — “Microarray expression data are measured as intensities, which need to be background corrected and normalized before any statistical analysis can be conducted. limma includes a range of background correction and normalization procedures suitable for different types of DNA microarrays or protein arrays. Notable are the maximum likelihood implementation of the normal-exponential convolution model for background correction ( 22 ) and the implementation of loess curves and normalization using quantitative weights. The guiding principle in the pre-processing steps is to preserve information, avoiding missing values or inflated variances ( 23 ). Normalized intensities are offset from zero before transforming to the log-scale to avoid missing values or large variances. Offsets in a range of moderate values have been shown to achieve an effective compromise between noise and bias ( 24 ).” | content_tags.findings[8] |
| finding | Permutation inference is poorly suited to small complex correlated experiments with unequal precision. — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | content_tags.findings[9] |
| finding | MDS leading-fold-change distances expose sample structure and batch effects before model construction. — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | content_tags.findings[10] |
| finding | Quality weights handle reliability variation gradually and improve DE power. — “limma is the only package that allows variations in quality to be handled in a graduated way via quantitative weights. Both observation-level ( 2 , 10 , 43 ) and sample-specific weights ( 11 ) can be used in an analysis. For microarray data, the arrayWeights function estimates relative array variances, which are converted to weights that can be used in the linear model analysis to down-weight observations from less reliable arrays. Probe and array weights can be easily combined by multiplying them together and, when used appropriately, have been demonstrated to increase power to detect DE ( 11 , 43 ).” | content_tags.findings[11] |
| finding | Robust shrinkage reduces false positives from tiny variances without over-shrinking unusually large variances. — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | content_tags.findings[12] |
| finding | Treat can prioritize fold changes that are biologically as well as statistically significant. — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | content_tags.findings[13] |
| finding | DiffSplice makes large-scale differential-exon-usage surveys faster than alternatives. — “The linear model framework of limma is extended to test very easily for differential splicing events when exon-level expression data are available. The data can be either from an exon microarray or from RNA-seq data summarized at the exon level. In either case, the approach is based on fitting linear models to the exon-level expression data. The approach can relate differential exon usage to continuous as well as categorical predictors or to any contrast in a linear model. The test is conducted by the diffSplice function and results are displayed using plotSplice and topSplice. The plotExons function is also useful for exploring exon expression for individual genes. This approach is considerably faster than alternative approaches to differential splicing, making large-scale surveys of differential exon usage feasible.” | content_tags.findings[14] |
| finding | Genas estimates error-free biological correlation and can reveal shared pathways missed by significance-cutoff overlap. — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | content_tags.findings[15] |
| finding | Camera corrects optimistic gene-set P-values for inter-gene dependence at a cost in power. — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | content_tags.findings[16] |
| finding | Roast preserves directional prior annotation and works in small complex experiments. — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | content_tags.findings[17] |
| finding | Most limma estimation procedures finish in seconds and virtually all in under a minute on a desktop. — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | content_tags.findings[18] |
| finding | More than 120 Bioconductor packages used limma as of March 2014. — “The limma software is freely available online as part of the Bioconductor project ( http://www.bioconductor.org ). More than 120 other Bioconductor packages make use of limma (as of March 2014). The limma package is used as a building block or as the underlying computational engine by a number of software projects designed to provide user-interfaces for gene expression data analysis including limmaGUI ( 66 ), affylmGUI ( 67 ), WebArray ( 68 ), RACE ( 69 ), CarmaWEB ( 70 ), Goulphar ( 71 ), MAGMA ( 72 ), Asterias ( 73 ), GenePattern ( 74 ), GEO2R ( http://www.ncbi.nlm.nih.gov/geo/geo2r ), the EBI expression atlas ( 75 ), Guide ( 76 ) and Degust ( http://www.vicbioinformatics.com/degust ).” | content_tags.findings[19] |
| quantity | small number — typical number of biological replicates — “Despite the complexity, gene expression studies often involve only a small number of biological replicates. The small but complex nature of gene expression studies poses challenging statistical problems and motivates the use of a number of specialized statistical techniques in order to get the most out of each data set. We have developed the limma software over the past decade to provide a framework for analysing gene expression experiments from beginning to end in a flexible and statistically rigorous way.” | content_tags.quantities[0] |
| quantity | two important directions — RNA-seq/splicing support and higher-order signature analysis — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | content_tags.quantities[1] |
| quantity | two group comparisons — simplest experimental-design case — “Gene expression technologies are used frequently in molecular biology research to gain a snapshot of transcriptional activity in different tissues or populations of cells. These profiles are then compared to identify gene expression changes associated with a treatment condition or phenotype of interest. Gene expression studies may be randomized designed experiments in which a biological system is perturbed, for example by a gene knock-out or by applying a specified stressor. Such experiments are amongst the most powerful tools in functional genomics, providing insights into normal cellular processes as well as disease pathogenesis. Or they may be observational studies in which different phenotypes are compared, diseased and normal tissue for example or cells from different populations. Such studies are common in cancer research and in the study of cell development. In either case, the study design can range from simple two group comparisons to complex set-ups with several experimental factors varying over multiple levels. Researchers might be interested for example in whether a particular gene facilitates or blocks the action of a particular drug, in which case knock-down and wild-type samples both with and without drug treatment would be profiled. Observational studies may involve multiple batch effects and covariates that must be accounted for in the analysis.” | content_tags.quantities[2] |
| quantity | 0.05 — FDR threshold used in the displayed DE example — “Figure 4 shows example DE summary plots. Panels (A) and (B) were generated using the two-colour microarray data from GEO series GSE2593 . Intensities were background corrected and normalized as previously described ( 13 ). Panel (A) shows a volcano plot for the comparison of samples with RUNX1 over-expressed versus wild-type samples, while panel (B) shows a Venn diagram of differentially expressed probes for each of the three over-expressed genes versus wild-type. Probes with false discovery rate less than 0.05 were considered to be differentially expressed. Panel (C) uses RNA-seq data from GEO series GSE52870 . The data were analysed as described in Figure 5 of Liu et al . ( 7 ).” | content_tags.quantities[3] |
| quantity | over 120 — other Bioconductor packages using limma as of March 2014 — “The limma software is freely available online as part of the Bioconductor project ( http://www.bioconductor.org ). More than 120 other Bioconductor packages make use of limma (as of March 2014). The limma package is used as a building block or as the underlying computational engine by a number of software projects designed to provide user-interfaces for gene expression data analysis including limmaGUI ( 66 ), affylmGUI ( 67 ), WebArray ( 68 ), RACE ( 69 ), CarmaWEB ( 70 ), Goulphar ( 71 ), MAGMA ( 72 ), Asterias ( 73 ), GenePattern ( 74 ), GEO2R ( http://www.ncbi.nlm.nih.gov/geo/geo2r ), the EBI expression atlas ( 75 ), Guide ( 76 ) and Degust ( http://www.vicbioinformatics.com/degust ).” | content_tags.quantities[4] |
| quantity | over 120 pages — length of the limma user guide — “Third, the package comes with an extensive user's guide of over 120 pages, available from the drop-down menu in Windows or alternatively launched by the limmaUsersGuide command. The user's guide gives detailed advice on how to analyse a variety of common study designs. It also includes 10 fully worked case studies for which full data and code are provided.” | content_tags.quantities[5] |
| quantity | 10 — fully worked case studies with data and code — “Third, the package comes with an extensive user's guide of over 120 pages, available from the drop-down menu in Windows or alternatively launched by the limmaUsersGuide command. The user's guide gives detailed advice on how to analyse a variety of common study designs. It also includes 10 fully worked case studies for which full data and code are provided.” | content_tags.quantities[6] |
| quantity | a few seconds — typical estimation time on a standard desktop — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | content_tags.quantities[7] |
| quantity | less than a minute — upper time for virtually all estimation procedures — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | content_tags.quantities[8] |
| quantity | linear — memory scaling in genes and samples — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | content_tags.quantities[9] |
| variable | gene expression values (row-wise response matrix) — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | content_tags.variables[0] |
| variable | genes (matrix rows and units of gene-wise modelling) — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | content_tags.variables[1] |
| variable | RNA samples (matrix columns and experimental observations) — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | content_tags.variables[2] |
| variable | residual variances (gene-wise quantities moderated by empirical Bayes) — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | content_tags.variables[3] |
| variable | mean-variance trend (global abundance-precision relationship) — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | content_tags.variables[4] |
| variable | precision weights (voom-derived observation reliability) — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | content_tags.variables[5] |
| variable | read counts (RNA-seq input matrix values) — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | content_tags.variables[6] |
| variable | logCPM (log2 counts-per-million transformed response) — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | content_tags.variables[7] |
| variable | batch effects (nuisance variation adjusted in design matrices) — “Linear models permit very general analyses. Researchers can adjust for the effects of multiple experimental factors or can adjust for batch effects. The linear model might include time course effects or regression splines. The linear model could even include the expression values themselves of one or more genes as covariates, allowing researchers to test for inter-gene dependencies. Linear models allow researchers to test very flexible hypotheses, not just simple comparisons between groups but also interaction effects or more complex customized comparisons.” | content_tags.variables[8] |
| variable | leading fold change (sample-distance quantity in MDS) — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | content_tags.variables[9] |
| variable | false discovery rate (multiple-testing error criterion) — “limma provides a number of options to adjust tests for multiple testing. Users can control either the family-wise type I error rate or the false discovery rate ( 46 ). As well as the usual control for multiple testing across multiple genes, limma is the only software package to provide methods for error rate control across multiple contrasts and genes simultaneously. For individual tests, multiple testing can be applied using the topTable function. The decideTests function gives access to the full range of options.” | content_tags.variables[10] |
| variable | inter-gene correlation (camera input used for variance inflation) — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | content_tags.variables[11] |
| definition | limma — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | content_tags.definitions[0] |
| definition | parametric empirical Bayes — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | content_tags.definitions[2] |
| definition | posterior variance estimator — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | content_tags.definitions[3] |
| definition | quantitative weights — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | content_tags.definitions[4] |
| definition | voom — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | content_tags.definitions[5] |
| definition | background correction — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | content_tags.definitions[6] |
| definition | normalization — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | content_tags.definitions[7] |
| definition | leading fold change — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | content_tags.definitions[8] |
| definition | biological correlation — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | content_tags.definitions[9] |
| definition | gene set analyses — “Gene set analyses assess the overall significance of a set of co-regulated genes. Each gene set is chosen to represent a particular molecular pathway or some other biological process of interest. Gene sets are defined by gene annotation external to the current expression study, for example from Gene Ontology (GO) database ( 52 ) or from previous expression studies. For gene sets defined by previous studies, the genes may optionally be annotated with the direction and magnitude of expression changes in the earlier experiment. In this way, a gene set may contain genes both positively and negatively associated with the molecular pathway that it represents.” | content_tags.definitions[10] |
| definition | rotation test — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | content_tags.definitions[11] |
| definition | GSEA — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | content_tags.definitions[12] |
| definition | RGList — “For raw intensity data, the classes ‘RGList’ and ‘EListRaw’ are used to store two-colour and single-channel data, respectively. These objects are often created using the function read.maimages and contain the raw values from the image analysis output files along with probe annotation information.” | content_tags.definitions[13] |
| definition | EList — “Normalized data are stored in ‘MAList’ or ‘EList’ objects. Normalized two-colour data are converted from red and green intensities, R and G , into M and A -values, which hold the log-ratio and average log-intensity values for each spot. Single channel data are background corrected and log 2 transformed and stored in an ‘EList’ object. For RNA-seq data, the voom transformed matrix of gene/exon counts is also stored in an ‘EList’ object.” | content_tags.definitions[14] |
| definition | MArrayLM — “The next major classes store output from a DE analysis. ‘MArrayLM’ objects store the result of fitting gene-wise linear models to the normalized intensities or log-ratios. Objects of this class are created by the lmFit and eBayes functions. After running decideTests, an object of class ‘TestResults’ stores the results of testing a set of contrasts equal to zero for each probe/gene.” | content_tags.definitions[15] |
Concepts
terms this paper carries, each with the sentence it was cut from — the layer that joins one paper to another. Terms with no quote yet are on CONCEPTS.html.
| concept | limma — not defined in this paper's own words; integrated R/Bioconductor expression-analysis package — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concepts[0] |
| concept | R/Bioconductor software package — not defined in this paper's own words; software class assigned to limma — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concepts[1] |
| concept | gene expression experiments — not defined in this paper's own words; experiments generating transcriptional profiles — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concepts[2] |
| concept | information borrowing — not defined in this paper's own words; cross-feature estimation used for small samples — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concepts[3] |
| concept | differential expression — not defined in this paper's own words; gene-level change between conditions — “limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. It contains rich features for handling complex experimental designs and for information borrowing to overcome the problem of small sample sizes. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limma have been significantly expanded in two important directions. First, the package can now perform both differential expression and differential splicing analyses of RNA sequencing (RNA-seq) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures. This provides enhanced possibilities for biological interpretation of gene expression differences. This article reviews the philosophy and design of the limma package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.” | concepts[4] |
| concept | differential splicing — not defined in this paper's own words; condition-associated exon-usage change — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | concepts[5] |
| concept | RNA sequencing — not defined in this paper's own words; count-based expression technology — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | concepts[6] |
| concept | RNA-seq — not defined in this paper's own words; abbreviation for RNA sequencing — “Recently, the capabilities of limma have expanded significantly in two important directions. First, the package can now perform both differential expression (DE) and differential splicing analyses of RNA sequencing (RNA-seq) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seq as well. These capabilities allow users to analyse both RNA-seq and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genes or in terms of higher-order expression signatures ( 7 ). This provides enhanced possibilities for biological interpretation of gene expression differences.” | concepts[7] |
| concept | Bioconductor — not defined in this paper's own words; R open-source statistical-genomics project — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | concepts[8] |
| concept | R-based open-source software development project — not defined in this paper's own words; class assigned to Bioconductor — “The limma package is a core component of Bioconductor, an R-based open-source software development project in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalization, the same analysis pipeline is used for data from all technologies.” | concepts[9] |
| concept | linear model — not defined in this paper's own words; row-wise model for expression features — “limma integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrix of expression values, where each row represents a gene, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sample. On the one hand, it fits a linear model to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genes and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of gene expression signatures. Figure 1 depicts the linear model and highlights the statistical principles employed in a typical limma analysis.” | concepts[10] |
| concept | linear models — not defined in this paper's own words; whole-experiment expression models — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | concepts[11] |
| concept | piece-meal comparisons between pairs of treatments — not defined in this paper's own words; alternative avoided by integrated modelling — “The hallmark of the limma approach is the use of linear models to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatments. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.” | concepts[12] |
| concept | global parameters — not defined in this paper's own words; shared estimates linking gene-wise models — “A separate model is fitted for each gene, but the gene-wise models can be linked by global parameters or global hyper-parameters. The use of global parameters is a simple means of sharing information between genes that can be used even for the smallest experiments, because the global parameters can be estimated from the entire data set involving all the genes at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same gene, or correlations between related RNA samples, or variations in quality between the RNA samples.” | concepts[13] |
| concept | global hyper-parameters — not defined in this paper's own words; shared higher-level estimates — “A separate model is fitted for each gene, but the gene-wise models can be linked by global parameters or global hyper-parameters. The use of global parameters is a simple means of sharing information between genes that can be used even for the smallest experiments, because the global parameters can be estimated from the entire data set involving all the genes at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same gene, or correlations between related RNA samples, or variations in quality between the RNA samples.” | concepts[14] |
| concept | parametric empirical Bayes — not defined in this paper's own words; parallel-model information-borrowing class — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | concepts[15] |
| concept | class of statistical methods — not defined in this paper's own words; class assigned to parametric empirical Bayes — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | concepts[16] |
| concept | empirical Bayes procedures — not defined in this paper's own words; flexible mean-variance and robust variance modelling procedures — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concepts[17] |
| concept | residual variances — not defined in this paper's own words; gene-wise variability estimates under moderation — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | concepts[18] |
| concept | posterior variance estimators — not defined in this paper's own words; empirical-Bayes small-sample variance estimates — “The highly parallel nature of gene expression experiments lends itself to a particular class of statistical methods, called parametric empirical Bayes, that borrow information between genes in a dynamic way ( 14 , 15 ). The fact that the same linear model is fitted to each gene allows us to borrow strength between genes in order to moderate the residual variances ( 16 ). The estimated variance for each gene then becomes a compromise between the gene-wise estimator, obtained from the data for that gene alone, and the global variability across all genes, estimated by pooling the ensemble of all genes. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limma package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimators ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.” | concepts[19] |
| concept | mean-variance trend — not defined in this paper's own words; abundance-dependent reliability pattern — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concepts[20] |
| concept | statistical power and accuracy — not defined in this paper's own words; outcomes improved by flexible global modelling — “In recent years, the empirical Bayes procedures of limma have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trend ( 10 , 17 , 18 ). This is important because many gene expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genes. This allows a sophisticated robust empirical Bayes procedure in which hyper-variable genes are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracy by improving the modelling of the global characteristics of the data in a more flexible way.” | concepts[21] |
| concept | quantitative weights — not defined in this paper's own words; graduated reliability modifiers — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concepts[22] |
| concept | power to detect differentially expressed genes — not defined in this paper's own words; detection capability increased by weights — “Another unique feature of limma is the ability to incorporate quantitative weights into all levels of the statistical analysis, from normalization to linear modelling and gene set testing. Weights can be applied to genes or to RNA samples or to individual expression values. Weights can be used to give more emphasis to control probes during normalization, or can be used to down-weight measurements or samples that are less reliable in a gene expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genes, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).” | concepts[23] |
| concept | voom function — not defined in this paper's own words; count mean-variance to precision-weight transformation — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | concepts[24] |
| concept | precision weights — not defined in this paper's own words; observation reliability supplied to linear models — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | concepts[25] |
| concept | negative binomial-based software packages — not defined in this paper's own words; RNA-seq comparators for voom — “All the downstream analysis features of limma are available for RNA-seq and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seq data require specialized software based on the negative binomial or similar distributions ( 20 ). limma however is able to analyse RNA-seq read counts with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trend is converted by the voom function into precision weights, which are incorporated into the analysis of log-transformed RNA-seq counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packages but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seq and microarray experiments within limma . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.” | concepts[26] |
| concept | heteroscedasticity — not defined in this paper's own words; unequal expression variability — “Expression values often show some degree of heteroscedasticity, either because there is a relationship between abundance and measurement precision, or because some treatment conditions are more heterogeneous than others. For example, tumours might be more variable than normal tissue. Concern for such effects has prompted some researchers to filter out low intensity observations or to use Welch's t -test for DE between two groups instead of classical pooled t -tests. The use of weights and the ability to model global parameters allow limma to incorporate unequal variances in a number of ways. One way is through estimating a mean-variance trend, which can either be incorporated into the empirical Bayes procedure as mentioned above or used to generate observation weights ( 10 ). A recent development is the ability to estimate precision weights associated with treatment groups or more generally with any given set of covariates. More generally again, the mean-variance trend can be estimated in a treatment-specific way, combining the two types of heteroscedasticity mentioned above. These approaches allow limma to model unequal variances even for experiments with a small number of RNA samples. Importantly, they accommodate unequal variances without compromising the linear modelling and empirical Bayes framework of the package.” | concepts[27] |
| concept | expression signature — not defined in this paper's own words; gene-set representation of a process or cell type — “In recent years, the linear modelling capabilities of limma have been extended to higher-level expression signature analyses involving co-regulated sets of genes. The idea is to use a set of genes, together with their log-fold-changes, to represent the transcriptional signature of a biological process or cell type. One way that this is done is by rotation technology, which permits statistical significance to be tested for sets of genes for any linear model contrast ( 55 ). A particular feature of rotation tests is the ability to incorporate prior information about the direction and strength with which each gene is expected to contribute to the statistical signature. In this way, limma provides a uniquely flexible means to relate new expression data sets to previous results collated from earlier experiments, taking into account for example the fold-change and direction of change for each gene in the earlier experiment.” | concepts[28] |
| concept | background correction — not defined in this paper's own words; removal of ambient probe intensity — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | concepts[29] |
| concept | Removing this non-specific signal — not defined in this paper's own words; operation named background correction — “When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signal from the foreground intensity of each probe is called background correction and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limma backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.” | concepts[30] |
| concept | normalization — not defined in this paper's own words; technical-effect removal and scale alignment — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | concepts[31] |
| concept | systematic effects — not defined in this paper's own words; assay differences removed by normalization — “Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalization is to remove systematic effects due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalization methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalization or cyclic loess ( 28 ). Both functions provide a range of different normalization methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalization methods for two-colour arrays ( 36 , 37 ). limma is the only software to allow the use of quantitative weights in loess normalization ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genes. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genes may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalization).” | concepts[32] |
| concept | Bland–Altman plot — not defined in this paper's own words; difference-versus-mean agreement plot — “Measuring expression in multiple RNA samples produces columns of correlated expression values, which are highly correlated because they are measured on the same set of genes or genomic features. It has long been established in the biomedical literature that the level of agreement between correlated variables can be usefully examined by plotting differences versus means. Such a plot is called a Bland–Altman plot ( 25 ) or a Tukey mean-difference plot ( 26 ). Indeed the concept of DE can be viewed as a measure of disagreement between expression measures for the same genes in different samples. Mean-difference plots were introduced to the two-colour microarray literature by Dudoit et al . ( 27 ) and to the single-channel literature by Bolstad et al . ( 28 ), who called them MA-plots. limma generalized the concept of an MA-plot in two ways. First, the idea was extended to apply to sets of single-channel expression values. In this case, the plot is used to compare each sample to the average of all other samples. A virtual array is constructed by averaging the log-expression value for all the samples other than the sample of interest, and then a mean-difference plot is made between the single array and the virtual array. Second, the idea was extended to apply to the fitted model objects. In this case, the plot compares the log-fold-changes for a chosen contrast versus the average log-expression values of each gene across all samples. In effect, this plots a coefficient of the linear model versus an overall mean intercept parameter. These ideas were part of the original limma package submitted to Bioconductor in 2003.” | concepts[33] |
| concept | MA-plots — not defined in this paper's own words; microarray mean-difference plots — “Measuring expression in multiple RNA samples produces columns of correlated expression values, which are highly correlated because they are measured on the same set of genes or genomic features. It has long been established in the biomedical literature that the level of agreement between correlated variables can be usefully examined by plotting differences versus means. Such a plot is called a Bland–Altman plot ( 25 ) or a Tukey mean-difference plot ( 26 ). Indeed the concept of DE can be viewed as a measure of disagreement between expression measures for the same genes in different samples. Mean-difference plots were introduced to the two-colour microarray literature by Dudoit et al . ( 27 ) and to the single-channel literature by Bolstad et al . ( 28 ), who called them MA-plots. limma generalized the concept of an MA-plot in two ways. First, the idea was extended to apply to sets of single-channel expression values. In this case, the plot is used to compare each sample to the average of all other samples. A virtual array is constructed by averaging the log-expression value for all the samples other than the sample of interest, and then a mean-difference plot is made between the single array and the virtual array. Second, the idea was extended to apply to the fitted model objects. In this case, the plot compares the log-fold-changes for a chosen contrast versus the average log-expression values of each gene across all samples. In effect, this plots a coefficient of the linear model versus an overall mean intercept parameter. These ideas were part of the original limma package submitted to Bioconductor in 2003.” | concepts[34] |
| concept | permutation or re-sampling-based inference — not defined in this paper's own words; inference family excluded by limma — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | concepts[35] |
| concept | parametric distributions — not defined in this paper's own words; variance models preferred for small complex studies — “It is worth mentioning what limma does not do, which is permutation or re-sampling-based inference. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributions, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.” | concepts[36] |
| concept | read counts — not defined in this paper's own words; RNA-seq input values — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | concepts[37] |
| concept | log 2 counts per million — not defined in this paper's own words; voom-transformed count scale — “limma accepts RNA-seq data in the form of a matrix of read counts, with rows for genomic features and columns for RNA samples. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genes or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genes for simplicity of terminology. The read counts are processed by the voom function in limma to convert them into log 2 counts per million (logCPM) with associated precision weights. The logCPM values can be normalized between samples by the voom function or can be pre-normalized by adding normalization factors within edgeR .” | concepts[38] |
| concept | multi-dimensional scaling — not defined in this paper's own words; unsupervised sample-profile display — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concepts[39] |
| concept | leading fold change — not defined in this paper's own words; MDS distance between samples — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concepts[40] |
| concept | root-mean-square average of the log-fold-changes — not defined in this paper's own words; operational definition of leading fold change — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concepts[41] |
| concept | batch effects — not defined in this paper's own words; unexpected systematic sample patterns — “After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genes. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scaling to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold change’, which is defined as the root-mean-square average of the log-fold-changes for the genes best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effects, that should be adjusted for during the linear model analysis. This helps guide the construction of the design matrix used for the linear modelling below.” | concepts[42] |
| concept | consensus correlation — not defined in this paper's own words; shared intrablock correlation estimate — “limma includes a unique strategy for incorporating the fact that observations or samples may be correlated. The strategy is similar to fitting a random effects model, with the difference that all genes are constrained to share the same intrablock correlation. The duplicateCorrelation function is used to estimate the consensus correlation. The correlation structure is then incorporated into the linear model fit and hence into all tests for DE. Originally the idea was used to estimate the correlation between replicate copies of the same probe on a microarray ( 44 ). The correlation strategy preserves more information than simply averaging the replicate probe copies. More generally, the same idea is also used to model the correlation between related RNA samples, for example repeated measures on the same individual or RNA samples collected at the same time.” | concepts[43] |
| concept | robustified shrinkage strategy — not defined in this paper's own words; gene-specific variance moderation — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | concepts[44] |
| concept | false positives — not defined in this paper's own words; spurious DE calls reduced by moderation — “An empirical Bayes framework to borrow information between genes when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positives for genes with very small variances and improves power to detect DE for genes with larger variances. limma includes a robustified shrinkage strategy that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genes are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genes, whilst negating the effects of outliers.” | concepts[45] |
| concept | log 2 -fold-change — not defined in this paper's own words; effect-size quantity tested by treat — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concepts[46] |
| concept | greater than a threshold — not defined in this paper's own words; minimum effect-size alternative — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concepts[47] |
| concept | different to zero — not defined in this paper's own words; traditional null departure contrasted with a threshold — “When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-change is greater than a threshold rather than merely different to zero ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.” | concepts[48] |
| concept | false discovery rate — not defined in this paper's own words; multiple-testing error rate — “limma provides a number of options to adjust tests for multiple testing. Users can control either the family-wise type I error rate or the false discovery rate ( 46 ). As well as the usual control for multiple testing across multiple genes, limma is the only software package to provide methods for error rate control across multiple contrasts and genes simultaneously. For individual tests, multiple testing can be applied using the topTable function. The decideTests function gives access to the full range of options.” | concepts[49] |
| concept | propTrueNull function — not defined in this paper's own words; estimator of undetected true DE — “The propTrueNull function estimates the number of truly differentially expressed genes that remain to be identified. Mathematically, it estimates the proportion of true null hypotheses in a collection of hypothesis tests given a vector of P -values. In a gene expression study, it estimates the proportion of non-differentially expressed genes, out of all tested, for any contrast in the linear model. The function implements a number of different methods for estimating the proportion of true nulls, ranging from quick and simple to more computationally demanding. The default is based on averaging local false discovery rates across the P -values ( 19 ). Other methods are the histogram method of ( 48 , 49 ), the convex decreasing density estimate of ( 50 ) and a very simple estimate based on averaging the P -values.” | concepts[50] |
| concept | proportion of true null hypotheses — not defined in this paper's own words; quantity estimated from P-values — “The propTrueNull function estimates the number of truly differentially expressed genes that remain to be identified. Mathematically, it estimates the proportion of true null hypotheses in a collection of hypothesis tests given a vector of P -values. In a gene expression study, it estimates the proportion of non-differentially expressed genes, out of all tested, for any contrast in the linear model. The function implements a number of different methods for estimating the proportion of true nulls, ranging from quick and simple to more computationally demanding. The default is based on averaging local false discovery rates across the P -values ( 19 ). Other methods are the histogram method of ( 48 , 49 ), the convex decreasing density estimate of ( 50 ) and a very simple estimate based on averaging the P -values.” | concepts[51] |
| concept | genas function — not defined in this paper's own words; bias-adjusted expression-profile association estimator — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | concepts[52] |
| concept | biological correlation — not defined in this paper's own words; error-free correlation between contrast fold changes — “Gene expression experiments typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genes from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genes, has little chance of success in situations where power to detect differentially expressed genes is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas function ( 19 ) addresses these problems. It tests whether two different contrasts in a linear model affect the same genes in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlation between the log 2 -fold-changes of two different contrasts. By biological correlation, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limma . This method is particularly powerful for gaining insight into commonly affected gene pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genes from each contrast alone.” | concepts[53] |
| concept | gene set analyses — not defined in this paper's own words; tests of co-regulated gene groups — “Gene set analyses assess the overall significance of a set of co-regulated genes. Each gene set is chosen to represent a particular molecular pathway or some other biological process of interest. Gene sets are defined by gene annotation external to the current expression study, for example from Gene Ontology (GO) database ( 52 ) or from previous expression studies. For gene sets defined by previous studies, the genes may optionally be annotated with the direction and magnitude of expression changes in the earlier experiment. In this way, a gene set may contain genes both positively and negatively associated with the molecular pathway that it represents.” | concepts[54] |
| concept | molecular pathway — not defined in this paper's own words; biological process represented by a gene set — “Gene set analyses assess the overall significance of a set of co-regulated genes. Each gene set is chosen to represent a particular molecular pathway or some other biological process of interest. Gene sets are defined by gene annotation external to the current expression study, for example from Gene Ontology (GO) database ( 52 ) or from previous expression studies. For gene sets defined by previous studies, the genes may optionally be annotated with the direction and magnitude of expression changes in the earlier experiment. In this way, a gene set may contain genes both positively and negatively associated with the molecular pathway that it represents.” | concepts[55] |
| concept | rank-based tests — not defined in this paper's own words; threshold-free ordered-gene tests — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | concepts[56] |
| concept | optimistic P -values — not defined in this paper's own words; anti-conservative result of independence assumptions — “Unlike goana, the other gene set options do not require a significance threshold to be applied to identify differentially expressed genes. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based tests ( 13 ). These tests assess whether the specified set of genes is more highly ranked in an ordered list of all genes than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each gene is conditionally independent of other genes and hence give optimistic P -values ( 55 ).” | concepts[57] |
| concept | camera function — not defined in this paper's own words; competitive inter-gene-correlation-adjusted test — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | concepts[58] |
| concept | variance inflation factor — not defined in this paper's own words; camera dependence adjustment — “More sophisticated competitive tests that take into account dependence between the genes in the linear modelling framework are implemented in the camera function ( 56 ). Camera computes a variance inflation factor from the inter-gene correlation and uses this to adjust the variance of the summary statistics. This avoids optimistic P -values in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).” | concepts[59] |
| concept | rotation tests — not defined in this paper's own words; residual-space gene-set simulations — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | concepts[60] |
| concept | smoothed version of permutation — not defined in this paper's own words; interpretation assigned to rotation — “The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genes is differentially expressed ( 55 ). These functions use rotation tests, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutation that is suitable for linear models. Rotation occurs only in the residual space of the linear model in such a way that any coefficients in the linear model other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experiments. Of all gene set tests, roast has the unique ability to take into account directional annotation information about genes in the set. It is able to accommodate genes both positively and negatively associated with a specified pathway, as well as the magnitude of change of each gene. It is therefore especially useful for finding similarities in gene expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genes, for example exon-level expression analyses to test whether any exon of a given gene is differentially expressed.” | concepts[61] |
| concept | Gene set enrichment analysis — not defined in this paper's own words; database-scale co-regulated-set correlation — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | concepts[62] |
| concept | GSEA — not defined in this paper's own words; abbreviation for Gene set enrichment analysis — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | concepts[63] |
| concept | hybrid approach — not defined in this paper's own words; competitive-plus-permutation class assigned to GSEA — “Gene set enrichment analysis (GSEA) is an approach that correlates a large database of co-regulated gene sets with respect to a microarray or RNA-seq data set ( 62 , 63 ). GSEA is a hybrid approach: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limma implements a GSEA approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of gene sets and with any linear model. The limma authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).” | concepts[64] |
| concept | object-oriented paradigm — not defined in this paper's own words; limma's consistent user-interface design — “A simple but appropriate object-oriented paradigm provides users with a consistent analysis interface that is very easy from a user point of view. Many limma functions are generic or operate appropriately on objects of different classes. limma defines a number of classes that have been tailored to handle both microarray and RNA-seq data. The philosophy has been to define simple list-based data objects that can be easily explored and manipulated by users, in the same style as familiar, long-standing core functions in R such as lm and glm.” | concepts[65] |
| concept | RGList — not defined in this paper's own words; raw two-colour data class — “For raw intensity data, the classes ‘RGList’ and ‘EListRaw’ are used to store two-colour and single-channel data, respectively. These objects are often created using the function read.maimages and contain the raw values from the image analysis output files along with probe annotation information.” | concepts[66] |
| concept | EListRaw — not defined in this paper's own words; raw single-channel data class — “For raw intensity data, the classes ‘RGList’ and ‘EListRaw’ are used to store two-colour and single-channel data, respectively. These objects are often created using the function read.maimages and contain the raw values from the image analysis output files along with probe annotation information.” | concepts[67] |
| concept | MAList — not defined in this paper's own words; normalized two-colour data class — “Normalized data are stored in ‘MAList’ or ‘EList’ objects. Normalized two-colour data are converted from red and green intensities, R and G , into M and A -values, which hold the log-ratio and average log-intensity values for each spot. Single channel data are background corrected and log 2 transformed and stored in an ‘EList’ object. For RNA-seq data, the voom transformed matrix of gene/exon counts is also stored in an ‘EList’ object.” | concepts[68] |
| concept | EList — not defined in this paper's own words; normalized single-channel and voom data class — “Normalized data are stored in ‘MAList’ or ‘EList’ objects. Normalized two-colour data are converted from red and green intensities, R and G , into M and A -values, which hold the log-ratio and average log-intensity values for each spot. Single channel data are background corrected and log 2 transformed and stored in an ‘EList’ object. For RNA-seq data, the voom transformed matrix of gene/exon counts is also stored in an ‘EList’ object.” | concepts[69] |
| concept | MArrayLM — not defined in this paper's own words; fitted linear-model result class — “The next major classes store output from a DE analysis. ‘MArrayLM’ objects store the result of fitting gene-wise linear models to the normalized intensities or log-ratios. Objects of this class are created by the lmFit and eBayes functions. After running decideTests, an object of class ‘TestResults’ stores the results of testing a set of contrasts equal to zero for each probe/gene.” | concepts[70] |
| concept | TestResults — not defined in this paper's own words; contrast-test result class — “The next major classes store output from a DE analysis. ‘MArrayLM’ objects store the result of fitting gene-wise linear models to the normalized intensities or log-ratios. Objects of this class are created by the lmFit and eBayes functions. After running decideTests, an object of class ‘TestResults’ stores the results of testing a set of contrasts equal to zero for each probe/gene.” | concepts[71] |
| concept | matrix — not defined in this paper's own words; lossy coercion target for limma classes — “Furthermore all of these classes may be coerced to be of class matrix using as.matrix, although this entails loss of information. Fitted model objects of class ‘MArrayLM’ can be coerced to class data.frame using as.data.frame in R. The first five classes belong to the virtual class ‘LargeDataObject’ for which a show method is defined to display the leading rows of a large vector, matrix or data.frame.” | concepts[72] |
| concept | C code — not defined in this paper's own words; implementation used for intensive steps — “The limma package is implemented primarily in R ( 65 ) and includes some C code to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genes and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.” | concepts[73] |
| concept | Bioconductor project — not defined in this paper's own words; distribution host for limma — “The limma software is freely available online as part of the Bioconductor project ( http://www.bioconductor.org ). More than 120 other Bioconductor packages make use of limma (as of March 2014). The limma package is used as a building block or as the underlying computational engine by a number of software projects designed to provide user-interfaces for gene expression data analysis including limmaGUI ( 66 ), affylmGUI ( 67 ), WebArray ( 68 ), RACE ( 69 ), CarmaWEB ( 70 ), Goulphar ( 71 ), MAGMA ( 72 ), Asterias ( 73 ), GenePattern ( 74 ), GEO2R ( http://www.ncbi.nlm.nih.gov/geo/geo2r ), the EBI expression atlas ( 75 ), Guide ( 76 ) and Degust ( http://www.vicbioinformatics.com/degust ).” | concepts[74] |
| concept | voom methodology — not defined in this paper's own words; extension unlocking RNA-seq analysis — “Although originally developed with microarray data in mind, the development of the voom methodology unlocks the majority of analysis methods for use on RNA-seq data, such as random effects modelling and gene set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seq) and experimental design.” | concepts[75] |
| concept | reproducible research — not defined in this paper's own words; goal of code-data report compendia — “Being R-based, reports of limma analyses can be compiled using Sweave ( 77 ) or knitr ( 78 ) and provided along with the raw data in a compendium to promote reproducible research in genomics ( 79 ).” | concepts[76] |
Full text — tagged
limmaterm powers differential expressionterm analyses for RNA-sequencing and microarray studies Abstract limmaterm is an R/Bioconductor software packageterm that provides an integrated solution for analysing data from gene expression experimentsterm. It contains rich features for handling complex experimental designs and for information borrowingterm to overcome the problem of small sample sizes. Over the past decade, limmaterm has been a popular choice for geneterm·var discovery through differential expressionterm analyses of microarray and high-throughput PCR data. The package contains particularly strong facilities for reading, normalizing and exploring such data. Recently, the capabilities of limmaterm have been significantly expanded in two important directions. First, the package can now perform both differential expressionterm and differential splicingterm analyses of RNA sequencingterm (RNA-seqterm) data. All the downstream analysis tools previously restricted to microarray data are now available for RNA-seqterm as well. These capabilities allow users to analyse both RNA-seqterm and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genesterm·var or in terms of higher-order expression signaturesterm. This provides enhanced possibilities for biological interpretation of geneterm·var expression differences. This article reviews the philosophy and design of the limmaterm package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described.
INTRODUCTION Geneterm·var expression technologies are used frequently in molecular biology research to gain a snapshot of transcriptional activity in different tissues or populations of cells. These profiles are then compared to identify geneterm·var expression changes associated with a treatment condition or phenotype of interest. Geneterm·var expression studies may be randomized designed experiments in which a biological system is perturbed, for example by a geneterm·var knock-out or by applying a specified stressor. Such experiments are amongst the most powerful tools in functional genomics, providing insights into normal cellular processes as well as disease pathogenesis. Or they may be observational studies in which different phenotypes are compared, diseased and normal tissue for example or cells from different populations. Such studies are common in cancer research and in the study of cell development. In either case, the study design can range from simple two group comparisons to complex set-ups with several experimental factors varying over multiple levels. Researchers might be interested for example in whether a particular geneterm·var facilitates or blocks the action of a particular drug, in which case knock-down and wild-type samples both with and without drug treatment would be profiled. Observational studies may involve multiple batch effectsterm and covariates that must be accounted for in the analysis.
Despite the complexity, geneterm·var expression studies often involve only a small number of biological replicates. The small but complex nature of geneterm·var expression studies poses challenging statistical problems and motivates the use of a number of specialized statistical techniques in order to get the most out of each data set. We have developed the limmaterm software over the past decade to provide a framework for analysing gene expression experimentsterm from beginning to end in a flexible and statistically rigorous way.
The limmaterm package is a core component of Bioconductorterm, an R-based open-source software development projectterm in statistical genomics ( 1 , 2 ). It has proven a popular choice for the analysis of data from experiments involving microarrays ( 3 , 4 ), high-throughput polymerase chain reaction (PCR) ( 5 ), protein arrays ( 6 ) and other platforms. The package is designed in such a way that, after initial pre-processing and normalizationterm, the same analysis pipeline is used for data from all technologies.
Recently, the capabilities of limmaterm have expanded significantly in two important directions. First, the package can now perform both differential expressionterm (DE) and differential splicingterm analyses of RNA sequencingterm (RNA-seqterm) data ( 7 , 8 ). All the downstream analysis tools previously restricted to microarray data are now available for RNA-seqterm as well. These capabilities allow users to analyse both RNA-seqterm and microarray data with very similar pipelines. Second, the package is now able to go past the traditional gene-wise expression analyses in a variety of ways, analysing expression profiles in terms of co-regulated sets of genesterm·var or in terms of higher-order expression signaturesterm ( 7 ). This provides enhanced possibilities for biological interpretation of geneterm·var expression differences.
This article reviews the philosophy and design of the limmaterm package, summarizing both new and historical features, with an emphasis on recent enhancements and features that have not been previously described. The article outlines limmaterm 's functionality at each of the main steps in a geneterm·var expression analysis, from data import, pre-processing, quality assessment and normalizationterm, through to linear modelling, DE and geneterm·var signature analyses.
MATERIALS AND METHODS Figure 3 shows example diagnostic plots. Panel (A) shows RNA-seqterm data from Pickrell et al . ( 9 ) that has been analysed as described by LawpersName et al . ( 10 ). Panels (B) and (C) display the two-colour microarray quality control data set presented by RitchiepersName et al . ( 11 ). Panel (B) displays background corrected but non-normalized intensities from one typical array. Panel (C) was generated from a subset of 30 of the control arrays after print-tip loess normalizationterm ( 12 ).
Example diagnostic plots produced by limmaterm . ( A ) Plot of variability versus count size for RNA-seqterm data, generated by voomterm·def with plot=TRUE. This plot shows that technical variability decreases with count size. Total variability asymptotes to biological variability as count sizes increases. ( B ) Mean-difference plot produced by the plotMA function for a two-colour microarray. The plot highlights negative (NC), constant (DR) and differentially expressed (D03, D10, U03, U10) spike-in controls. Regular probes are non-highlighted. ( C ) Multidimensional scalingterm (MDS) plot of a set of 30 microarrays, generated by plotMDS. All arrays are biologically identical and the plot reveals strong batch effectsterm. Distances represent leading log2-fold changes between samples.
Figure 4 shows example DE summary plots. Panels (A) and (B) were generated using the two-colour microarray data from GEO series GSE2593 . Intensities were background corrected and normalized as previously described ( 13 ). Panel (A) shows a volcano plot for the comparison of samples with RUNX1 over-expressed versus wild-type samples, while panel (B) shows a Venn diagram of differentially expressed probes for each of the three over-expressed genesterm·var versus wild-type. Probes with false discovery rateterm less than 0.05 were considered to be differentially expressed. Panel (C) uses RNA-seqterm data from GEO series GSE52870 . The data were analysed as described in Figure 5 of Liu et al . ( 7 ).
Example plots displaying results from DE and gene set analysesterm. ( A ) Volcano plot showing fold changes and posterior odds of DE for a particular comparison (RUNX1 over-expression versus wild-type in this case), generated by volcanoplot. Probes with P < 0.00001 are highlighted in red. ( B ) Venn diagram showing overlap in the number of DE genesterm·var for three comparisons from the same study as (A), generated by the vennDiagram function. ( C ) Geneterm·var set enrichment plot produced by barcodeplot. The central bar orders differentially expressed genesterm·var by significance from up to down upon Pax5 restoration in an RNA-seqterm experiment ( 7 ). The vertical bars mark genesterm·var that are induced (red) or repressed (blue) upon the transition from large cycling pre-B cells to small resting pre-B cells during normal B cell development according to the published literature ( 47 ). The plot shows a strong positive concordance between Pax5 restoration and the large to small cell transition. The roast function can be used to assign statistical significance to this correlation.
RESULTS Statistical principles limmaterm integrates a number of statistical principles in a way that is effective for large-scale expression studies. It operates on a matrixterm of expression values, where each row represents a geneterm·var, or some other genomic feature relevant to the current study, and each column corresponds to an RNA sampleterm·var. On the one hand, it fits a linear modelterm to each row of data and takes advantage of the flexibility of such models in various ways, for example to handle complex experimental designs and to test very flexible hypotheses. On the other hand, it leverages the highly parallel nature of genomic data to borrow strength between the gene-wise models, allowing for different levels of variability between genesterm·var and between samples, and making statistical conclusions more reliable when the number of samples is small. All the features of the statistical models can be accessed not just for gene-wise expression analyses but also for higher level analyses of geneterm·var expression signaturesterm. Figure 1 depicts the linear modelterm and highlights the statistical principles employed in a typical limmaterm analysis.
Schematic of the major components that are central to any limmaterm analysis. For each geneterm·var g , we have a vector of gene expression valuesterm·var ( y g ) and a design matrixterm X that relates these values to some coefficients of interest (β g ). The limmaterm package includes statistical methods that (i) facilitate information borrowingterm using empirical Bayes methods to obtain posterior variance estimatorsterm ( ), (ii) incorporate observation weights ( w gj where j refers to sample) to allow for variations in data quality, (iii) allow variance modelling to accommodate technical or biological heterogeneity that may be present and (iv) pre-processing methods such as variance stabilization to reduce noise. These methods all help improve inference at both the geneterm·var and geneterm·var set level in small experiments.
Linear modelsterm analyse complete experiments together The hallmark of the limmaterm approach is the use of linear modelsterm to analyse entire experiments as an integrated whole rather than making piece-meal comparisons between pairs of treatmentsterm. This has the effect of sharing information between samples. Analysing the data as a whole also allows us to model correlations that may exist between samples due to repeated measures or other causes. This kind of analysis would not be feasible were the data partitioned into subsets and analysed as a series of pairwise comparisons.
Linear modelsterm permit very general analyses. Researchers can adjust for the effects of multiple experimental factors or can adjust for batch effectsterm. The linear modelterm might include time course effects or regression splines. The linear modelterm could even include the expression values themselves of one or more genesterm·var as covariates, allowing researchers to test for inter-gene dependencies. Linear modelsterm allow researchers to test very flexible hypotheses, not just simple comparisons between groups but also interaction effects or more complex customized comparisons.
Shared global parametersterm link gene-wise models A separate model is fitted for each geneterm·var, but the gene-wise models can be linked by global parametersterm or global hyper-parametersterm. The use of global parametersterm is a simple means of sharing information between genesterm·var that can be used even for the smallest experiments, because the global parametersterm can be estimated from the entire data set involving all the genesterm·var at once. This strategy allows the gene-wise models to incorporate such things as correlations between duplicate probes for the same geneterm·var, or correlations between related RNA samplesterm·var, or variations in quality between the RNA samplesterm·var.
Empirical Bayes borrows information between genesterm·var The highly parallel nature of gene expression experimentsterm lends itself to a particular class of statistical methodsterm, called parametric empirical Bayesterm, that borrow information between genesterm·var in a dynamic way ( 14 , 15 ). The fact that the same linear modelterm is fitted to each geneterm·var allows us to borrow strength between genesterm·var in order to moderate the residual variancesterm ( 16 ). The estimated variance for each geneterm·var then becomes a compromise between the gene-wise estimator, obtained from the data for that geneterm·var alone, and the global variability across all genesterm·var, estimated by pooling the ensemble of all genesterm·var. This has the effect of increasing the effective degrees of freedom with which the gene-wise variances are estimated. It was an innovation of the limmaterm package to show that exact small-sample inference could be conducted using the empirical Bayes posterior variance estimatorsterm ( 16 ). This approach has proven particularly advantageous in experiments with small sample sizes, ensuring that inference is reliable and stable even when the number of replicates is small.
In recent years, the empirical Bayes proceduresterm of limmaterm have been enhanced in two important ways. First, the global variance estimate can now incorporate a mean-variance trendterm ( 10 , 17 , 18 ). This is important because many geneterm·var expression technologies produce data that are less reliable at lower intensities or abundances. Second, the relative weighting of the gene-wise and global variance estimators no longer needs to be the same for all genesterm·var. This allows a sophisticated robust empirical Bayes procedureterm in which hyper-variable genesterm·var are identified and treated separately ( 18 , 19 ). Both of these enhancements improve statistical power and accuracyterm by improving the modelling of the global characteristics of the data in a more flexible way.
Quantitative weightsterm allow for unequal quality Another unique feature of limmaterm is the ability to incorporate quantitative weightsterm into all levels of the statistical analysis, from normalizationterm to linear modelling and geneterm·var set testing. Weights can be applied to genesterm·var or to RNA samplesterm·var or to individual expression values. Weights can be used to give more emphasis to control probes during normalizationterm, or can be used to down-weight measurements or samples that are less reliable in a geneterm·var expression analysis. The weights can be preset based on external quality information, or may be estimated from the expression data itself. The use of weights increases power to detect differentially expressed genesterm, and having a model based approach avoids the need for ad hoc decisions about which observations or samples to filter out ( 11 ).
RNA-seqterm and sequence data All the downstream analysis features of limmaterm are available for RNA-seqterm and other sequence count data, as well as for data from microarrays and other platforms. Traditionally, RNA-seqterm data require specialized software based on the negative binomial or similar distributions ( 20 ). limmaterm however is able to analyse RNA-seqterm read countsterm with high precision by converting counts to the log-scale and estimating the mean-variance relationship empirically (Figure 3A ). The mean-variance trendterm is converted by the voom functionterm into precision weightsterm, which are incorporated into the analysis of log-transformed RNA-seqterm counts using the same linear modelling commands as for microarrays. The resulting pipeline gives comparable performance to the best of the negative binomial-based software packagesterm but with greater speed and reliability for large data sets ( 10 , 21 ). Additionally, and conveniently, only minimal pipeline changes are required when switching between analyses for RNA-seqterm and microarray experiments within limmaterm . This also means that the same statistical tests with the same format of results and graphical displays are available for both data types.
Variance models allow for unequal variability Expression values often show some degree of heteroscedasticityterm, either because there is a relationship between abundance and measurement precision, or because some treatment conditions are more heterogeneous than others. For example, tumours might be more variable than normal tissue. Concern for such effects has prompted some researchers to filter out low intensity observations or to use Welch's t -test for DE between two groups instead of classical pooled t -tests. The use of weights and the ability to model global parametersterm allow limmaterm to incorporate unequal variances in a number of ways. One way is through estimating a mean-variance trendterm, which can either be incorporated into the empirical Bayes procedureterm as mentioned above or used to generate observation weights ( 10 ). A recent development is the ability to estimate precision weightsterm associated with treatment groups or more generally with any given set of covariates. More generally again, the mean-variance trendterm can be estimated in a treatment-specific way, combining the two types of heteroscedasticityterm mentioned above. These approaches allow limmaterm to model unequal variances even for experiments with a small number of RNA samplesterm·var. Importantly, they accommodate unequal variances without compromising the linear modelling and empirical Bayes framework of the package.
Using sets of genesterm·var to represent higher-level expression signaturesterm In recent years, the linear modelling capabilities of limmaterm have been extended to higher-level expression signatureterm analyses involving co-regulated sets of genesterm·var. The idea is to use a set of genesterm·var, together with their log-fold-changes, to represent the transcriptional signature of a biological process or cell type. One way that this is done is by rotation technology, which permits statistical significance to be tested for sets of genesterm·var for any linear modelterm contrast ( 55 ). A particular feature of rotation teststerm is the ability to incorporate prior information about the direction and strength with which each geneterm·var is expected to contribute to the statistical signature. In this way, limmaterm provides a uniquely flexible means to relate new expression data sets to previous results collated from earlier experiments, taking into account for example the fold-change and direction of change for each geneterm·var in the earlier experiment.
A closely related statistical approach implemented in limmaterm is to fit global covariance models, either to estimate correlations between genesterm·var or to estimate the relatedness between the DE profiles resulting from difference comparisons. These new analyses are described briefly later in this article.
Pre-processing methods preserve information Microarray expression data are measured as intensities, which need to be background corrected and normalized before any statistical analysis can be conducted. limmaterm includes a range of background correctionterm and normalizationterm procedures suitable for different types of DNA microarrays or protein arrays. Notable are the maximum likelihood implementation of the normal-exponential convolution model for background correctionterm ( 22 ) and the implementation of loess curves and normalizationterm using quantitative weightsterm. The guiding principle in the pre-processing steps is to preserve information, avoiding missing values or inflated variances ( 23 ). Normalized intensities are offset from zero before transforming to the log-scale to avoid missing values or large variances. Offsets in a range of moderate values have been shown to achieve an effective compromise between noise and bias ( 24 ).
Mean-difference plots Measuring expression in multiple RNA samplesterm·var produces columns of correlated expression values, which are highly correlated because they are measured on the same set of genesterm·var or genomic features. It has long been established in the biomedical literature that the level of agreement between correlated variables can be usefully examined by plotting differences versus means. Such a plot is called a Bland–Altman plotterm ( 25 ) or a Tukey mean-difference plot ( 26 ). Indeed the concept of DE can be viewed as a measure of disagreement between expression measures for the same genesterm·var in different samples. Mean-difference plots were introduced to the two-colour microarray literature by Dudoit et al . ( 27 ) and to the single-channel literature by Bolstad et al . ( 28 ), who called them MA-plotsterm. limmaterm generalized the concept of an MA-plotterm in two ways. First, the idea was extended to apply to sets of single-channel expression values. In this case, the plot is used to compare each sample to the average of all other samples. A virtual array is constructed by averaging the log-expression value for all the samples other than the sample of interest, and then a mean-difference plot is made between the single array and the virtual array. Second, the idea was extended to apply to the fitted model objects. In this case, the plot compares the log-fold-changes for a chosen contrast versus the average log-expression values of each geneterm·var across all samples. In effect, this plots a coefficient of the linear modelterm versus an overall mean intercept parameter. These ideas were part of the original limmaterm package submitted to Bioconductorterm in 2003.
Parametric modelling versus permutation methods It is worth mentioning what limmaterm does not do, which is permutation or re-sampling-based inferenceterm. Permutation is frequently useful in large-scale studies when the aim is to compare two groups. However permutation has a number of disadvantages that make it unattractive for assessing DE in experiments with complex designs. If permutation is applied only to samples involved in two treatment conditions to be compared, then the typically small number of replicates is a severe limitation that will result in low power to detect differences. If permutation is applied to all the samples in a multi-factor experiment, then the composite null hypothesis being testing is an uninteresting one and the power to reject it may be highly dependent on the existence of DE between treatment conditions that are not of primary interest. In other words, permutation cannot be tuned to test specific null hypotheses of interest in a designed experiment. Even more importantly, permutation assumes that all samples are independent and identically distributed under the null hypothesis, and these assumptions are frequently, usually perhaps, unrealistic. In addition, permutation is potentially misleading when the samples are correlated or of unequal precision. In other words, permutation is unable to accommodate blocking structures or quality weights. In small, complex experiments, the potential compromises involved in modelling expression values using parametric distributionsterm, which can never be perfectly correct, are outweighed by the gains in precision and accuracy by modelling the variance structure more realistically.
Pre-processing RNA-seqterm and other sequencing data Figure 2 provides an overview of the functions available at each stage of a geneterm·var expression analysis. The first step is to import expression data into the R session.
The limmaterm workflow. The diagram shows the main steps in a geneterm·var expression analysis, along with individual functions that might be used and the corresponding classes used to store data or results. Online documentation pages are available both for each individual function and for each major step.
limmaterm accepts RNA-seqterm data in the form of a matrixterm of read countsterm, with rows for genomic features and columns for RNA samplesterm·var. Alternatively it can accept a DGEList object from the edgeR package. The genomic regions are often genesterm·var or exons, but could in principle be any genomic feature of interest. In this article, the regions will usually be called genesterm·var for simplicity of terminology. The read countsterm are processed by the voom functionterm in limmaterm to convert them into log 2 counts per millionterm (logCPMterm·var) with associated precision weightsterm. The logCPMterm·var values can be normalized between samples by the voom functionterm or can be pre-normalized by adding normalizationterm factors within edgeR .
Raw read countsterm are assembled outside limmaterm using tools such as featureCounts ( 29 ), HTSeq-counts ( 30 ) or RSEM ( 31 ). The authors of this article find the Subread ( 32 ) and featureCounts pipeline particularly convenient because it is fast, accurate ( 8 ) and can be run from the R prompt using the Rsubread package. The data input to limmaterm should be counts, rather than popular expression summaries such as reads-per-kilobase-per-million (RPKM), so that limmaterm can estimate the appropriate mean-variance relationship. The voomterm·def output can be converted to RPKM values for convenience of interpretation, by subtracting log-gene-lengths, but this should be done after running voomterm·def rather than before.
After running voomterm·def, downstream analysis for RNA-seqterm data is the same as for any other technology. For example, RNA-seqterm data can be explored using boxplots or mean-difference plots, similarly to single-channel microarray data. More detail about this is given in the following sections.
Preprocessing microarray data Reading or importing data For DNA or protein microarrays, importing expression data often involves reading output files created by an image analysis program. Alternatively, a data frame of expression values may be read from a file or data might be directly imported as an R object.
The main limmaterm function to read image output files is read.maimages. This function directly supports formats written by many different image analysis programs including GenePix, Agilent Feature Extraction, ArrayVision, BlueFuse, ImaGene, QuantArray and SPOT (Table 1 ). It also supports the Stanford Microarray Database format. Output in other formats can be read if the appropriate column names are supplied. Two-colour and single-channel data are both supported. Illumina BeadChips need special treatment: output from Illumina's GenomeStudio can be read by read.ilmn if exported as a text file or by read.idat if in binary format.
Standard microarray data formats handled by limmaterm Probe annotation is read automatically if contained in the image output files, or can be read separately and added to the data object. readGAL supports the GenePix Geneterm·var Array List format. read.maimages includes the ability to generate spot quality weights according to any user-specified rule based on any information found in the image output files.
limmaterm includes many possibilities for using or highlighting different types of control probes. The functions readSpotTypes and controlStatus are provided to conveniently classify probes based on text found in the input files. The status of each probe is automatically carried through to appropriate downstream functions.
The function readTargets is provided to read information about the RNA samplesterm·var or targets . This information typically includes information about the treatment conditions and experimental design.
limmaterm can accept data objects containing expression data from other Bioconductorterm packages. It can accept marrayNorm objects from the marray package, PLMset objects from the affyPLM package, vsn objects from the vsn package or objects of any class inheriting from ExpressionSet. Alternatively, expression data can be supplied as a numeric matrixterm. Expression values can be image intensities or normalized log-expression values.
Background correctionterm When array images are read, it is usual to read both foreground and background intensities for each probe. The background intensities can be used to derive an estimate of the ambient intensity affecting each probe. Removing this non-specific signalterm from the foreground intensity of each probe is called background correctionterm and it is typically the first step in processing microarray images. Simply subtracting background from foreground intensities is too heavy-handed ( 23 ). The limmaterm backgroundCorrect function offers a range of more sophisticated alternatives, most unique to the package. These include a method based on a convolution of normal distributions ( 33 ) and a normal-exponential ( normexp ) convolution ( 23 ) with different options for parameter estimation ( 22 ). The plotFB function plots foreground against background intensities for each array and is useful for choosing an appropriate correction method.
Illumina arrays again benefit from special treatment. The nec function implements normexp background correctionterm for Illumina BeadChips making special use of the control probes that are specific to these arrays ( 24 ).
The propexpr function compares intensities to those of negative control probes to estimate the total proportion of probes on each array that correspond to expressed genesterm·var ( 34 ). This provides an estimate of the size of the transcriptome in each sample and is useful for deciding how many probes to filter from downstream analyses.
Normalizationterm Before meaningful comparisons can be made between treatment conditions in a designed experiment, it is critical that the expression values are normalized so that all the samples are as far as possible on the same measurement scale. The purpose of normalizationterm is to remove systematic effectsterm due to technical differences between the assays unassociated with the biological differences of interest. Different technology platforms introduce different biases and so require different normalizationterm methods. The normalizeWithinArrays function normalizes data from two-colour microarrays by aligning the two channels for each array. A popular method is to remove intensity-dependent dye-biases and spatial artefacts from M -values (log-intensity ratios) using locally weighted regression (loess) ( 35 ). The normalizeBetweenArrays function aligns expression values between samples for one-colour microarrays and other single channel platforms using methods such as quantile normalizationterm or cyclic loess ( 28 ). Both functions provide a range of different normalizationterm methods suitable for different platforms. normalizeBetweenArrays also implements separate channel normalizationterm methods for two-colour arrays ( 36 , 37 ). limmaterm is the only software to allow the use of quantitative weightsterm in loess normalizationterm ( 38 ), giving it the ability to downweigh less reliable probes or to give higher priority to control probes or house-keeping genesterm·var. The latter ability has been exploited for normalizing assays when the proportion of differentially expressed genesterm·var may be high, for example boutique arrays ( 39 ), miRNA arrays ( 40 ), PCR arrays, protein arrays or protein mass spectrometry. Other enhancements include the ability to replace the loess curve with a spline curve that has high robustness breakdown properties, and the ability to apply empirical Bayes moderation to the spline curves for multiple regions within the same array (robust-spline normalizationterm).
The neqc function implements quantile normalizationterm for Illumina BeadChips making special use of the control probes specific to these arrays ( 24 ).
All the between-arrays normalizationterm methods are accessible for RNA-seqterm data from within the voom functionterm. Alternatively, voomterm·def has the ability to respect normalizationterm factors computed outside of limmaterm by methods such as trimmed mean of M -values ( 41 ) or conditional quantile normalizationterm ( 42 ).
Graphical exploration of data quality Diagnostic plots allow the user to visually inspect data from a designed experiment in order to identify potential quality problems, such as degraded samples, or problems that arise due to array handling or sample processing. Such displays may also reveal systematic biases that should be removed prior to downstream analysis. Figure 3 presents examples from three different plotting functions.
Plots for individual arrays include the foreground–background plots mentioned above (plotFG), image plots that can reveal inconsistencies across the array surface ( imageplot) and mean-difference plots that show intensity-dependent trends in the log-ratios of two-colour arrays (plotMA, Figure 3B ). The plotMA function can show similar plots for single channel data. In this case, the mean-difference plot is constructed by comparing the log-expression values for that sample compared with the mean of all other samples. The plotMA function makes it simple to highlight particular subsets of probes or genesterm·var, for example control probes. Control probes are automatically highlighted if they have previously been identified using controlStatus (Figure 3B ).
The distribution of expression values can be compared between samples using box plots or density plots (plotDensities). The latter is particularly useful when considering separate channel analyses of two-colour arrays.
Finding differentially expressed genesterm·var Exploration of sample relationships After the pre-processing steps described above, the next major analysis stage is to identify differentially expressed genesterm·var. It is advisable to begin the DE analysis with a plot that visualizes the relative differences in transcriptional profile between the samples. The plotMDS function uses multi-dimensional scalingterm to plot differences in expression profiles between different samples (Figure 3C ). Distances between samples on the plot represent ‘leading fold changeterm’, which is defined as the root-mean-square average of the log-fold-changesterm for the genesterm·var best distinguishing each pair of samples. This provides a type of unsupervised clustering of the samples. It is useful for examining how different are the profiles produced by different experimental factors and for identifying unexpected patterns, such as batch effectsterm, that should be adjusted for during the linear modelterm analysis. This helps guide the construction of the design matrixterm used for the linear modelling below.
The plotRLDF function provides a supervised plot of the samples that shows whether the expression data can distinguish a set of known groups. The function implements a regularized version of linear discriminant functions.
The removeBatchEffect function can be used to remove systematic variation due to batches or other covariates prior to plotting the data so that the effect of treatments can be better seen.
Linear modelling The core component of the limmaterm package is the ability to fit gene-wise linear modelsterm to geneterm·var expression data in order to assess DE ( 16 ). The basic idea is to estimate log-ratios (for two-channel data) or log-intensities (for single-channel data) between two or more target RNA samplesterm·var simultaneously.
Each analysis begins with a matrixterm of expression levels, with probes/genesterm·var/exons in the rows and different samples (biological/technical replicates) in the columns. The linear modelling is performed in a row-wise fashion, with regression coefficients and standard errors either directly estimating the comparisons of interest or via contrasts. Test-statistics are obtained for geneterm·var ranking that can be further summarized at the geneterm·var set level to perform geneterm·var signature/pathway-level ranking.
The flexibility of the linear modelling approach allows almost any experimental design to be handled. Experiments with two or more groups, factorial and time-course designs, and internal controls such as dye-swaps can all be modelled and summarized using the lmFit function. Where appropriate, nuisance variables such as batch and dye effects can also be modelled. Models can be fit robustly or by least squares. Once a linear modelterm is fitted, the makeContrasts function can be used to form a contrast matrixterm. The fitted model object and contrast matrixterm are used by contrasts.fit to compute log 2 -fold-changesterm and t -statistics for the contrasts of interest. This allows all possible pairwise comparisons between treatments to be made.
The plotSA function provides a useful diagnostic plot of the linear modelterm fit, plotting gene-wise residual standard deviations against average log-expression. This allows mean-variance trendsterm to be readily identified, should they exist.
Quality weights and heteroscedasticityterm limmaterm is the only package that allows variations in quality to be handled in a graduated way via quantitative weightsterm. Both observation-level ( 2 , 10 , 43 ) and sample-specific weights ( 11 ) can be used in an analysis. For microarray data, the arrayWeights function estimates relative array variances, which are converted to weights that can be used in the linear modelterm analysis to down-weight observations from less reliable arrays. Probe and array weights can be easily combined by multiplying them together and, when used appropriately, have been demonstrated to increase power to detect DE ( 11 , 43 ).
For RNA-seqterm data, the voomWithQualityWeights function combines observation-level and sample-specific weights for use in the subsequent linear modelling.
Blocking and random effects limmaterm includes a unique strategy for incorporating the fact that observations or samples may be correlated. The strategy is similar to fitting a random effects model, with the difference that all genesterm·var are constrained to share the same intrablock correlation. The duplicateCorrelation function is used to estimate the consensus correlationterm. The correlation structure is then incorporated into the linear modelterm fit and hence into all tests for DE. Originally the idea was used to estimate the correlation between replicate copies of the same probe on a microarray ( 44 ). The correlation strategy preserves more information than simply averaging the replicate probe copies. More generally, the same idea is also used to model the correlation between related RNA samplesterm·var, for example repeated measures on the same individual or RNA samplesterm·var collected at the same time.
Separate channel analysis of two-colour microarrays Two-colour microarrays are traditionally analysed in terms of log-ratios between the two channels hybridized to each probe. limmaterm also provides the possibility of analysing two-colour microarrays as if they were single channel microarrays with two separate samples hybridized to each physical array. This provides a very powerful type of analysis in which intensities can be directly compared between microarrays. The pairing of the red and green channels from each array is kept track of by estimating the correlation between the two channels hybridized to each probe ( 37 ). This type of separate channel analysis uses the intraspotCorrelation and lmscFit functions.
All the linear modelterm fits, whether using lmFit or lmscFit, produce a fitted model object with the same structure. The same fitted model applies regardless of whether correlations have been estimated, whether robust regression or least squares has been used, or whether quality weights have been included. This consistency allows the same framework for DE to be used for all experimental designs and platform technologies.
Testing for DE An empirical Bayes framework to borrow information between genesterm·var when estimating the variances is implemented in the eBayes function. Gene-wise variances are squeezed towards the common or trended variance, which reduces the number of false positivesterm for genesterm·var with very small variances and improves power to detect DE for genesterm·var with larger variances. limmaterm includes a robustified shrinkage strategyterm that allows for gene-wise shrinkage factors to be estimated ( 18 ). This ensures unusually large variances are not squeezed too heavily, reducing the chance that they will appear statistically significant, while more consistently expressed genesterm·var are squeezed more severely towards the common variance. This robust strategy offers the benefits of shrinkage to the majority of the genesterm·var, whilst negating the effects of outliers.
For each coefficient in the linear modelterm or contrast, empirical Bayes moderated t -statistics and their associated P -values are generally used to assess the significance of the observed expression changes. T -statistics can also be translated into Bayesian log-odds of DE. Moderated F -statistics that combine the t -statistics for all contrasts into an overall test of significance for each geneterm·var can also be used.
When one has a particular cut-off for log-fold-change in mind, the treat function can be used to test whether the log 2 -fold-changeterm is greater than a thresholdterm rather than merely different to zeroterm ( 45 ). This can be effective for prioritizing results that are biologically as well as statistically significant.
limmaterm provides a number of options to adjust tests for multiple testing. Users can control either the family-wise type I error rate or the false discovery rateterm ( 46 ). As well as the usual control for multiple testing across multiple genesterm·var, limmaterm is the only software package to provide methods for error rate control across multiple contrasts and genesterm·var simultaneously. For individual tests, multiple testing can be applied using the topTable function. The decideTests function gives access to the full range of options.
To visualize the results of a DE analysis for single or multiple contrasts, limmaterm provides a number of plotting options. Figure 4 shows three such displays: a volcano plot showing the DE results from a single condition, a Venn diagram showing the number of differentially expressed genesterm·var in multiple experimental conditions and a barcode enrichment plot highlighting a particular geneterm·var signature in a DE analysis ranked by moderated t -statistics.
Another useful plot is produced by plotMA, which plots estimated log-fold-changes against mean log-expression for each geneterm·var. This allows the magnitude of changes to be visualized in the context of overall expression level; see for example Figure 5C of Liu et al . ( 7 ).
Testing for differential splicingterm The linear modelterm framework of limmaterm is extended to test very easily for differential splicingterm events when exon-level expression data are available. The data can be either from an exon microarray or from RNA-seqterm data summarized at the exon level. In either case, the approach is based on fitting linear modelsterm to the exon-level expression data. The approach can relate differential exon usage to continuous as well as categorical predictors or to any contrast in a linear modelterm. The test is conducted by the diffSplice function and results are displayed using plotSplice and topSplice. The plotExons function is also useful for exploring exon expression for individual genesterm·var. This approach is considerably faster than alternative approaches to differential splicingterm, making large-scale surveys of differential exon usage feasible.
Higher-level analyses Linear modelling of geneterm·var expression data provides the ideal platform from which to attack larger functional genomic questions related to gene-wise independence or interaction and the decomposition of geneterm·var signatures into distinct molecular pathwaysterm. This section describes higher-level analyses involving multiple genesterm·var.
Estimating the proportion of true null hypothesesterm First we consider DE from a genomic point of view. The DE analyses described so far identify individual differentially expressed genesterm·var according to an FDR criterion. However, for most studies, there are likely to be false negatives: truly differentially expressed genesterm·var that are not detected as differentially expressed because the study did not have enough statistical power to identify them with confidence.
The propTrueNull functionterm estimates the number of truly differentially expressed genesterm·var that remain to be identified. Mathematically, it estimates the proportion of true null hypothesesterm in a collection of hypothesis tests given a vector of P -values. In a geneterm·var expression study, it estimates the proportion of non-differentially expressed genesterm·var, out of all tested, for any contrast in the linear modelterm. The function implements a number of different methods for estimating the proportion of true nulls, ranging from quick and simple to more computationally demanding. The default is based on averaging local false discovery ratesterm across the P -values ( 19 ). Other methods are the histogram method of ( 48 , 49 ), the convex decreasing density estimate of ( 50 ) and a very simple estimate based on averaging the P -values.
Genuine association of geneterm·var expression profiles Gene expression experimentsterm typically involve a number of different treatment conditions. A question that often arises is this: to what extent do two different treatments produce similar or different expression profiles? One way to address this question is to count the overlap in differentially expressed genesterm·var from the two treatments, as in Figure 4B. This approach however is often too crude. It is very sensitive to the significance cut-off used to identify differentially expressed genesterm·var, has little chance of success in situations where power to detect differentially expressed genesterm is relatively low and is subject to technical biases when both treatments are compared back to the same control samples. The genas functionterm ( 19 ) addresses these problems. It tests whether two different contrasts in a linear modelterm affect the same genesterm·var in similar or different ways, adjusting for biases, without needing to apply a significance cut-off for assessing DE. More technically, it estimates the true biological correlationterm between the log 2 -fold-changesterm of two different contrasts. By biological correlationterm, we mean the correlation that would exist between log-fold-changes if they could be measured perfectly without any statistical error. genas is based on a bivariate generalization of the empirical Bayes model that is used to assess DE in limmaterm . This method is particularly powerful for gaining insight into commonly affected geneterm·var pathways when the changes are small but consistent. For instance, applying genas to a microarray study looking at the relationship between polycomb repressor complex (PRC) 1 and PRC2 facilitated the discovery of the opposing roles of these two complexes ( 51 ). This relationship would have been missed if the analysis had been restricted to the statistically significant genesterm·var from each contrast alone.
Geneterm·var set testing Gene set analysesterm assess the overall significance of a set of co-regulated genesterm·var. Each geneterm·var set is chosen to represent a particular molecular pathwayterm or some other biological process of interest. Geneterm·var sets are defined by geneterm·var annotation external to the current expression study, for example from Geneterm·var Ontology (GO) database ( 52 ) or from previous expression studies. For geneterm·var sets defined by previous studies, the genesterm·var may optionally be annotated with the direction and magnitude of expression changes in the earlier experiment. In this way, a geneterm·var set may contain genesterm·var both positively and negatively associated with the molecular pathwayterm that it represents.
limmaterm contains a range of options for geneterm·var set testing via the goana, geneSetTest, camera, roast and romer functions. The goana function provides a traditional GO overlap analysis but with the added ability to adjust for geneterm·var length or abundance biases in RNA-seqterm DE detection. goana uses generalized hypergeometric tests to test for enrichment of GO terms in the list of differentially expressed genesterm·var ( 53 ). It operates directly on the fitted model object and extracts differentially expressed genesterm·var automatically.
Unlike goana, the other geneterm·var set options do not require a significance threshold to be applied to identify differentially expressed genesterm·var. The simplest approach is implemented in the geneSetTest and wilcoxGST functions, which perform rank-based teststerm ( 13 ). These tests assess whether the specified set of genesterm·var is more highly ranked in an ordered list of all genesterm·var than would be expected by chance. These tests have been found to give an effective ranking of biologically significant pathways ( 54 ), but they implicitly assume that the expression level of each geneterm·var is conditionally independent of other genesterm·var and hence give optimistic P -valuesterm ( 55 ).
More sophisticated competitive tests that take into account dependence between the genesterm·var in the linear modelling framework are implemented in the camera functionterm ( 56 ). Camera computes a variance inflation factorterm from the inter-gene correlationterm·var and uses this to adjust the variance of the summary statistics. This avoids optimistic P -valuesterm in the test results, but also reduces statistical power. Camera has been used successfully in a number of biomedical projects ( 57 ).
The roast and mroast functions implement a self-contained test of whether any worthwhile subset of the specified set of genesterm·var is differentially expressed ( 55 ). These functions use rotation teststerm, a specialized simulation technique for multivariate normal models ( 58 ). Rotation can be viewed as a smoothed version of permutationterm that is suitable for linear modelsterm. Rotation occurs only in the residual space of the linear modelterm in such a way that any coefficients in the linear modelterm other than the contrast being tested are held constant ( 55 , 58 ). Roast can provide good statistical power in small complex gene expression experimentsterm. Of all geneterm·var set tests, roast has the unique ability to take into account directional annotation information about genesterm·var in the set. It is able to accommodate genesterm·var both positively and negatively associated with a specified pathway, as well as the magnitude of change of each geneterm·var. It is therefore especially useful for finding similarities in geneterm·var expression patterns between different expression studies ( 57 , 59 – 61 ) (Figure 4C ). Other potential applications for roast include those where the set might not be made up of genesterm·var, for example exon-level expression analyses to test whether any exon of a given geneterm·var is differentially expressed.
Gene set enrichment analysisterm (GSEAterm) is an approach that correlates a large database of co-regulated geneterm·var sets with respect to a microarray or RNA-seqterm data set ( 62 , 63 ). GSEAterm is a hybrid approachterm: it is competitive in that different sets are pitted against one another, but significance is evaluated by permutation of sample labels. The romer function in limmaterm implements a GSEAterm approach that is based on rotation instead of permutation. Like camera and mroast, it can be used with a battery of geneterm·var sets and with any linear modelterm. The limmaterm authors maintain mouse and human versions of the Molecular Signatures Database collections ( 64 ) in R binary format that can be conveniently used with camera, mroast or romer ( http://bioinf.wehi.edu.au/software/MSigDB ).
The use of geneterm·var sets require that geneterm·var symbols and annotation be matched between different databases and studies. Geneterm·var symbols change over time, so limmaterm includes the alias2Symbol and alias2SymbolTable functions to map geneterm·var symbol aliases to current official geneterm·var symbols. The id2indices function matches geneterm·var identifiers in a collection of geneterm·var sets to those in the expression matrixterm, in a format suitable for input into camera, mroast or romer.
A geneterm·var set test can be visualized using the barcodeplot function (Figure 4C ). Genesterm·var are ranked according to their DE results in the current study, then genesterm·var from the a priori set are highlighted by vertical bars, with a smoother line showing the relative enrichment of the geneterm·var set amongst high and low ranked genesterm·var. Barcodeplot is similar to the set location plot introduced by Subramanian et al . ( 63 ), with a number of enhancements, especially the ability to incorporate genesterm·var with positive and negative prior directions. Barcodeplot can optionally display varying weights for different genesterm·var, for example log-fold-changes from a previous experiment.
User-interface Object-oriented programming A simple but appropriate object-oriented paradigmterm provides users with a consistent analysis interface that is very easy from a user point of view. Many limmaterm functions are generic or operate appropriately on objects of different classes. limmaterm defines a number of classes that have been tailored to handle both microarray and RNA-seqterm data. The philosophy has been to define simple list-based data objects that can be easily explored and manipulated by users, in the same style as familiar, long-standing core functions in R such as lm and glm.
For raw intensity data, the classes ‘RGListterm’ and ‘EListRawterm’ are used to store two-colour and single-channel data, respectively. These objects are often created using the function read.maimages and contain the raw values from the image analysis output files along with probe annotation information.
Normalized data are stored in ‘MAListterm’ or ‘EListterm’ objects. Normalized two-colour data are converted from red and green intensities, R and G , into M and A -values, which hold the log-ratio and average log-intensity values for each spot. Single channel data are background corrected and log 2 transformed and stored in an ‘EListterm’ object. For RNA-seqterm data, the voomterm·def transformed matrixterm of geneterm·var/exon counts is also stored in an ‘EListterm’ object.
The next major classes store output from a DE analysis. ‘MArrayLMterm’ objects store the result of fitting gene-wise linear modelsterm to the normalized intensities or log-ratios. Objects of this class are created by the lmFit and eBayes functions. After running decideTests, an object of class ‘TestResultsterm’ stores the results of testing a set of contrasts equal to zero for each probe/geneterm·var.
All of these data classes obey many analogies with matrices. In the case of ‘RGListterm’, ‘MAListterm’, ‘EListRawterm’ and ‘EListterm’, rows correspond to probes/genesterm·var and columns to different samples. In the case of ‘MarrayLMterm’ and ‘TestResultsterm’ rows correspond to unique probes/genesterm·var and the columns to linear modelterm coefficients or contrasts. The standard R functions summary, dim, length, ncol, nrow, dimnames, rownames, colnames have methods for each of these classes. Objects of any of these classes may also be subsetted. Multiple data objects may be combined by rows to add extra probes, or by columns to add extra arrays.
Furthermore all of these classes may be coerced to be of class matrixterm using as.matrixterm, although this entails loss of information. Fitted model objects of class ‘MArrayLMterm’ can be coerced to class data.frame using as.data.frame in R. The first five classes belong to the virtual class ‘LargeDataObject’ for which a show method is defined to display the leading rows of a large vector, matrixterm or data.frame.
Computational efficiency The limmaterm package is implemented primarily in R ( 65 ) and includes some C codeterm to speed up computationally intensive steps. At every stage, effort has been expended to achieve high numerical reliability and efficiency. The memory requirements are linear in the number of genesterm·var and the number of samples. Most estimation procedures finish in a few seconds on a standard desktop computer and virtually all in less than a minute.
Availability The limmaterm software is freely available online as part of the Bioconductor projectterm ( http://www.bioconductorterm.org ). More than 120 other Bioconductorterm packages make use of limmaterm (as of March 2014). The limmaterm package is used as a building block or as the underlying computational engine by a number of software projects designed to provide user-interfaces for geneterm·var expression data analysis including limmaGUI ( 66 ), affylmGUI ( 67 ), WebArray ( 68 ), RACE ( 69 ), CarmaWEB ( 70 ), Goulphar ( 71 ), MAGMA ( 72 ), Asterias ( 73 ), GenePattern ( 74 ), GEO2R ( http://www.ncbi.nlm.nih.gov/geo/geo2r ), the EBI expression atlas ( 75 ), Guide ( 76 ) and Degust ( http://www.vicbioinformatics.com/degust ).
Documentation The limmaterm package provides three levels of documentation. First, each function has its own documentation page that concisely but completely specifies the input data, options and output format of the function. Similarly, each data class has a documentation page explaining all the required and optional components of objects of that class. Care is taken to adhere to the same standards and style that users will be familiar with from help pages in the base R packages.
Second, a series of more general subject help pages serve to link together functions and classes that are used for related purposes. The subject pages cover the topics of (i) introduction, (ii) classes, (iii) reading data, (iv) background correctionterm, (v) normalizationterm, (vi) linear modelsterm, (vii) individual channel analysis of two-colour data, (viii) hypothesis testing for linear modelsterm, (ix) diagnostics and quality assessment, (x) geneterm·var set tests and (xi) RNA-seqterm.
Third, the package comes with an extensive user's guide of over 120 pages, available from the drop-down menu in Windows or alternatively launched by the limmaUsersGuide command. The user's guide gives detailed advice on how to analyse a variety of common study designs. It also includes 10 fully worked case studies for which full data and code are provided.
Users who need more help or advice are invited to post questions to the Bioconductorterm support site ( https://support.bioconductorterm.org ). Questions are usually answered promptly, either by the authors or by other members of the Bioconductorterm community. The support site archives answers to many common questions, including many queries about experimental design and setting up appropriate design matrices.
DISCUSSION This article has summarized the current features of the widely used, open source limmaterm package for geneterm·var expression analysis. This software provides an integrated data analysis solution, using advanced computational algorithms to deliver reliable performance on large data sets and object-oriented ideas to represent expression data and simplify the user interface. New functionality is continually being added as model refinements and new use cases arise.
Although originally developed with microarray data in mind, the development of the voom methodologyterm unlocks the majority of analysis methods for use on RNA-seqterm data, such as random effects modelling and geneterm·var set testing. As with any data analysis problem, the appropriate combination of methods to use will depend upon the biological question, platform used (microarray/RNA-seqterm) and experimental design.
Being R-based, reports of limmaterm analyses can be compiled using Sweave ( 77 ) or knitr ( 78 ) and provided along with the raw data in a compendium to promote reproducible researchterm in genomics ( 79 ).
Applications of limmaterm 's linear modelling strategy beyond the intended analysis of geneterm·var expression data have been made in a variety of applications, including the analysis of data from Nuclear Magnetic Resonance spectroscopy, PCR (including Nanostring), quantitative proteomics ( 80 ), DNA methylation arrays and comparative ChIP-seq ( 81 ).
As the cost of collecting genome-wide profiles continues to fall, we expect the popularity of this approach to continue to grow, with new applications in the analysis of single cell geneterm·var expression data, CRISPR/Cas9 knock-out screens and methylation analysis ( 82 ).
Acknowledgments We are grateful to our colleagues who have been actively involved in the limmaterm project over the years, including James Wettenhall, Jeremy Silver, Davis McCarthy, Natalie Thorne, Aaron Lun, Alicia Oshlack, Ken Simpson, Yang Liao, Yunshun Chen and Carolyn de Graaf. The basic design of the RGListterm and MAListterm classes for two-colour microarrays was based on similar objects defined by the sma package written by Yee-Hwa (Jean) Yang. The limmaterm package has benefited from many other people, too many to list here, who have made suggestions, reported bugs or contributed code.
FUNDING National Health and Medical Research Council (NHMRC) Project Grant [1050661 to M.E.R., G.K.S.; 1023454 to G.K.S., M.E.R., W.S.]; NHMRC Program Grant [1054618 to G.K.S.]; Victorian State Government Operational Infrastructure Support and Australian Government NHMRC IRIISS. Funding for open access charge: NHMRC Program Grant 1054618.
Conflict of interest statement . None declared.
REFERENCES