Introduction to the LIMMA Package

Description

LIMMA is a library for the analysis of gene expression microarray data, especially the use of linear models for analysing designed experiments and the assessment of differential expression. LIMMA provides the ability to analyse comparisons between many RNA targets simultaneously in arbitrary complicated designed experiments. Empirical Bayesian methods are used to provide stable results even when the number of arrays is small. The linear model and differential expression functions apply to all gene expression technologies, including microarrays, RNA-seq and quantitative PCR.

Details

There are three types of documentation available:

• The list("LIMMA User's Guide") can be reached through the "User Guides and Package Vignettes" links at the top of the LIMMA contents page. The function limmaUsersGuide gives the file location of the User's Guide. list()

• An overview of limma functions grouped by purpose is contained in the numbered chapters at the foot of the LIMMA package index page, of which this page is the first. list()

• The LIMMA contents page gives an alphabetical index of detailed help topics. list()

The function changeLog displays the record of changes to the package.

Seealso

02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth, with contributions from many colleagues

References

Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. http://projecteuclid.org/euclid.aoas/1469199900

Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43, e47. http://nar.oxfordjournals.org/content/43/7/e47

Law, CW, Chen, Y, Shi, W, and Smyth, GK (2014). Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology , Volume 3, Article 3. http://www.statsci.org/smyth/pubs/ebayes.pdf

Topic: Classes Defined by this Package

Description

This package defines the following data classes. list(" ", list(list(" ", list(list("RGList")), " "), list(" ", " A class used to store raw intensities as they are read in from an image analysis output file, ", " usually by ", list("read.maimages"), ".")), " ", " ", list(list(" ", list(list("MAList")), " "), list(" ", " Intensities converted to M-values and A-values, i.e., to with-spot and whole-spot contrasts on the log-scale. ", " Usually created from an ", list("RGList"), " using ", list("MA.RG"), " or ", list("normalizeWithinArrays"),

".

", " Objects of this class contain one row for each spot. ", " There may be more than one spot and therefore more than one row for each probe.")), " ", " ", list(list(" ", list(list("EListRaw")), " "), list(" ", " A class to store raw intensities for one-channel microarray data. ", " May or may not be background corrected. ", " Usually created by ", list("read.maimages"), ".")), " ", " ", list(list(" ", list(list("EList")), " "), list(" ", " A class to store normalized log2 expression values for one-channel microarray data. ",

"   Usually created by ", list("normalizeBetweenArrays"), ".")), "

", " ", list(list(" ", list(list("MArrayLM")), " "), list(" ", " Store the result of fitting gene-wise linear models to the normalized intensities or log-ratios. ", " Usually created by ", list("lmFit"), ". ", " Objects of this class normally contain only one row for each unique probe.")), " ", " ", list(list(" ", list(list("TestResults")), " "), list(" ", " Store the results of testing a set of contrasts equal to zero for each probe. ",

"   Usually created by ", list("decideTests"), ".

", " Objects of this class normally contain one row for each unique probe.")), " ")

All these data classes obey many analogies with matrices. In the case of RGList , MAList , EListRaw and EList , rows correspond to spots or probes and columns to arrays. In the case of MarrayLM , rows correspond to unique probes and the columns to parameters or contrasts. The functions summary , dim , length , ncol , nrow , dimnames , rownames , colnames have methods for these classes. Objects of any of these classes may 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 actually 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 .

The first three classes belong to the virtual class LargeDataObject . A show method is defined for LargeDataOject s which uses the utility function printHead .

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Reading Microarray Data from Files

Description

This help page gives an overview of LIMMA functions used to read data from files.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Background Correction

Description

This page deals with background correction methods provided by the backgroundCorrect , kooperberg or neqc functions. Microarray data is typically background corrected by one of these functions before normalization and other downstream analysis.

backgroundCorrect works on matrices, EListRaw or RGList objects, and calls backgroundCorrect.matrix .

The movingmin method of backgroundCorrect uses utility functions ma3x3.matrix and ma3x3.spottedarray .

The normexp method of backgroundCorrect uses utility functions normexp.fit and normexp.signal .

kooperberg is a Bayesian background correction tool designed specifically for two-color GenePix data. It is computationally intensive and requires several additional columns from the GenePix data files. These can be read in using read.maimages and specifying the other.columns argument.

neqc is for single-color data. It performs normexp background correction and quantile normalization using control probes. It uses utility functions normexp.fit.control and normexp.signal . If robust=TRUE , then normexp.fit.control uses the function huber in the MASS package.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Normalization of Microarray Data

Description

This page gives an overview of the LIMMA functions available to normalize data from single-channel or two-colour microarrays. Smyth and Speed (2003) give an overview of the normalization techniques implemented in the functions for two-colour arrays.

Usually data from spotted microarrays will be normalized using normalizeWithinArrays . A minority of data will also be normalized using normalizeBetweenArrays if diagnostic plots suggest a difference in scale between the arrays.

In rare circumstances, data might be normalized using normalizeForPrintorder before using normalizeWithinArrays .

All the normalization routines take account of spot quality weights which might be set in the data objects. The weights can be temporarily modified using modifyWeights to, for example, remove ratio control spots from the normalization process.

If one is planning analysis of single-channel information from the microarrays rather than analysis of differential expression based on log-ratios, then the data should be normalized using a single channel-normalization technique. Single channel normalization uses further options of the normalizeBetweenArrays function. For more details see the LIMMA User's Guide which includes a section on single-channel normalization.

normalizeWithinArrays uses utility functions MA.RG , loessFit and normalizeRobustSpline .

normalizeBetweenArrays is the main normalization function for one-channel arrays, as well as an optional function for two-colour arrays. normalizeBetweenArrays uses utility functions normalizeMedianAbsValues , normalizeMedianAbsValues , normalizeQuantiles and normalizeCyclicLoess , none of which need to be called directly by users.

neqc is a between array normalization function customized for Illumina BeadChips.

The function normalizeVSN is also provided as a interface to the vsn package. It performs variance stabilizing normalization, an algorithm which includes background correction, within and between normalization together, and therefore doesn't fit into the paradigm of the other methods.

removeBatchEffect can be used to remove a batch effect, associated with hybridization time or some other technical variable, prior to unsupervised analysis.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

References

Smyth, G. K., and Speed, T. P. (2003). Normalization of cDNA microarray data. Methods 31, 265-273. http://www.statsci.org/smyth/pubs/normalize.pdf

Topic: Linear Models for Microarrays

Description

This page gives an overview of the LIMMA functions available to fit linear models and to interpret the results. This page covers models for two color arrays in terms of log-ratios or for single-channel arrays in terms of log-intensities. If you wish to fit models to the individual channel log-intensities from two colour arrays, see 07.SingleChannel .

The core of this package is the fitting of gene-wise linear models to microarray data. The basic idea is to estimate log-ratios between two or more target RNA samples simultaneously. See the LIMMA User's Guide for several case studies.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

References

Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology , 3 , No. 1, Article 3. http://www.statsci.org/smyth/pubs/ebayes.pdf

Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics 21(9), 2067-2075.

Topic: Individual Channel Analysis of Two-Color Microarrays

Description

This page gives an overview of the LIMMA functions fit linear models to two-color microarray data in terms of the log-intensities rather than log-ratios.

The function intraspotCorrelation estimates the intra-spot correlation between the two channels. The regression function lmscFit takes the correlation as an argument and fits linear models to the two-color data in terms of the individual log-intensities. The output of lmscFit is an MArrayLM object just the same as from lmFit , so inference proceeds in the same way as for log-ratios once the linear model is fitted. See 06.LinearModels .

The function targetsA2C converts two-color format target data frames to single channel format, i.e, converts from array-per-line to channel-per-line, to facilitate the formulation of the design matrix.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Hypothesis Testing for Linear Models

Description

LIMMA provides a number of functions for multiple testing across both contrasts and genes. The starting point is an MArrayLM object, called fit say, resulting from fitting a linear model and running eBayes and, optionally, contrasts.fit . See 06.LinearModels or 07.SingleChannel for details.

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Diagnostics and Quality Assessment

Description

This page gives an overview of the LIMMA functions available for microarray quality assessment and diagnostic plots.

This package provides an anova method which is designed for assessing the quality of an array series or of a normalization method. It is not designed to assess differential expression of individual genes. anova uses utility functions bwss and bwss.matrix .

The function arrayWeights estimates the empirical reliability of each array following a linear model fit.

Diagnostic plots can be produced by

list(" ", " ", list(list(" ", list(list("imageplot")), " "), list(" ", "Produces a spatial picture of any spot-specific measure from an array image. If the log-ratios are plotted, then this produces an in-silico representation of the well known false-color TIFF image of an array. ", list(list("imageplot3by2")), " will write imageplots to files, six plots to a page.")), " ", " ", list(list(" ", list(list("plotFB")), " "), list(" ", "Plots foreground versus background log-intensies.")), " ",

"

", list(list(" ", list(list("plotMD")), " "), list(" ", "Mean-difference plots. ", "Very versatile plot. ", "For two color arrays, this plots the M-values vs A-values. ", "For single channel technologies, this plots one column of log-expression values vs the average of the other columns. ", "For fitted model objects, this plots a log-fold-change versus average log-expression. ", list(list("mdplot")), " can also be useful for comparing two one-channel microarrays. ")), " ", " ", list(

    list(" ", list(list("plotMA")), " "), list("

", "MA-plots, essentially the same as mean-difference plots. ", list(list("plotMA3by2")), " will write MA-plots to files, six plots to a page. ")), " ", " ", list(list(" ", list(list("plotWithHighlights")), " "), list(" ", "Scatterplots with highlights. ", "This is the underlying engine for ", list("plotMD"), " and ", list("plotMA"), ". ")), " ", " ", list(list(" ", list(list("plotPrintTipLoess")), " "), list(" ", "Produces a grid of MA-plots, one for each print-tip group on an array, together with the corresponding lowess curve. ",

    "Intended to help visualize print-tip loess normalization.")), "

", " ", list(list(" ", list(list("plotPrintorder")), " "), list(" ", "For an array, produces a scatter plot of log-ratios or log-intensities by print order.")), " ", " ", list(list(" ", list(list("plotDensities")), " "), list(" ", "Individual channel densities for one or more arrays. ", "An essential plot to accompany between array normalization, especially quantile normalization.")), " ", " ", list(list(" ", list(list(

    "plotMDS")), " "), list("

", "Multidimensional scaling plot for a set of arrays. ", "Useful for visualizing the relationship between the set of samples.")), " ", " ", list(list(" ", list(list("plotSA")), " "), list(" ", "Sigma vs A plot. ", "After a linear model is fitted, this checks constancy of the variance with respect to intensity level.")), " ")

plotPrintTipLoess uses utility functions gridr and gridc . plotDensities uses utility function RG.MA .

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Gene Set Tests

Description

This page gives an overview of the LIMMA functions for gene set testing and pathway analysis.

list(" ", list(list(" ", list(list("roast")), " "), list(" ", " Self-contained gene set testing for one set.")), " ", " ", list(list(" ", list(list("mroast")), " "), list(" ", " Self-contained gene set testing for many sets.")), " ", " ", list(list(" ", list(list("fry")), " "), list(" ", " Fast approximation to ", list("mroast"), ", especially useful when heteroscedasticity of genes can be ignored.")), " ", " ", list(list(" ", list(list("camera")), " "), list(" ", " Competitive gene set testing.")),

"

", " ", list(list(" ", list(list("romer")), " and ", list(list("topRomer")), " "), list(" ", " Gene set enrichment analysis.")), " ", " ", list(list(" ", list(list("ids2indices")), " "), list(" ", " Convert gene sets consisting of vectors of gene identifiers into a list of indices suitable for use in the above functions.")), " ", " ", list(list(" ", list(list("alias2Symbol")), " and ", list(list("alias2SymbolTable")), " "), list(" ", " Convert gene symbols or aliases to current official symbols.")),

"

", " ", list(list(" ", list(list("geneSetTest")), " or ", list(list("wilcoxGST")), " "), list(" ", " Simple gene set testing based on gene or probe permutation.")), " ", " ", list(list(" ", list(list("barcodeplot")), " "), list(" ", " Enrichment plot of a gene set.")), " ", " ", list(list(" ", list(list("goana")), " and ", list(list("topGO"))), list(" ", " Gene ontology over-representation analysis of gene lists using Entrez Gene IDs. ", " ", list("goana"), " can work directly on a fitted model object or on one or more lists of genes.")),

"

", " ", list(list(" ", list(list("kegga")), " and ", list(list("topKEGG"))), list(" ", " KEGG pathway over-representation analysis of gene lists using Entrez Gene IDs. ", " ", list("kegga"), " can work directly on a fitted model object or on one or more lists of genes.")), " ")

Seealso

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

Author

Gordon Smyth

Topic: Analysis of RNA-seq Data

Description

This page gives an overview of LIMMA functions to analyze RNA-seq data.

list(" ", list(list(" ", list(list("voom")), " "), list(" ", " Transform RNA-seq or ChIP-seq counts to log counts per million (log-cpm) with associated precision weights. ", " After this tranformation, RNA-seq or ChIP-seq data can be analyzed using the same functions as would be used for microarray data.")), " ", " ", list(list(" ", list(list("voomWithQualityWeights")), " "), list(" ", " Combines the functionality of ", list("voom"), " and ", list("arrayWeights"), ".")), " ", " ", list(

list(" ", list(list("diffSplice")), " "), list("

", " Test for differential exon usage between experimental conditions.")), " ", " ", list(list(" ", list(list("topSplice")), " "), list(" ", " Show a data.frame of top results from ", list("diffSplice"), ".")), " ", " ", list(list(" ", list(list("plotSplice")), " "), list(" ", " Plot results from ", list("diffSplice"), ".")), " ", " ", list(list(" ", list(list("plotExons")), " "), list(" ", " Plot logFC for individual exons for a given gene.")),

"

")

Seealso

See also the edgeR package for normalization and data summaries of RNA-seq data, as well as for alternative differential expression methods based on the negative binomial distribution. voom accepts DGEList objects and normalization factors from edgeR.

01.Introduction , 02.Classes , 03.ReadingData , 04.Background , 05.Normalization , 06.LinearModels , 07.SingleChannel , 08.Tests , 09.Diagnostics , 10.GeneSetTests , 11.RNAseq

References

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43, e47. http://nar.oxfordjournals.org/content/43/7/e47

Expression List - class

Description

A list-based S4 classes for storing expression values (E-values), for example for a set of one-channel microarrays or a set of RNA-seq samples. EListRaw holds expression values on the raw scale. EList holds expression values on the log scale, usually after background correction and normalization.

EListRaw objects are often created by read.maimages , while EList objects are often created by normalizeBetweenArrays or by voom . Alternatively, an EList object can be created directly by new("EList",x) , where x is a list.

Seealso

02.Classes gives an overview of all the classes defined by this package.

ExpressionSet is a more formal class in the Biobase package used for the same purpose.

Author

Gordon Smyth

Large Data Object - class

Description

A virtual class including the data classes RGList , MAList and MArrayLM , all of which typically contain large quantities of numerical data in vector, matrices and data.frames.

Seealso

02.Classes gives an overview of all the classes defined by this package.

Author

Gordon Smyth

Examples

#  see normalizeBetweenArrays

Print Layout - class

Description

A list-based class for storing information about the process used to print spots on a microarray.

PrintLayout objects can be created using getLayout . The printer component of an RGList or MAList object is of this class.

Seealso

02.Classes gives an overview of all the classes defined by this package.

Author

Gordon Smyth

Examples

#  Settings for Swirl and ApoAI example data sets in User's Guide

printer <- list(ngrid.r=4, ngrid.c=4, nspot.r=22, nspot.c=24,
ndups=1, spacing=1, npins=16, start="topleft")

#  Typical settings at the Australian Genome Research Facility

#  Full pin set, duplicates side-by-side on same row
printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20,
ndups=2, spacing=1, npins=48, start="topright")

#  Half pin set, duplicates in top and lower half of slide
printer <- list(ngrid.r=12, ngrid.c=4, nspot.r=20, nspot.c=20,
ndups=2, spacing=9600, npins=24, start="topright")

Matrix of Test Results - class

Description

A matrix-based class for storing the results of simultanous tests. TestResults objects are usually created by decideTests .

Usage

list(list("summary"), list("TestResults"))(object, list())
list(list("labels"), list("TestResults"))(object, list())
list(list("levels"), list("TestResults"))(x)

Arguments

Seealso

02.Classes gives an overview of all the classes defined by this package. 08.Tests gives an overview of multiple testing.

Author

Gordon Smyth

Examples

#  Assume a data object y and a design matrix
fit <- lmFit(y, design)
fit <- eBayes(fit)
results <- decideTests(fit)
summary(results)

Convert Gene Aliases to Official Gene Symbols

Description

Maps gene alias names to official gene symbols.

Usage

alias2Symbol(alias, species = "Hs", expand.symbols = FALSE)
alias2SymbolTable(alias, species = "Hs")
alias2SymbolUsingNCBI(alias, gene.info.file,
                      required.columns = c("GeneID","Symbol","description"))

Arguments

Details

Aliases are mapped via NCBI Entrez Gene identity numbers using Bioconductor organism packages.

alias2Symbol maps a set of aliases to a set of symbols, without necessarily preserving order. The output vector may be longer or shorter than the original vector, because some aliases might not be found and some aliases may map to more than one symbol.

alias2SymbolTable returns of vector of the same length as the vector of aliases. If an alias maps to more than one symbol, then the one with the lowest Entrez ID number is returned. If an alias can't be mapped, then NA is returned.

species can be any character string XX for which an organism package org.XX.eg.db exists and is installed. The only requirement of the organism package is that it contains objects org.XX.egALIAS2EG and org.XX.egSYMBOL linking the aliases and symbols to Entrez Gene Ids. At the time of writing, the following organism packages are available from Bioconductor 3.6: list(list("lll"), list(" ", list(), " Package ", list(), " Species", list(), " ", list(), " org.Ag.eg.db ", list(), " Anopheles", list(), " ", list(), " org.Bt.eg.db ", list(), " Bovine", list(), " ", list(), " org.Ce.eg.db ", list(), " Worm", list(), " ", list(), " org.Cf.eg.db ", list(), " Canine", list(), " ", list(), " org.Dm.eg.db ", list(), " Fly", list(), " ", list(), " org.Dr.eg.db ", list(), " Zebrafish", list(), " ", list(), " org.EcK12.eg.db ", list(), " E coli strain K12",

list(), "

", list(), " org.EcSakai.eg.db ", list(), " E coli strain Sakai", list(), " ", list(), " org.Gg.eg.db ", list(), " Chicken", list(), " ", list(), " org.Hs.eg.db ", list(), " Human", list(), " ", list(), " org.Mm.eg.db ", list(), " Mouse", list(), " ", list(), " org.Mmu.eg.db ", list(), " Rhesus", list(), " ", list(), " org.Pt.eg.db ", list(), " Chimp", list(), " ", list(), " org.Rn.eg.db ", list(), " Rat", list(), " ", list(), " org.Ss.eg.db ", list(), " Pig", list(),

"

", list(), " org.Xl.eg.db ", list(), " Xenopus "))

alias2SymbolUsingNCBI is analogous to alias2SymbolTable but uses a gene-info file from NCBI instead of a Bioconductor organism package. It also gives the option of returning multiple columns from the gene-info file. NCBI gene-info files can be downloaded from ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO . For example, the human file is ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz and the mouse file is ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Mus_musculus.gene_info.gz .

Value

alias2Symbol and alias2SymbolTable produce a character vector of gene symbols. alias2SymbolTable returns a vector of the same length and order as alias , including NA values where no gene symbol was found. alias2Symbol returns an unordered vector that may be longer or shorter than alias .

alias2SymbolUsingNCBI returns a data.frame with rows corresponding to the entries of alias and columns as specified by required.columns .

Seealso

This function is often used to assist gene set testing, see 10.GeneSetTests .

Author

Gordon Smyth and Yifang Hu

Examples

alias2Symbol(c("PUMA","NOXA","BIM"), species="Hs")
alias2Symbol("RS1", expand=TRUE)

ANOVA Table - method

Description

Analysis of variance method for objects of class MAList . Produces an ANOVA table useful for quality assessment by decomposing between and within gene sums of squares for a series of replicate arrays. This method produces a single ANOVA Table rather than one for each gene and is not used to identify differentially expressed genes.

Seealso

MAList-class , bwss.matrix , anova .

An overview of quality assessment and diagnostic functions in LIMMA is given by 09.Diagnostics .

Note

This function does not give valid results in the presence of missing M-values.

Author

Gordon Smyth

Array Quality Weights

Description

Estimate relative quality weights for each array or group in a multi-array experiment.

Usage

arrayWeights(object, design = NULL, weights = NULL,
     var.design = NULL, var.group = NULL, prior.n = 10,
     method = "auto", maxiter = 50, tol = 1e-5, trace = FALSE)

Arguments

Details

The relative reliability of each array is estimated by measuring how well the expression values for that array follow the linear model. Arrays that tend to have larger residuals are assigned lower weights.

The basic method is described by Ritchie et al (2006) and the extension to custom variance models by Liu et al (2015). A weighted linear model is fitted to the expression values for each gene. The variance model is fitted to the squared residuals from the linear model fit and is updated either by full REML scoring iterations ( method="reml" ) or using an efficient gene-by-gene update algorithm ( method="genebygene" ). The final estimates of these array variances are converted to weights. The gene-by-gene algorithm is described by Ritchie et al (2006) while the REML algorithm is an adaption of the algorithm of Smyth (2002).

For stability, the array weights are squeezed slightly towards equality. This is done by adding a prior likelihood corresponding to unit array weights equivalent to prior.n genes. The gene-by-gene algorithm is started from the prior genes while the REML algorithm adds the prior to the log-likelihood derivatives.

By default, arrayWeights chooses between the REML and gene-by-gene algorithms automatically ( method="auto" ). REML is chosen if there are no prior weights or missing values and otherwise the gene-by-gene algorithm is used.

The input object is interpreted as for lmFit and getEAWP . In particular, the arguments design and weights will be extracted from the data object if available and do not normally need to be set explicitly in the call; if any of these are set in the call then they will over-ride the slots or components in the data object .

Value

A numeric vector of array weights, which multiply to 1.

Seealso

arrayWeightsQuick , voomWithQualityWeights

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Matthew Ritchie and Gordon Smyth

References

Liu, R., Holik, A. Z., Su, S., Jansz, N., Chen, K., Leong, H. S., Blewitt, M. E., Asselin-Labat, M.-L., Smyth, G. K., Ritchie, M. E. (2015). Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses. Nucleic Acids Research 43, e97. http://nar.oxfordjournals.org/content/43/15/e97

Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7 , 261. http://www.biomedcentral.com/1471-2105/7/261

Smyth, G. K. (2002). An efficient algorithm for REML in heteroscedastic regression. Journal of Computational and Graphical Statistics 11 , 836-847. http://www.statsci.org/smyth/pubs/remlalgo.pdf

Examples

ngenes <- 1000
narrays <- 6
y <- matrix(rnorm(ngenes*narrays), ngenes, narrays)
var.group <- c(1,1,1,2,2,2)
y[,var.group==1] <- 2*y[,var.group==1]
arrayWeights(y, var.group=var.group)

Array Quality Weights

Description

Estimates relative quality weights for each array in a multi-array experiment with replication.

Usage

arrayWeightsQuick(y, fit)

Arguments

Details

Estimates the relative reliability of each array by measuring how well the expression values for that array follow the linear model.

This is a quick and dirty version of arrayWeights .

Value

Numeric vector of weights of length ncol(fit) .

Seealso

See arrayWeights . An overview of LIMMA functions for reading data is given in 03.ReadingData .

Author

Gordon Smyth

References

Ritchie, M. E., Diyagama, D., Neilson, van Laar, R., J., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7, 261. http://www.biomedcentral.com/1471-2105/7/261

Examples

fit <- lmFit(y, design)
arrayWeightsQuick(y, fit)

asMatrixWeights

Description

Convert probe-weights or array-weights to a matrix of weights.

Usage

asMatrixWeights(weights, dim)

Arguments

Details

This function converts a vector or probe-weights or a vector of array-weights to a matrix of the correct size. Probe-weights are repeated across rows while array-weights are repeated down the columns. If weights has length equal to the number of probes, it is assumed to contain probe-weights. If it has length equal to the number of arrays, it is assumed to contain array-weights. If the number of probes is equal to the number of arrays, then weights is assumed to contain array-weights if it is a row-vector of the correct size, i.e., if it is a matrix with one row.

This function is used internally by the linear model fitting functions in limma.

Value

Numeric matrix of dimension dim .

Seealso

modifyWeights .

An overview of functions in LIMMA used for fitting linear models is given in 06.LinearModels .

Author

Gordon Smyth

Examples

asMatrixWeights(1:3,c(4,3))
asMatrixWeights(1:4,c(4,3))

Turn a Microarray Linear Model Object into a Dataframe

Description

Turn a MArrayLM object into a data.frame .

Usage

list(list("as.data.frame"), list("MArrayLM"))(x, row.names = NULL, optional = FALSE, list())

Arguments

Details

This method combines all the components of x which have a row for each probe on the array into a data.frame .

Value

A data.frame.

Seealso

as.data.frame in the base package.

02.Classes gives an overview of data classes used in LIMMA. 06.LinearModels gives an overview of linear model functions in LIMMA.

Author

Gordon Smyth

Convert marrayNorm Object to an MAList Object

Description

Convert marrayNorm Object to an MAList Object

Usage

as.MAList(object)

Arguments

Details

The marrayNorm class is defined in the marray package. This function converts a normalized two color microarray data object created by the marray package into the corresponding limma data object.

Note that such conversion is not necessary to access the limma linear modelling functions, because lmFit will operate on a marrayNorm data object directly.

Value

Object of class MAList

Seealso

02.Classes gives an overview of all the classes defined by this package.

The marrayNorm class is defined in the marray package.

Author

Gordon Smyth

Turn a Microarray Data Object into a Matrix

Description

Turn a microarray data object into a numeric matrix by extracting the expression values.

Usage

list(list("as.matrix"), list("MAList"))(x,list())

Arguments

Details

These methods extract the matrix of log-ratios, for MAList or marrayNorm objects, or the matrix of expression values for other expression objects such as EList or ExressionSet . For MArrayLM objects, the matrix of fitted coefficients is extracted.

These methods involve loss of information, so the original data object is not recoverable.

Value

A numeric matrix.

Seealso

as.matrix in the base package or exprs in the Biobase package.

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Area Under Receiver Operating Curve

Description

Compute exact area under the ROC for empirical data.

Usage

auROC(truth, stat=NULL)

Arguments

Details

A receiver operating curve (ROC) is a plot of sensitivity (true positive rate) versus 1-specificity (false positive rate) for a statistical test or binary classifier. The area under the ROC is a well accepted measure of test performance. It is equivalent to the probability that a randomly chosen pair of cases is corrected ranked.

Here we consider a test statistic stat , with larger values being more significant, and a vector truth indicating whether the alternative hypothesis is in fact true. truth==TRUE or truth==1 indicates a true discovery and truth=FALSE or truth=0 indicates a false discovery. Correct ranking here means that truth[i] is greater than or equal to truth[j] when stat[i] is greater than stat[j] . The function computes the exact area under the empirical ROC curve defined by truth when ordered by stat .

If stat contains ties, then auROC returns the average area under the ROC for all possible orderings of truth for tied stat values.

The area under the curve is undefined if truth is all TRUE or all FALSE or if truth or stat contain missing values.

Value

Numeric value between 0 and 1 giving area under the curve, 1 being perfect and 0 being the minimum.

Author

Gordon Smyth

Examples

auROC(c(1,1,0,0,0))
truth <- rbinom(30,size=1,prob=0.2)
stat <- rchisq(30,df=2)
auROC(truth,stat)

Average Over Replicate Arrays

Description

Condense a microarray data object so that technical replicate arrays are replaced with (weighted) averages.

Usage

list(list("avearrays"), list("default"))(x, ID=colnames(x), weights=NULL)
list(list("avearrays"), list("MAList"))(x, ID=colnames(x), weights=x$weights)
list(list("avearrays"), list("EList"))(x, ID=colnames(x), weights=x$weights)

Arguments

Details

A new data object is computed in which technical replicate arrays are replaced by their (weighted) averages.

For an MAList object, the components M and A are both averaged in this way, as are weights and any matrices found in object$other .

EList objects are similar, except that the E component is averaged instead of M and A .

If x is of mode "character" , then the replicate values are assumed to be equal and the first is taken as the average.

Value

A data object of the same class as x with a column for each unique value of ID .

Seealso

avereps .

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Examples

x <- matrix(rnorm(8*3),8,3)
colnames(x) <- c("a","a","b")
avearrays(x)

Average Over Duplicate Spots

Description

Condense a microarray data object so that values for within-array replicate spots are replaced with their average.

Usage

list(list("avedups"), list("default"))(x, ndups=2, spacing=1, weights=NULL)
list(list("avedups"), list("MAList"))(x, ndups=x$printer$ndups, spacing=x$printer$spacing, weights=x$weights)
list(list("avedups"), list("EList"))(x, ndups=x$printer$ndups, spacing=x$printer$spacing, weights=x$weights)

Arguments

Details

A new data object is computed in which each probe is represented by the (weighted) average of its duplicate spots. For an MAList object, the components M and A are both averaged in this way. For an EList object, the component E is averaged in this way.

If x is of mode "character" , then the duplicate values are assumed to be equal and the first is taken as the average.

Value

A data object of the same class as x with 1/ndups as many rows.

Seealso

avereps .

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Average Over Irregular Replicate Probes

Description

Condense a microarray data object so that values for within-array replicate probes are replaced with their average.

Usage

list(list("avereps"), list("default"))(x, ID=rownames(x), list())
list(list("avereps"), list("MAList"))(x, ID=NULL, list())
list(list("avereps"), list("EList"))(x, ID=NULL, list())

Arguments

Details

A new data object is computed in which each probe ID is represented by the average of its replicate spots or features.

For an MAList object, the components M and A are both averaged in this way, as are weights and any matrices found in object$other . For an MAList object, ID defaults to MA$genes$ID is that exists, otherwise to rownames(MA$M) .

EList objects are similar, except that the E component is averaged instead of M and A .

If x is of mode "character" , then the replicate values are assumed to be equal and the first is taken as the average.

Value

A data object of the same class as x with a row for each unique value of ID .

Seealso

avedups , avearrays . Also rowsum in the base package.

02.Classes gives an overview of data classes used in LIMMA.

Note

This function should only be applied to normalized log-expression values, and not to raw unlogged expression values. It will generate an error message if applied to RGList or EListRaw objects.

Author

Gordon Smyth

Examples

x <- matrix(rnorm(8*3),8,3)
colnames(x) <- c("S1","S2","S3")
rownames(x) <- c("b","a","a","c","c","b","b","b")
avereps(x)

Correct Intensities for Background

Description

Background correct microarray expression intensities.

Usage

backgroundCorrect(RG, method="auto", offset=0, printer=RG$printer,
                  normexp.method="saddle", verbose=TRUE)
backgroundCorrect.matrix(E, Eb=NULL, method="auto", offset=0, printer=NULL,
                  normexp.method="saddle", verbose=TRUE)

Arguments

Details

This function implements the background correction methods reviewed or developed in Ritchie et al (2007) and Silver at al (2009). Ritchie et al (2007) recommend method="normexp" whenever RG contains local background estimates. Silver et al (2009) shows that either normexp.method="mle" or normexp.method="saddle" are excellent options for normexp. If RG contains morphological background estimates instead (available from SPOT or GenePix image analysis software), then method="subtract" performs well.

If method="none" then no correction is done, i.e., the background intensities are treated as zero. If method="subtract" then the background intensities are subtracted from the foreground intensities. This is the traditional background correction method, but is not necessarily recommended. If method="movingmin" then the background estimates are replaced with the minimums of the backgrounds of the spot and its eight neighbors, i.e., the background is replaced by a moving minimum of 3x3 grids of spots.

The remaining methods are all designed to produce positive corrected intensities. If method="half" then any intensity which is less than 0.5 after background subtraction is reset to be equal to 0.5. If method="minimum" then any intensity which is zero or negative after background subtraction is set equal to half the minimum of the positive corrected intensities for that array. If method="edwards" a log-linear interpolation method is used to adjust lower intensities as in Edwards (2003). If method="normexp" a convolution of normal and exponential distributions is fitted to the foreground intensities using the background intensities as a covariate, and the expected signal given the observed foreground becomes the corrected intensity. This results in a smooth monotonic transformation of the background subtracted intensities such that all the corrected intensities are positive.

The normexp method is available in a number of variants depending on how the model parameters are estimated, and these are selected by normexp.method . Here "saddle" gives the saddle-point approximation to maximum likelihood from Ritchie et al (2007) and improved by Silver et al (2009), "mle" gives exact maximum likelihood from Silver at al (2009), "rma" gives the background correction algorithm from the RMA-algorithm for Affymetrix microarray data as implemented in the affy package, and "rma75" gives the RMA-75 method from McGee and Chen (2006). In practice "mle" performs well and is nearly as fast as "saddle" , but "saddle" is the default for backward compatibility. See normexp.fit for more details.

The offset can be used to add a constant to the intensities before log-transforming, so that the log-ratios are shrunk towards zero at the lower intensities. This may eliminate or reverse the usual 'fanning' of log-ratios at low intensities associated with local background subtraction.

Background correction (background subtraction) is also performed by the normalizeWithinArrays method for RGList objects, so it is not necessary to call backgroundCorrect directly unless one wants to use a method other than simple subtraction. Calling backgroundCorrect before normalizeWithinArrays will over-ride the default background correction.

Value

A matrix, EListRaw or RGList object in which foreground intensities have been background corrected and any components containing background intensities have been removed.

Seealso

kooperberg , neqc .

An overview of background correction functions is given in 04.Background .

Author

Gordon Smyth

References

Edwards, D. E. (2003). Non-linear normalization and background correction in one-channel cDNA microarray studies Bioinformatics 19, 825-833.

McGee, M., and Chen, Z. (2006). Parameter estimation for the exponential-normal convolution model for background correction of Affymetrix GeneChip data. Stat Appl Genet Mol Biol , Volume 5, Article 24.

Ritchie, M. E., Silver, J., Oshlack, A., Silver, J., Holmes, M., Diyagama, D., Holloway, A., and Smyth, G. K. (2007). A comparison of background correction methods for two-colour microarrays. Bioinformatics 23, 2700-2707. http://bioinformatics.oxfordjournals.org/content/23/20/2700

Silver, J., Ritchie, M. E., and Smyth, G. K. (2009). Microarray background correction: maximum likelihood estimation for the normal-exponential convolution model. Biostatistics 10, 352-363. http://biostatistics.oxfordjournals.org/content/10/2/352

Examples

RG <- new("RGList", list(R=c(1,2,3,4),G=c(1,2,3,4),Rb=c(2,2,2,2),Gb=c(2,2,2,2)))
backgroundCorrect(RG)
backgroundCorrect(RG, method="half")
backgroundCorrect(RG, method="minimum")
backgroundCorrect(RG, offset=5)

Barcode Enrichment Plot

Description

Display the enrichment of one or two gene sets in a ranked gene list.

Usage

barcodeplot(statistics, index = NULL, index2 = NULL, gene.weights = NULL,
            weights.label = "Weight", labels = c("Down","Up"),
            quantiles = c(-1,1)*sqrt(2), col.bars = NULL, alpha = 0.4,
            worm = TRUE, span.worm = 0.45, xlab = "Statistic", list())

Arguments

Details

This function plots the positions of one or two gene sets in a ranked list of statistics. If there are two sets, then one is considered to be the positive set and the other the down set. For example, the first set and second sets often correspond to genes that are expected to be up- or down-regulated respectively. The function can optionally display varying weights for different genes, for example log-fold-changes from a previous experiment.

The statistics are ranked left to right from smallest to largest. The ranked statistics are represented by a shaded bar or bed, and the positions of the specified subsets are marked by vertical bars, forming a pattern like a barcode. An enrichment worm optionally shows the relative enrichment of the vertical bars in each part of the plot. The worm is computed by the tricubeMovingAverage function.

Barcode plots are often used in conjunction with gene set tests, and show the enrichment of gene sets amongst high or low ranked genes. They were inspired by the set location plot of Subramanian et al (2005), with a number of enhancements, especially the ability to plot positive and negative sets simultaneously. Barcode plots first appeared in the literature in Lim et al (2009). More recent examples can be seen in Liu et al (2014), Sheikh et al (2015), Witkowski et al (2015) and Ng et al (2015).

The function can be used with any of four different calling sequences:

• index is specified, but not index2 or gene.weights . Single direction plot.

• index and index2 are specified. Two directional plot.

• index and gene.weights are specified. gene.weights must have same length as statistics[index] . Plot will be two-directional if gene.weights contains positive and negative values.

• gene.weights is specified by not index or index2 . gene.weights must have same length as statistics . Plot will be two-directional if gene.weights contains positive and negative values.

Value

No value is returned but a plot is produced as a side effect.

Seealso

tricubeMovingAverage , roast , camera , romer , geneSetTest

There is a topic page on 10.GeneSetTests .

Author

Yifang Hu, Gordon Smyth and Di Wu

References

Ng, AP, Hu, Y, Metcalf, D, Hyland, CD, Ierino, H, Phipson, B, Wu, D, Baldwin, TM, Kauppi, M, Kiu, H, Di, Rago, L, Hilton, DJ, Smyth, GK, Alexander, WS (2015). Early lineage priming by trisomy of Erg leads to myeloproliferation in a down syndrome model. PLOS Genetics 11, e1005211. http://www.ncbi.nlm.nih.gov/pubmed/25973911

Lim E, Vaillant F, Wu D, Forrest NC, Pal B, Hart AH, Asselin-Labat ML, Gyorki DE, Ward T, Partanen A, Feleppa F, Huschtscha LI, Thorne HJ; kConFab; Fox SB, Yan M, French JD, Brown MA, Smyth GK, Visvader JE, and Lindeman GJ (2009). Aberrant luminal progenitors as the candidate target population for basal tumor development in BRCA1 mutation carriers. Nature Medicine 15, 907-913.

Liu, GJ, Cimmino, L, Jude, JG, Hu, Y, Witkowski, MT, McKenzie, MD, Kartal-Kaess, M, Best, SA, Tuohey, L, Liao, Y, Shi, W, Mullighan, CG, Farrar, MA, Nutt, SL, Smyth, GK, Zuber, J, and Dickins, RA (2014). Pax5 loss imposes a reversible differentiation block in B progenitor acute lymphoblastic leukemia. Genes & Development 28, 1337-1350. http://www.ncbi.nlm.nih.gov/pubmed/24939936

Sheikh, B, Lee, S, El-saafin, F, Vanyai, H, Hu, Y, Pang, SHM, Grabow, S, Strasser, A, Nutt, SL, Alexander, WS, Smyth, GK, Voss, AK, and Thomas, T (2015). MOZ regulates B cell progenitors in mice, consequently, Moz haploinsufficiency dramatically retards MYC-induced lymphoma development. Blood 125, 1910-1921. http://www.ncbi.nlm.nih.gov/pubmed/25605372

Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, and Mesirov JP (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA 102, 15545-15550.

Witkowski, MT, Cimmino, L, Hu, Y, Trimarchi, T, Tagoh, H, McKenzie, MD, Best, SA, Tuohey, L, Willson, TA, Nutt, SL, Meinrad Busslinger, M, Aifantis, I, Smyth, GK, and Dickins, RA (2015). Activated Notch counteracts Ikaros tumor suppression in mouse and human T cell acute lymphoblastic leukemia. Leukemia 29, 1301-1311. http://www.ncbi.nlm.nih.gov/pubmed/25655195

Examples

stat <- rnorm(100)
sel <- 1:10
sel2 <- 11:20
stat[sel] <- stat[sel]+1
stat[sel2] <- stat[sel2]-1

# One directional
barcodeplot(stat, index = sel)

# Two directional
barcodeplot(stat, index = sel, index2 = sel2)

# Second set can be indicated by negative weights
barcodeplot(stat, index = c(sel,sel2), gene.weights = c(rep(1,10), rep(-1,10)))

# Two directional with unequal weights
w <- rep(0,100)
w[sel] <- runif(10)
w[sel2] <- -runif(10)
barcodeplot(stat, gene.weights = w, weights.label = "logFC")

# One directional with unequal weights
w <- rep(0,100)
w[sel2] <- -runif(10)
barcodeplot(stat, gene.weights = w, weights.label = "logFC", col.bars = "dodgerblue")

Bead Count Weights for Illumina BeadChips

Description

Estimates weights which account for biological variation and technical variation resulting from varying bead numbers.

Usage

beadCountWeights(y, x, design = NULL, bead.stdev = NULL, bead.stderr = NULL,
                 nbeads = NULL, array.cv = TRUE, scale = FALSE)

Arguments

Details

This function estimates optimum weights using the bead statistics for each probe for an Illumina expression BeadChip. It can be used with any Illumina expression BeadChip, but is most likely to be useful with HumanHT-12 BeadChips.

Arguments x and y are both required. x contains the raw expression values and y contains the corresponding log2 values for the same probes and the same arrays after background correction and normalization. x and y be any type of object that can be coerced to a matrix, with rows corresponding to probes and columns to arrays. x and y must contain the same rows and columns in the same order.

The reliability of the normalized expression value for each probe on each array is measured by estimating its technical and biological variability. The bead number weights are the inverse sum of the technical and biological variances.

The technical variance for each probe on each array is inversely proportional to the number of beads and is estimated using array-specific bead-level coefficients of variation.

Coefficients of variation are calculated using raw expression values.

The biological variance for each probe across the arrays are estimated using a Newton iteration, with the assumption that the total residual deviance for each probe from lmFit is inversely proportional to the sum of the technical variance and biological variance.

Only one of bead.stdev or bead.stderr needs to be set. If bead.stdev is not provided, then it will be computed as bead.stderr * sqrt(nbeads) .

If arguments bead.stdev and nbeads are not set explicitly in the call, then they will be extracted from y$other$BEAD_STDEV and y$other$Avg_NBEADS . An EList object containing these components can be created by read.idat or read.ilmn , see the example code below.

Value

A list object with the following components:

*

Seealso

read.ilmn , read.idat , neqc .

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Charity Law and Gordon Smyth

References

Law, CW (2013). Precision weights for gene expression analysis . PhD Thesis. University of Melbourne, Australia. http://repository.unimelb.edu.au/10187/17598

Examples

z <- read.ilmn(files="probesummaryprofile.txt",
ctrfiles="controlprobesummary.txt",
other.columns=c("BEAD_STDEV","Avg_NBEADS"))
y <- neqc(z)
x <- z[z$genes$Status=="regular",]
bcw <- beadCountWeights(y,x,design)
fit <- lmFit(y,design,weights=bcw$weights)
fit <- eBayes(fit)

Block Diagonal Matrix

Description

Form a block diagonal matrix from the given blocks.

Usage

blockDiag(list())

Arguments

Details

This function is sometimes useful for constructing a design matrix for a disconnected two-color microarray experiment in conjunction with modelMatrix .

Value

A block diagonal matrix with dimensions equal to the sum of the input dimensions

Seealso

modelMatrix

Author

Gordon Smyth

Examples

a <- matrix(1,3,2)
b <- matrix(2,2,2)
blockDiag(a,b)

Between and within sums of squares

Description

Sums of squares between and within groups. Allows for missing values.

Usage

bwss(x,group)

Arguments

Details

This is equivalent to one-way analysis of variance.

Value

A list with components

*

Seealso

bwss.matrix

Author

Gordon Smyth

Between and within sums of squares for matrix

Description

Sums of squares between and within the columns of a matrix. Allows for missing values. This function is called by the anova method for MAList objects.

Usage

bwss.matrix(x)

Arguments

Details

This is equivalent to a one-way analysis of variance where the columns of the matrix are the groups. If x is a matrix then bwss.matrix(x) is the same as bwss(x,col(x)) except for speed of execution.

Value

A list with components

*

Seealso

bwss , anova.MAList

Author

Gordon Smyth

Competitive Gene Set Test Accounting for Inter-gene Correlation

Description

Test whether a set of genes is highly ranked relative to other genes in terms of differential expression, accounting for inter-gene correlation.

Usage

list(list("camera"), list("default"))(y, index, design, contrast = ncol(design), weights = NULL,
       use.ranks = FALSE, allow.neg.cor=FALSE, inter.gene.cor=0.01, trend.var = FALSE,
       sort = TRUE, list())
list(list("cameraPR"), list("default"))(statistic, index, use.ranks = FALSE, inter.gene.cor=0.01, sort = TRUE, list())
interGeneCorrelation(y, design)

Arguments

Details

camera and interGeneCorrelation implement methods proposed by Wu and Smyth (2012). camera performs a competitive test in the sense defined by Goeman and Buhlmann (2007). It tests whether the genes in the set are highly ranked in terms of differential expression relative to genes not in the set. It has similar aims to geneSetTest but accounts for inter-gene correlation. See roast for an analogous self-contained gene set test.

The function can be used for any microarray experiment which can be represented by a linear model. The design matrix for the experiment is specified as for the lmFit function, and the contrast of interest is specified as for the contrasts.fit function. This allows users to focus on differential expression for any coefficient or contrast in a linear model by giving the vector of test statistic values.

camera estimates p-values after adjusting the variance of test statistics by an estimated variance inflation factor. The inflation factor depends on estimated genewise correlation and the number of genes in the gene set.

By default, camera uses interGeneCorrelation to estimate the mean pair-wise correlation within each set of genes. camera can alternatively be used with a preset correlation specified by inter.gene.cor that is shared by all sets. This usually works best with a small value, say inter.gene.cor=0.01 .

If interGeneCorrelation=NA , then camera will estimate the inter-gene correlation for each set. In this mode, camera gives rigorous error rate control for all sample sizes and all gene sets. However, in this mode, highly co-regulated gene sets that are biological interpretable may not always be ranked at the top of the list.

With interGeneCorrelation=0.01 , camera will rank biologically interpetable sets more highly. This gives a useful compromise between strict error rate control and interpretable gene set rankings.

cameraPR is a "pre-ranked" version of camera where the genes are pre-ranked according to a pre-computed statistic.

Value

camera and cameraPR return a data.frame with a row for each set and the following columns:

interGeneCorrelation returns a list with components:

*

Seealso

getEAWP

rankSumTestWithCorrelation , geneSetTest , roast , fry , romer , ids2indices .

There is a topic page on 10.GeneSetTests .

Note

The default settings for inter.gene.cor and allow.neg.cor were changed to the current values in limma 3.29.6. Previously, the default was to estimate an inter-gene correlation for each set. To reproduce the earlier default, use allow.neg.cor=TRUE and inter.gene.cor=NA .

Author

Di Wu and Gordon Smyth

References

Wu, D, and Smyth, GK (2012). Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Research 40, e133. http://nar.oxfordjournals.org/content/40/17/e133

Goeman, JJ, and Buhlmann, P (2007). Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23, 980-987.

Examples

y <- matrix(rnorm(1000*6),1000,6)
design <- cbind(Intercept=1,Group=c(0,0,0,1,1,1))

# First set of 20 genes are genuinely differentially expressed
index1 <- 1:20
y[index1,4:6] <- y[index1,4:6]+1

# Second set of 20 genes are not DE
index2 <- 21:40

camera(y, index1, design)
camera(y, index2, design)

camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=NA)
camera(y, list(set1=index1,set2=index2), design, inter.gene.cor=0.01)

# Pre-ranked version
fit <- eBayes(lmFit(y, design))
cameraPR(fit$t[,2], list(set1=index1,set2=index2))

Combine RGList, MAList, EList or EListRaw Objects

Description

Combine a set of RGList , MAList , EList or EListRaw objects.

Usage

list(list("cbind"), list("RGList"))(list(), deparse.level=1)
list(list("rbind"), list("RGList"))(list(), deparse.level=1)

Arguments

Details

cbind combines data objects assuming the same probes in the same order but different arrays. rbind combines data objects assuming equivalent arrays, i.e., the same RNA targets, but different probes.

For cbind , the matrices of expression data from the individual objects are cbinded. The data.frames of target information, if they exist, are rbinded. The combined data object will preserve any additional components or attributes found in the first object to be combined. For rbind , the matrices of expression data are rbinded while the target information, in any, is unchanged.

Value

An RGList , MAList , EList or EListRaw object holding data from all the arrays and all genes from the individual objects.

Seealso

cbind in the base package.

03.ReadingData gives an overview of data input and manipulation functions in LIMMA.

Author

Gordon Smyth

Examples

M <- A <- matrix(11:14,4,2)
rownames(M) <- rownames(A) <- c("a","b","c","d")
colnames(M) <- colnames(A) <- c("A1","A2")
MA1 <- new("MAList",list(M=M,A=A))

M <- A <- matrix(21:24,4,2)
rownames(M) <- rownames(A) <- c("a","b","c","d")
colnames(M) <- colnames(A) <- c("B1","B2")
MA2 <- new("MAList",list(M=M,A=A))

cbind(MA1,MA2)

Limma Change Log

Description

Write as text the most recent changes from the limma package changelog.

Usage

changeLog(n=20)

Arguments

Value

No value is produced, but a number of lines of text are written to standard output.

Seealso

01.Introduction

Author

Gordon Smyth

Examples

changeLog()

Convert Individual Channel Design Matrix to M-A Format

Description

Convert a design matrix in terms of individual channels to ones in terms of M-values or A-values for two-color microarray data.

Usage

designI2M(design)
designI2A(design)

Arguments

Details

If design is a model matrix suitable for modelling individual log-intensities for two color microarray data, then designI2M computes the corresponding model matrix for modelling M-values (log-ratios) and designI2A computes the model matrix for modelling A-values (average log-intensities).

Note that the matrices designI2M(design) or designI2A(design) may be singular if not all of the coefficients are estimable from the M or A-values. In that case there will be columns containing entirely zeros.

Value

numeric model matrix with half as many rows as design

Seealso

model.matrix in the stats package.

An overview of individual channel linear model functions in limma is given by 07.SingleChannel .

Author

Gordon Smyth

Examples

X <- cbind(1,c(1,1,1,1,0,0,0,0),c(0,0,0,0,1,1,1,1))
designI2M(X)
designI2A(X)

Genewise Nested F-Tests

Description

For each gene, classify a series of related t-statistics as significantly up or down using nested F-tests.

Usage

classifyTestsF(object, cor.matrix = NULL, df = Inf, p.value = 0.01, fstat.only = FALSE)

Arguments

Details

classifyTestsF implements the "nestedF" multiple testing option offered by decideTests . Users should generally use decideTests rather than calling classifyTests directly because, by itself, classifyTests does not incorporate any multiple testing adjustment across genes. Instead it simply adjusts for multiple testing across contrasts for each gene individually.

classifyTestsF uses a nested F-test approach giving particular attention to correctly classifying genes that have two or more significant t-statistics, i.e., which are differentially expressed in two or more conditions. For each row of tstat , the overall F-statistics is constructed from the t-statistics as for FStat . At least one constrast will be classified as significant if and only if the overall F-statistic is significant. If the overall F-statistic is significant, then the function makes a best choice as to which t-statistics contributed to this result. The methodology is based on the principle that any t-statistic should be called significant if the F-test is still significant for that row when all the larger t-statistics are set to the same absolute size as the t-statistic in question.

Compared to conventional multiple testing methods, the nested F-test approach achieves better consistency between related contrasts. (For example, if B is judged to be different from C, then at least one of B or C should be different to A.) The approach was first used by Michaud et al (2008). The nested F-test approach provides weak control of the family-wise error rate, i.e., it correctly controls the type I error rate of calling any contrast as significant if all the null hypotheses are true. In other words, it provides error rate control at the overall F-test level but does not provide strict error rate control at the individual contrast level.

Usually object is a limma linear model fitted object, from which a matrix of t-statistics can be extracted, but it can also be a numeric matrix of t-statistics. In either case, rows correspond to genes and columns to coefficients or contrasts. If object is a matrix, then it may be necessary to supply values for cor.matrix and df . The cor.matrix is the same as the correlation matrix of the coefficients from which the t-statistics were calculated and df is the degrees of freedom of the t-statistics. All statistics for the same gene must have the same degrees of freedom.

If fstat.only=TRUE , the classifyTestsF just returns the vector of overall F-statistics for each gene.

Value

If fstat.only=FALSE , then an object of class TestResults is returned. This is essentially a numeric matrix with elements -1 , 0 or 1 depending on whether each t-statistic is classified as significantly negative, not significant or significantly positive respectively.

If fstat.only=TRUE , then a numeric vector of F-statistics is returned with attributes df1 and df2 giving the corresponding degrees of freedom.

Seealso

An overview of multiple testing functions is given in 08.Tests .

Author

Gordon Smyth

References

Michaud, J, Simpson, KM, Escher, R, Buchet-Poyau, K, Beissbarth, T, Carmichael, C, Ritchie, ME, Schutz, F, Cannon, P, Liu, M, Shen, X, Ito, Y, Raskind, WH, Horwitz, MS, Osato, M, Turner, DR, Speed, TP, Kavallaris, M, Smyth, GK, and Scott, HS (2008). Integrative analysis of RUNX1 downstream pathways and target genes. BMC Genomics 9, 363.

Examples

TStat <- matrix(c(0,10,0, 0,5,0, -4,-4,4, 2,2,2), 4, 3, byrow=TRUE)
colnames(TStat) <- paste0("Contrast",1:3)
rownames(TStat) <- paste0("Gene",1:4)
classifyTestsF(TStat, df=20)
FStat <- classifyTestsF(TStat, df=20, fstat.only=TRUE)
P <- pf(FStat, df1=attr(FStat,"df1"), df2=attr(FStat,"df2"), lower.tail=FALSE)
data.frame(F.Statistic=FStat,P.Value=P)

Reform a Design Matrix to that Contrasts Become Coefficients

Description

Reform a design matrix so that one or more coefficients from the new matrix correspond to specified contrasts of coefficients from the old matrix.

Usage

contrastAsCoef(design, contrast=NULL, first=TRUE)

Arguments

Details

If the contrasts contained in the columns of contrast are not linearly dependent, then superfluous columns are dropped until the remaining matrix has full column rank. The number of retained contrasts is stored in qr$rank and the retained columns are given by qr$pivot .

Value

A list with components

*

Seealso

model.matrix in the stats package.

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Gordon Smyth

Examples

design <- cbind(1,c(0,0,1,1,0,0),c(0,0,0,0,1,1))
cont <- c(0,-1,1)
design2 <- contrastAsCoef(design, cont)$design

#  Original coef[3]-coef[2] becomes coef[1]
y <- rnorm(6)
fit1 <- lm(y~0+design)
fit2 <- lm(y~0+design2)
coef(fit1)
coef(fit1)%*% cont
coef(fit2)

Compute Contrasts from Linear Model Fit

Description

Given a linear model fit to microarray data, compute estimated coefficients and standard errors for a given set of contrasts.

Usage

contrasts.fit(fit, contrasts=NULL, coefficients=NULL)

Arguments

Details

This function accepts input from any of the functions lmFit , lm.series , mrlm , gls.series or lmscFit . The function re-orientates the fitted model object from the coefficients of the original design matrix to any set of contrasts of the original coefficients. The coefficients, unscaled standard deviations and correlation matrix are re-calculated in terms of the contrasts.

The idea of this function is to fit a full-rank model using lmFit or equivalent, then use contrasts.fit to obtain coefficients and standard errors for any number of contrasts of the coefficients of the original model. Unlike the design matrix input to lmFit , which normally has one column for each treatment in the experiment, the matrix contrasts may have any number of columns and these are not required to be linearly independent. Methods of assessing differential expression, such as eBayes or classifyTestsF , can then be applied to fitted model object.

The coefficients argument provides a simpler way to specify the contrasts matrix when the desired contrasts are just a subset of the original coefficients.

Value

An list object of the same class as fit , usually MArrayLM . This is a list with components

• Most other components found in fit are passed through unchanged, but t , p.value , lods , F and F.p.value will all be removed.

Seealso

An overview of linear model functions in limma is given by 06.LinearModels .

Note

For efficiency reasons, this function does not re-factorize the design matrix for each probe. A consequence is that, if the design matrix is non-orthogonal and the original fit included precision weights or missing values, then the unscaled standard deviations produced by this function are approximate rather than exact. The approximation is usually acceptable. If not, then the issue can be avoided by redefining the design matrix to fit the contrasts directly.

Even with precision weights, the results from contrasts.fit are always exact if the coefficients being compared are statistically independent. This will always be true, for example, if the original fit was a oneway model and the group-means (no-intercept) parametrization was used for the design matrix.

Author

Gordon Smyth

Examples

#  Simulate gene expression data: 6 microarrays and 100 genes
#  with one gene differentially expressed in first 3 arrays
M <- matrix(rnorm(100*6,sd=0.3),100,6)
M[1,1:3] <- M[1,1:3] + 2
#  Design matrix corresponds to oneway layout, columns are orthogonal
design <- cbind(First3Arrays=c(1,1,1,0,0,0),Last3Arrays=c(0,0,0,1,1,1))
fit <- lmFit(M,design=design)
#  Would like to consider original two estimates plus difference between first 3 and last 3 arrays
contrast.matrix <- cbind(First3=c(1,0),Last3=c(0,1),"Last3-First3"=c(-1,1))
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2)
#  Large values of eb$t indicate differential expression
results <- decideTests(fit2, method="nestedF")
vennCounts(results)

Set Status of each Spot from List of Spot Types

Description

Determine the type (or status) of each spot in the gene list.

Usage

controlStatus(types, genes, spottypecol="SpotType", regexpcol, verbose=TRUE)

Arguments

Details

This function constructs a vector of status codes by searching for patterns in the gene list. The data frame genes contains gene IDs and should have as many rows as there are spots on the microarrays. Such a data frame is often read using readGAL . The data frame types has as many rows as you want to distinguish types of spots in the gene list. This data frame should contain a column or columns, the regexpcol columns, which have the same names as columns in genes and which contain patterns to match in the gene list. Another column, the spottypecol , contains the names of the spot types. Any other columns are assumed to contain plotting parameters, such as colors or symbols, to be associated with the spot types.

The patterns in the regexpcol columns are simplified regular expressions. For example, AA* means any string starting with AA , *AA means any code ending with AA , AA means exactly these two letters, *AA* means any string containing AA , AA. means AA followed by exactly one other character and AA. means exactly AA followed by a period and no other characters. Any other regular expressions are allowed but the codes ^ for beginning of string and $ for end of string should not be included.

Note that the patterns are matched sequentially from first to last, so more general patterns should be included first. For example, it is often a good idea to include a default spot-type as the first line in types with pattern * for all regexpcol columns and default plotting parameters.

Value

Character vector specifying the type (or status) of each spot on the array. Attributes contain plotting parameters associated with each spot type.

Seealso

An overview of LIMMA functions for reading data is given in 03.ReadingData .

Author

Gordon Smyth

Examples

genes <- data.frame(
ID=c("Control","Control","Control","Control","AA1","AA2","AA3","AA4"),
Name=c("Ratio 1","Ratio 2","House keeping 1","House keeping 2",
"Gene 1","Gene 2","Gene 3","Gene 4"))
types <- data.frame(
SpotType=c("Gene","Ratio","Housekeeping"),
ID=c("*","Control","Control"),
Name=c("*","Ratio*","House keeping*"),
col=c("black","red","blue"))
status <- controlStatus(types,genes)

Heatmap of gene expression values

Description

Create a heatmap of a matrix of log-expression values.

Usage

coolmap(x, cluster.by="de pattern", col=NULL,
        linkage.row="complete", linkage.col="complete", show.dendrogram="both", ...)

Arguments

Details

This function calls the heatmap.2 function in the ggplots package with sensible argument settings for genomic log-expression data. The default settings for heatmap.2 are often not ideal for expression data, and overriding the defaults requires explicit calls to hclust and as.dendrogram as well as prior standardization of the data values. The coolmap function implements our preferred defaults for the two most common types of heatmaps. When clustering by relative expression ( cluster.by="de pattern" ), it implements a row standardization that takes account of NA values and standard deviations that might be zero.

Value

A plot is created on the current graphics device. A list is also invisibly returned, see heatmap.2 for details.

Seealso

heatmap.2 , hclust , dist .

An overview of diagnostic functions available in LIMMA is given in 09.Diagnostics .

Author

Gordon Smyth

Examples

# Simulate gene expression data for 50 genes and 6 microarrays.
# Samples are in two groups
# First 50 probes are differentially expressed in second group
ngenes <- 50
sd <- 0.3*sqrt(4/rchisq(ngenes,df=4))
x <- matrix(rnorm(ngenes*6,sd=sd),ngenes,6)
rownames(x) <- paste("Gene",1:ngenes)
x <- x + seq(from=0, to=16, length=ngenes)
x[,4:6] <- x[,4:6] + 2
coolmap(x)

Cumulative Overlap Analysis of Ordered Lists

Description

Test whether the leading members of ordered lists significantly overlap.

Usage

cumOverlap(ol1, ol2)

Arguments

Details

The function compares the top n members of each list, for every possible n , and conducts an hypergeometric test for overlap. The function returns the value of n giving the smallest p-value.

The p-values are adjusted for multiple testing in a similar way to Bonferroni's method, but starting from the top of th e ranked list instead of from the smallest p-values. This approach is designed to be sensitive to contexts where the number of Ids involved in the significant overlap are a small proportion of the total.

The vectors ol1 and ol2 do not need to be of the same length, but only values in common between the two vectors will be used in the calculation.

This method was described in Chapter 4 of Wu (2011).

Value

List containing the following components:

*

Author

Gordon Smyth and Di Wu

References

Wu, D (2011). Finding hidden relationships between gene expression profiles with application to breast cancer biology. PhD thesis, University of Melbourne. http://hdl.handle.net/11343/36278

Examples

GeneIds <- paste0("Gene",1:50)
ol1 <- GeneIds
ol2 <- c(sample(GeneIds[1:5]), sample(GeneIds[6:50]))
coa <- cumOverlap(ol1, ol2)
coa$p.min
coa$id.overlap

Multiple Testing Across Genes and Contrasts

Description

Identify which genes are significantly differentially expressed for each contrast from a fit object containing p-values and test statistics. A number of different multiple testing strategies are offered that adjust for multiple testing down the genes as well as across contrasts for each gene.

Usage

list(list("decideTests"), list("MArrayLM"))(object, method = "separate", adjust.method = "BH", p.value = 0.05,
            lfc = 0, list())
list(list("decideTests"), list("default"))(object, method = "separate", adjust.method = "BH", p.value = 0.05,
            lfc = 0, coefficients = NULL, cor.matrix = NULL, tstat = NULL, df = Inf,
            genewise.p.value = NULL, list())

Arguments

Details

This function can be applied to a matrix of p-values but is more often applied to an MArrayLM fit object produced by eBayes or treat . In either case, rows of object correspond to genes and columns to coefficients or contrasts.

This function applies a multiple testing procedure and a significance level cutoff to the statistics contained in object . It implements a number of multiple testing procedures for determining whether each statistic should be considered significantly different from zero.

method="separate" will apply multiple testing adjustments to each column of p-values separately. Setting method="separate" is equivalent to using topTable separately for each coefficient in the linear model fit and will identify the same probes as significantly differentially expressed if adjust.method is the same. method="global" will treat the entire matrix of t-statistics as a single vector of unrelated tests. method="hierarchical" adjusts down genes and then across contrasts. method="nestedF" adjusts down genes and then uses classifyTestsF to classify contrasts as significant or not for the selected genes.

The default method="separate" and adjust.method="BH" settings are appropriate for most analyses. method="global" is useful when it is important that the same t-statistic cutoff should correspond to statistical significance for all the contrasts. The "nestedF" method was proposed by Michaud et al (2008) and achieves better consistency between contrasts than the other methods. It provides formal error rate control at the gene level but not for individual contrasts. Please see the limma User's Guide for a discussion of the statistical properties of the various methods.

Value

An object of class TestResults . This is essentially a numeric matrix with elements -1 , 0 or 1 depending on whether each t-statistic is classified as significantly negative, not significant or significantly positive.

If lfc>0 then contrasts are judged significant only when the log2-fold change is at least this large in absolute value. For example, one might choose lfc=log2(1.5) to restrict to 50% changes or lfc=1 for 2-fold changes. In this case, contrasts must satisfy both the p-value and the fold-change cutoff to be judged significant.

Seealso

An overview of multiple testing functions is given in 08.Tests .

Note

Although this function enables users to set p-value and lfc cutoffs simultaneously, this combination criterion is not generally recommended. Unless the fold changes and p-values are very highly correlated, the addition of a fold change cutoff can increase the family-wise error rate or false discovery rate above the nominal level. Users wanting to use fold change thresholding are recommended to use treat instead of eBayes and to leave lfc at the default value when using decideTests .

Author

Gordon Smyth

References

Michaud, J, Simpson, KM, Escher, R, Buchet-Poyau, K, Beissbarth, T, Carmichael, C, Ritchie, ME, Schutz, F, Cannon, P, Liu, M, Shen, X, Ito, Y, Raskind, WH, Horwitz, MS, Osato, M, Turner, DR, Speed, TP, Kavallaris, M, Smyth, GK, and Scott, HS (2008). Integrative analysis of RUNX1 downstream pathways and target genes. BMC Genomics 9, 363.

Detection P-Values from Negative Controls

Description

Compute the proportion of negative controls greater than each observed expression value. Particularly useful for Illumina BeadChips.

Usage

list(list("detectionPValues"), list("EListRaw"))(x, status = NULL, list())
list(list("detectionPValues"), list("default"))(x, status, negctrl = "negative", list())

Arguments

Details

The rows of x for which status == negctrl are assumed to correspond to negative control probes.

For each column of x , the detection p-values are defined as (N.eq/2 + N.gt) / N.neg , where N.gt is the number of negative controls with expression greater than the observed value, N.eq is the number of negative controls with expression equal to the observed value, and N.neg is the total number of negative controls.

When used on Illumina BeadChip data, this function produces essentially the same detection p-values as returned by Illumina's GenomeStudio software.

Value

numeric matrix of same dimensions as x containing detection p-values.

Seealso

An overview of LIMMA functions to read expression data is given in 03.ReadingData .

read.idat reads Illumina BeadChip expression data from binary IDAT files.

neqc performs normexp background correction and quantile normalization aided by control probes.

Author

Gordon Smyth

References

Shi, W, de Graaf, C, Kinkel, S, Achtman, A, Baldwin, T, Schofield, L, Scott, H, Hilton, D, Smyth, GK (2010). Estimating the proportion of microarray probes expressed in an RNA sample. Nucleic Acids Research 38(7), 2168-2176. https://www.ncbi.nlm.nih.gov/pubmed/20056656

Examples

# Read Illumina binary IDAT files
x <- read.idat(idat, bgx)
x$genes$DectionPValue <- detectionPValues(x)
y <- neqc(x)

Test for Differential Splicing

Description

Given a linear model fit at the exon level, test for differences in exon retention between experimental conditions.

Usage

diffSplice(fit, geneid, exonid=NULL, robust=FALSE, verbose=TRUE)

Arguments

Details

This function tests for differential exon usage for each gene and for each column of fit .

Testing for differential exon usage is equivalent to testing whether the log-fold-changes in the fit differ between exons for the same gene. Two different tests are provided. The first is an F-test for differences between the log-fold-changes. The other is a series of t-tests in which each exon is compared to the average of all other exons for the same gene. The exon-level t-tests are converted into a genewise test by adjusting the p-values for the same gene by Simes method. The minimum adjusted p-value is then used for each gene.

This function can be used on data from an exon microarray or can be used in conjunction with voom for exon-level RNA-seq counts.

Value

An object of class MArrayLM containing both exon level and gene level tests. Results are sorted by geneid and by exonid within gene.

*

Seealso

topSplice , plotSplice

A summary of functions available in LIMMA for RNA-seq analysis is given in 11.RNAseq .

Author

Gordon Smyth and Charity Law

Examples

v <- voom(dge,design)
fit <- lmFit(v,design)
ex <- diffSplice(fit,geneid="EntrezID")
topSplice(ex)
plotSplice(ex)

Retrieve the Dimensions of an RGList, MAList or MArrayLM Object

Description

Retrieve the number of rows (genes) and columns (arrays) for an RGList, MAList or MArrayLM object.

Usage

list(list("dim"), list("RGList"))(x)

Arguments

Details

Microarray data objects share many analogies with ordinary matrices in which the rows correspond to spots or genes and the columns to arrays. These methods allow one to extract the size of microarray data objects in the same way that one would do for ordinary matrices.

A consequence is that row and column commands nrow(x) , ncol(x) and so on also work.

Value

Numeric vector of length 2. The first element is the number of rows (genes) and the second is the number of columns (arrays).

Seealso

dim in the base package.

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Examples

M <- A <- matrix(11:14,4,2)
rownames(M) <- rownames(A) <- c("a","b","c","d")
colnames(M) <- colnames(A) <- c("A1","A2")
MA <- new("MAList",list(M=M,A=A))
dim(M)
ncol(M)
nrow(M)

Retrieve the Dimension Names of an RGList, MAList, EList, EListRaw or MArrayLM Object

Description

Retrieve the dimension names of a microarray data object.

Usage

list(list("dimnames"), list("RGList"))(x)
list(list("dimnames"), list("RGList"))(x) <- value

Arguments

Details

The dimension names of a microarray object are the same as those of the most important matrix component of that object.

A consequence is that rownames and colnames will work as expected.

Value

Either NULL or a list of length 2. If a list, its components are either NULL or a character vector the length of the appropriate dimension of x .

Seealso

dimnames in the base package.

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Correlation Between Duplicates

Description

Estimate the correlation between duplicate spots (regularly spaced replicate spots on the same array) or between technical replicates from a series of arrays.

Usage

duplicateCorrelation(object, design=NULL, ndups=2, spacing=1, block=NULL,
                     trim=0.15, weights=NULL)

Arguments

Details

When block=NULL , this function estimates the correlation between duplicate spots (regularly spaced within-array replicate spots). If block is not null, this function estimates the correlation between repeated observations on the blocking variable. Typically the blocks are biological replicates and the repeated observations are technical replicates. In either case, the correlation is estimated by fitting a mixed linear model by REML individually for each gene. The function also returns a consensus correlation, which is a robust average of the individual correlations, which can be used as input for functions lmFit or gls.series .

At this time it is not possible to estimate correlations between duplicate spots and between technical replicates simultaneously. If block is not null, then the function will set ndups=1 , which is equivalent to ignoring duplicate spots.

For this function to return statistically useful results, there must be at least two more arrays than the number of coefficients to be estimated, i.e., two more than the column rank of design .

The function may take long time to execute as it fits a mixed linear model for each gene for an iterative algorithm. It is not uncommon for the function to return a small number of warning messages that correlation estimates cannot be computed for some individual genes. This is not a serious concern providing that there are only a few such warnings and the total number of genes is large. The consensus estimator computed by this function will not be materially affected by a small number of genes.

Value

A list with components

*

Seealso

These functions use mixedModel2Fit from the statmod package.

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Gordon Smyth

References

Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics 21(9), 2067-2075. [ http://bioinformatics.oxfordjournals.org/content/21/9/2067 ] [Preprint with corrections: http://www.statsci.org/smyth/pubs/dupcor.pdf ]

Examples

# Simulate gene expression data for 100 probes and 6 microarrays
# Microarray are in two groups
# First two probes are more highly expressed in second group
# Std deviations vary between genes with prior df=4
sd <- 0.3*sqrt(4/rchisq(100,df=4))
y <- matrix(rnorm(100*6,sd=sd),100,6)
rownames(y) <- paste("Gene",1:100)
y[1:2,4:6] <- y[1:2,4:6] + 2
design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
options(digits=3)

# Fit with correlated arrays
# Suppose each pair of arrays is a block
block <- c(1,1,2,2,3,3)
dupcor <- duplicateCorrelation(y,design,block=block)
dupcor$consensus.correlation
fit1 <- lmFit(y,design,block=block,correlation=dupcor$consensus)
fit1 <- eBayes(fit1)
topTable(fit1,coef=2)

# Fit with duplicate probes
# Suppose two side-by-side duplicates of each gene
rownames(y) <- paste("Gene",rep(1:50,each=2))
dupcor <- duplicateCorrelation(y,design,ndups=2)
dupcor$consensus.correlation
fit2 <- lmFit(y,design,ndups=2,correlation=dupcor$consensus)
dim(fit2)
fit2 <- eBayes(fit2)
topTable(fit2,coef=2)

Empirical Bayes Statistics for Differential Expression

Description

Given a microarray linear model fit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value.

Usage

eBayes(fit, proportion = 0.01, stdev.coef.lim = c(0.1,4),
       trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))
treat(fit, lfc = log2(1.2), trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))

Arguments

Details

These functions are used to rank genes in order of evidence for differential expression. They use an empirical Bayes method to squeeze the genewise-wise residual variances towards a common value (or towards a global trend) (Smyth, 2004; Phipson et al, 2016). The degrees of freedom for the individual variances are increased to reflect the extra information gained from the empirical Bayes moderation, resulting in increased statistical power to detect differential expression.

Theese functions accept as input an MArrayLM fitted model object fit produced by lmFit . The columns of fit define a set of contrasts which are to be tested equal to zero. The fitted model object may have been processed by contrasts.fit before being passed to eBayes to convert the coefficients of the original design matrix into an arbitrary number of contrasts.

The empirical Bayes moderated t-statistics test each individual contrast equal to zero. For each gene (row), the moderated F-statistic tests whether all the contrasts are zero. The F-statistic is an overall test computed from the set of t-statistics for that probe. This is exactly analogous the relationship between t-tests and F-statistics in conventional anova, except that the residual mean squares have been moderated between genes.

The estimates s2.prior and df.prior are computed by fitFDist . s2.post is the weighted average of s2.prior and sigma^2 with weights proportional to df.prior and df.residual respectively. The log-odds of differential expression lods was called the B-statistic by Loennstedt and Speed (2002). The F-statistics F are computed by classifyTestsF with fstat.only=TRUE .

eBayes does not compute ordinary t-statistics because they always have worse performance than the moderated versions. The ordinary (unmoderated) t-statistics can, however, can be easily extracted from the linear model output for comparison purposes---see the example code below.

treat computes empirical Bayes moderated-t p-values relative to a minimum meaningful fold-change threshold. Instead of testing for genes that have true log-fold-changes different from zero, it tests whether the true log2-fold-change is greater than lfc in absolute value (McCarthy and Smyth, 2009). In other words, it uses an interval null hypothesis, where the interval is [-lfc,lfc]. When the number of DE genes is large, treat is often useful for giving preference to larger fold-changes and for prioritizing genes that are biologically important. treat is concerned with p-values rather than posterior odds, so it does not compute the B-statistic lods . The idea of thresholding doesn't apply to F-statistics in a straightforward way, so moderated F-statistics are also not computed. When lfc=0 , treat is identical to eBayes , except that F-statistics and B-statistics are not computed. The lfc threshold is usually chosen relatively small, because significantly DE genes must all have fold changes substantially greater than the testing threshold. Typical values for lfc are log2(1.1) , log2(1.2) or log2(1.5) . The top genes chosen by treat can be examined using topTreat .

Note that the lfc testing threshold used by treat to the define the null hypothesis is not the same as a log2-fold-change cutoff, as the observed log2-fold-change needs to substantially larger than lfc for the gene to be called as significant. In practice, modest values for lfc such as log2(1.1) , log2(1.2) or log2(1.5) are usually the most useful. In practice, setting lfc=log2(1.2) or lfc=log2(1.5) will usually cause most differentially expressed genes to have estimated fold-changes of 2-fold or greater, depending on the sample size and precision of the experiment.

The use of eBayes or treat with trend=TRUE is known as the limma-trend method (Law et al, 2014; Phipson et al, 2016). With this option, an intensity-dependent trend is fitted to the prior variances s2.prior . Specifically, squeezeVar is called with the covariate equal to Amean , the average log2-intensity for each gene. The trend that is fitted can be examined by plotSA . limma-trend is useful for processing expression values that show a mean-variance relationship. This is often useful for microarray data, and it can also be applied to RNA-seq counts that have been converted to log2-counts per million (logCPM) values (Law et al, 2014). When applied to RNA-seq logCPM values, limma-trend give similar results to the voom method. The voom method incorporates the mean-variance trend into the precision weights, whereas limma-trend incorporates the trend into the empirical Bayes moderation. limma-trend is somewhat simpler than voom because it assumes that the sequencing depths (library sizes) are not wildly different between the samples and it applies the mean-variance trend on a genewise basis instead to individual observations. limma-trend is recommended for RNA-seq analysis when the library sizes are reasonably consistent (less than 3-fold difference from smallest to largest) because of its simplicity and speed.

If robust=TRUE then the robust empirical Bayes procedure of Phipson et al (2016) is used. This is frequently useful to protect the empirical Bayes procedure against hyper-variable or hypo-variable genes, especially when analysing RNA-seq data. See squeezeVar for more details.

Value

eBayes produces an object of class MArrayLM (see MArrayLM-class ) containing everything found in fit plus the following added components:

• treat a produces an MArrayLM object similar to eBayes but without lods , var.prior , F or F.p.value .

Seealso

squeezeVar , fitFDist , tmixture.matrix , plotSA .

An overview of linear model functions in limma is given by 06.LinearModels .

Note

The algorithm used by eBayes and treat with robust=TRUE was revised slightly in limma 3.27.6. The minimum df.prior returned may be slightly smaller than previously.

Author

Gordon Smyth and Davis McCarthy

References

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

Loennstedt, I., and Speed, T. P. (2002). Replicated microarray data. Statistica Sinica 12 , 31-46.

McCarthy, D. J., and Smyth, G. K. (2009). Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics 25, 765-771. http://bioinformatics.oxfordjournals.org/content/25/6/765

Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. http://projecteuclid.org/euclid.aoas/1469199900

Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, Article 3. http://www.statsci.org/smyth/pubs/ebayes.pdf

Examples

#  See also lmFit examples

#  Simulate gene expression data,
#  6 microarrays and 100 genes with one gene differentially expressed
set.seed(2016)
sigma2 <- 0.05 / rchisq(100, df=10) * 10
y <- matrix(rnorm(100*6,sd=sqrt(sigma2)),100,6)
design <- cbind(Intercept=1,Group=c(0,0,0,1,1,1))
y[1,4:6] <- y[1,4:6] + 1
fit <- lmFit(y,design)

#  Moderated t-statistic
fit <- eBayes(fit)
topTable(fit,coef=2)

#  Ordinary t-statistic
ordinary.t <- fit$coef[,2] / fit$stdev.unscaled[,2] / fit$sigma

# Treat relative to a 10% fold-change
tfit <- treat(fit,lfc=log2(1.1))
topTreat(tfit,coef=2)

Extract Log-Expression Matrix from MAList

Description

Extract the matrix of log-expression values from an MAList object.

Usage

exprs.MA(MA)

Arguments

Details

Converts M and A-values to log-expression values. The output matrix will have two columns for each array, in the order green, red for each array.

This contrasts with as.matrix.MAList which extracts the M-values only, or RG.MA which converts to expression values in RGList form.

Value

A numeric matrix with twice the columns of the input.

Seealso

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Fit Intercept to Vector of Gamma Distributed Variates

Description

Fit Intercept to Vector of Gamma Distributed Variates

Usage

fitGammaIntercept(y,offset=0,maxit=1000)

Arguments

Details

The values y are assumed to follow a gamma distribution with common shape parameter and with expected values given by x+offset . The function implements a globally convergent Newton iteration to estimate x .

Value

Numeric value giving intercept.

Seealso

This function is called by genas .

Author

Gordon Smyth and Belinda Phipson

References

Phipson, B. (2013). Empirical Bayes modelling of expression profiles and their associations . PhD Thesis. University of Melbourne, Australia.

Examples

offset <- runif(10)
x <- 9
mu <- x+offset
y <- rgamma(10,shape=20,scale=mu/20)
fitGammaIntercept(y,offset=offset)

Moment Estimation of Scaled F-Distribution

Description

Moment estimation of the parameters of a scaled F-distribution given one of the degrees of freedom. This function is called internally by eBayes and squeezeVar and is not usually called directly by a user.

Usage

fitFDist(x, df1, covariate=NULL)
fitFDistRobustly(x, df1, covariate=NULL, winsor.tail.p=c(0.05,0.1), trace=FALSE)

Arguments

Details

fitFDist implements an algorithm proposed by Smyth (2004) and Phipson et al (2016). It estimates scale and df2 under the assumption that x is distributed as scale times an F-distributed random variable on df1 and df2 degrees of freedom. The parameters are estimated using the method of moments, specifically from the mean and variance of the x values on the log-scale.

When covariate is supplied, a spline curve trend will be estimated for the x values and the estimation will be adjusted for this trend (Phipson et al, 2016).

fitFDistRobustly is similar to fitFDist except that it computes the moments of the Winsorized values of x , making it robust against left and right outliers. Larger values for winsor.tail.p produce more robustness but less efficiency. When covariate is supplied, a loess trend is estimated for the x values. The robust method is described by Phipson et al (2016).

As well as estimating the F-distribution for the bulk of the cases, i.e., with outliers discounted, fitFDistRobustly also returns an estimated F-distribution with reduced df2 that might be appropriate for each outlier case.

Value

fitFDist produces a list with the following components:

fitFDistRobustly returns the following components as well:

*

Seealso

This function is called by squeezeVar , which in turn is called by eBayes and treat .

This function calls trigammaInverse .

Note

The algorithm used by fitFDistRobustly was revised slightly in limma 3.27.6. The prob.outlier value, which is the lower bound for df2.shrunk , may be slightly smaller than previously.

Author

Gordon Smyth and Belinda Phipson

References

Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology , 3 , No. 1, Article 3. http://www.statsci.org/smyth/pubs/ebayes.pdf

Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. http://projecteuclid.org/euclid.aoas/1469199900

Examples

x <- rf(100,df1=8,df2=16)
fitFDist(x,df1=8)

Fit Mixture Model by Non-Linear Least Squares

Description

Fit Mixture Model by Non-Linear Least Squares

Usage

fitmixture(log2e, mixprop, niter = 4, trace = FALSE)

Arguments

Details

A mixture experiment is one in which two reference RNA sources are mixed in different proportions to create experimental samples. Mixture experiments have been used to evaluate genomic technologies and analysis methods (Holloway et al, 2006). This function uses all the data for each gene to estimate the expression level of the gene in each of two pure samples.

The function fits a nonlinear mixture model to the log2 expression values for each gene. The expected values of log2e for each gene are assumed to be of the form log2( mixprop*Y1 + (1-mixprop)*Y2 ) where Y1 and Y2 are the expression levels of the gene in the two reference samples being mixed. The mixprop values are the same for each gene but Y1 and Y2 are specific to the gene. The function returns the estimated values A=0.5*log2(Y1*Y2) and M=log2(Y2/Y1) for each gene.

The nonlinear estimation algorithm implemented in fitmixture uses a nested Gauss-Newton iteration (Smyth, 1996). It is fully vectorized so that the estimation is done for all genes simultaneously.

Value

List with three components:

*

Author

Gordon K Smyth

References

Holloway, A. J., Oshlack, A., Diyagama, D. S., Bowtell, D. D. L., and Smyth, G. K. (2006). Statistical analysis of an RNA titration series evaluates microarray precision and sensitivity on a whole-array basis. BMC Bioinformatics 7, Article 511. http://www.biomedcentral.com/1471-2105/7/511

Smyth, G. K. (1996). Partitioned algorithms for maximum likelihood and other nonlinear estimation. Statistics and Computing , 6, 201-216. http://www.statsci.org/smyth/pubs/partitio.pdf

Examples

ngenes <- 100
TrueY1 <- rexp(ngenes)
TrueY2 <- rexp(ngenes)
mixprop <- matrix(c(0,0.25,0.75,1),1,4)
TrueExpr <- TrueY1%*% mixprop + TrueY2 %*% (1-mixprop)

log2e <- log2(TrueExpr) + matrix(rnorm(ngenes*4),ngenes,4)*0.1
out <- fitmixture(log2e,mixprop)

# Plot true vs estimated log-ratios
plot(log2(TrueY1/TrueY2), out$M)

Fitted Values Method for MArrayLM Fits

Description

Obtains fitted values from a fitted microarray linear model object.

Usage

list(list("fitted"), list("MArrayLM"))(object, list())

Arguments

Value

A numeric matrix of fitted values.

Seealso

fitted

Author

Gordon Smyth

Genuine Association of Gene Expression Profiles

Description

Calculates biological correlation between two gene expression profiles.

Usage

genas(fit, coef=c(1,2), subset="all", plot=FALSE, alpha=0.4)

Arguments

Details

The function estimates the biological correlation between two different contrasts in a linear model. By biological correlation, we mean the correlation that would exist between the log2-fold changes (logFC) for the two contrasts, if measurement error could be eliminated and the true log-fold-changes were known. This function is motivated by the fact that different contrasts for a linear model are often strongly correlated in a technical sense. For example, the estimated logFC for multiple treatment conditions compared back to the same control group will be positively correlated even in the absence of any biological effect. This function aims to separate the biological from the technical components of the correlation. The method is explained briefly in Majewski et al (2010) and in full detail in Phipson (2013).

The subset argument specifies whether and how the fit object should be subsetted. Ideally, only genes that are truly differentially expressed for one or both of the contrasts should be used estimate the biological correlation. The default is "all" , which uses all genes in the fit object to estimate the biological correlation. The option "Fpval" chooses genes based on how many F-test p-values are estimated to be truly significant using the function propTrueNull . This should capture genes that display any evidence of differential expression in either of the two contrasts. The options "p.union" and "p.int" are based on the moderated t p-values from both contrasts. From the propTrueNull function an estimate of the number of p-values truly significant in either of the two contrasts can be obtained. "p.union" takes the union of these genes and "p.int" takes the intersection of these genes. The other options, "logFC" and "predFC" subsets on genes that attain a logFC or predFC at least as large as the 90th percentile of the log fold changes or predictive log fold changes on the absolute scale.

The plot option is a logical argument that specifies whether or not to plot a scatter plot of log-fold-changes for the two contrasts. The biological and technical correlations are overlaid on the scatterplot using semi-transparent ellipses. library(ellipse) is required to enable the plotting of ellipses.

Value

genas produces a list with the following components:

*

Seealso

lmFit , eBayes , contrasts.fit

Note

As present, genas assumes that technical correlations between coefficients are the same for all genes, and hence it only works with fit objects that were created without observation weights or missing values. It does not work with voom pipelines, because these involve observation weights.

Author

Belinda Phipson and Gordon Smyth

References

Majewski, IJ, Ritchie, ME, Phipson, B, Corbin, J, Pakusch, M, Ebert, A, Busslinger, M, Koseki, H, Hu, Y, Smyth, GK, Alexander, WS, Hilton, DJ, and Blewitt, ME (2010). Opposing roles of polycomb repressive complexes in hematopoietic stem and progenitor cells. Blood 116, 731-739. http://www.bloodjournal.org/content/116/5/731

Phipson, B. (2013). Empirical Bayes modelling of expression profiles and their associations . PhD Thesis. University of Melbourne, Australia. http://repository.unimelb.edu.au/10187/17614

Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43, e47. http://nar.oxfordjournals.org/content/43/7/e47

Examples

#  Simulate gene expression data

#  Three conditions (Control, A and B) and 1000 genes
ngene <- 1000
mu.A <- mu.B <- mu.ctrl <- rep(5,ngene)

#  200 genes are differentially expressed.
#  All are up in condition A and down in B
#  so the biological correlation is negative.
mu.A[1:200] <- mu.ctrl[1:200]+2
mu.B[1:200] <- mu.ctrl[1:200]-2

#  Two microarrays for each condition
mu <- cbind(mu.ctrl,mu.ctrl,mu.A,mu.A,mu.B,mu.B)
y <- matrix(rnorm(6000,mean=mu,sd=1),ngene,6)

# two experimental groups and one control group with two replicates each
group <- factor(c("Ctrl","Ctrl","A","A","B","B"), levels=c("Ctrl","A","B"))
design <- model.matrix(~group)

# fit a linear model
fit <- lmFit(y,design)
fit <- eBayes(fit)

#  Estimate biological correlation between the logFC profiles
#  for A-vs-Ctrl and B-vs-Ctrl
genas(fit, coef=c(2,3), plot=TRUE, subset="F")

Mean-rank Gene Set Test

Description

Test whether a set of genes is highly ranked relative to other genes in terms of a given statistic. Genes are assumed to be independent.

Usage

geneSetTest(index, statistics, alternative = "mixed", type= "auto",
            ranks.only = TRUE, nsim=9999)
wilcoxGST(index, statistics, list())

Arguments

Details

These functions compute a p-value to test the hypothesis that the indexed test set of genes tends to be more highly ranked in terms of some test statistic compared to randomly chosen genes. The statistic might be any statistic of interest, for example a t-statistic or F-statistic for differential expression. Like all gene set tests, these functions can be used to detect differential expression for a group of genes, even when the effects are too small or there is too little data to detect the genes individually.

wilcoxGST is a synonym for geneSetTest with ranks.only=TRUE . This version of the test procedure was developed by Michaud et al (2008), who called it mean-rank gene-set enrichment .

geneSetTest performs a competitive test in the sense that genes in the test set are compared to other genes (Goeman and Buhlmann, 2007). If the statistic is a genewise test statistic for differential expression, then geneSetTest tests whether genes in the set are more differentially expressed than genes not in the set. By contrast, a self-contained gene set test such as roast tests whether genes in the test set are differentially expressed, in an absolute sense, without regard to any other genes on the array.

Because it is based on permuting genes, geneSetTest assumes that the different genes (or probes) are statistically independent. (Strictly speaking, it assumes that the genes in the set are no more correlated on average than randomly chosen genes.) If inter-gene correlations are present, then a statistically significant result from geneSetTest indicates either that the set is highly ranked or that the genes in the set are positively correlated on average (Wu and Smyth, 2012). Unless gene sets with positive correlations are particularly of interest, it may be advisable to use camera or cameraPR instead to adjust the test for inter-gene correlations. Inter-gene correlations are likely to be present in differential expression experiments with biologically heterogeneous experimental units. On the other hand, the assumption of independence between genes should hold when the replicates are purely technical, i.e., when there is no biological variability between the replicate arrays in each experimental condition.

The statistics are usually a set of probe-wise statistics arising for some comparison from a microarray experiment. They may be t-statistics, meaning that the genewise null hypotheses would be rejected for large positive or negative values, or they may be F-statistics, meaning that only large values are significant. Any set of signed statistics, such as log-ratios, M-values or moderated t-statistics, are treated as t-like. Any set of unsigned statistics, such as F-statistics, posterior probabilities or chi-square tests are treated as F-like. If type="auto" then the statistics will be taken to be t-like if they take both positive and negative values and will be taken to be F-like if they are all of the same sign.

There are four possible alternatives to test for. alternative=="up" means the genes in the set tend to be up-regulated, with positive t-statistics. alternative=="down" means the genes in the set tend to be down-regulated, with negative t-statistics. alternative=="either" means the set is either up or down-regulated as a whole. alternative=="mixed" test whether the genes in the set tend to be differentially expressed, without regard for direction. In this case, the test will be significant if the set contains mostly large test statistics, even if some are positive and some are negative.

The latter three alternatives are appropriate when there is a prior expection that all the genes in the set will react in the same direction. The "mixed" alternative is appropriate if you know only that the genes are involved in the relevant pathways, possibly in different directions. The "mixed" is the only meaningful alternative with F-like statistics.

The test statistic used for the gene-set-test is the mean of the statistics in the set. If ranks.only is TRUE the only the ranks of the statistics are used. In this case the p-value is obtained from a Wilcoxon test. If ranks.only is FALSE , then the p-value is obtained by simulation using nsim random sets of genes.

Value

numeric value giving the estimated p-value.

Seealso

cameraPR , camera , roast , barcodeplot , wilcox.test .

There is a topic page on 10.GeneSetTests .

Note

Wu and Smyth (2012) show that geneSetTest does not does correct for inter-gene correlations and is more likely to assign small p-values to sets containing positive correlated genes. The function cameraPR is recommended as a alternative.

Author

Gordon Smyth and Di Wu

References

Wu, D, and Smyth, GK (2012). Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Research , doi: 10.1093/nar/gks461. http://nar.oxfordjournals.org/content/early/2012/05/24/nar.gks461.abstract

Goeman, JJ, and Buhlmann P (2007). Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23, 980-987.

Michaud, J, Simpson, KM, Escher, R, Buchet-Poyau, K, Beissbarth, T, Carmichael, C, Ritchie, ME, Schutz, F, Cannon, P, Liu, M, Shen, X, Ito, Y, Raskind, WH, Horwitz, MS, Osato, M, Turner, DR, Speed, TP, Kavallaris, M, Smyth, GK, and Scott, HS (2008). Integrative analysis of RUNX1 downstream pathways and target genes. BMC Genomics 9, 363. http://www.biomedcentral.com/1471-2164/9/363

Examples

stat <- rnorm(100)
sel <- 1:10; stat[sel] <- stat[sel]+1
wilcoxGST(sel,stat)

Extract Basic Data from Expression Data Objects

Description

Given an expression data object of any known class, get the expression values, weights, probe annotation and A-values that are needed for linear modelling. This function is called by the linear modelling functions in LIMMA.

Usage

getEAWP(object)

Arguments

Details

Rows correspond to probes and columns to RNA samples.

In the case of two-color microarray data objects ( MAList or marrayNorm ), Amean is the vector of row means of the matrix of A-values. For other data objects, Amean is the vector of row means of the matrix of expression values.

From April 2013, the rownames of the output exprs matrix are required to be unique. If object has no row names, then the output rownames of exprs are 1:nrow(object) . If object has row names but with duplicated names, then the rownames of exprs are set to 1:nrow(object) and the original row names are preserved in the ID column of probes .

object should be a normalized data object. getEAWP will return an error if object is a non-normalized data object such as RGList or EListRaw , because these do not contain log-expression values.

Value

A list with components

• exprs is the only required component. The other components will be NULL if not found in the input object.

Seealso

02.Classes gives an overview of data classes used in LIMMA.

Author

Gordon Smyth

Get Numerical Spacing

Description

Convert character to numerical spacing measure for within-array replicate spots.

Usage

getSpacing(spacing, layout)

Arguments

Details

"rows" means that duplicate spots are printed side-by-side by rows. These will be recorded in consecutive rows in the data object.

"columns" means that duplicate spots are printed side-by-sidy by columns. These will be separated in the data object by layout$nspot.r rows.

"subarrays" means that a number of sub-arrays, with identical probes in the same arrangement, are printed on each array. The spacing therefore will be the size of a sub-array.

"topbottom" is the same as "subarrays" when there are two sub-arrays.

Value

Integer giving spacing between replicate spots in the gene list.

Seealso

An overview of LIMMA functions for reading data is given in 03.ReadingData .

Author

Gordon Smyth

Examples

getSpacing("columns",list(ngrid.r=2,ngrid.c=2,nspot.r=20,nspot.c=19))
getSpacing("rows",list(ngrid.r=2,ngrid.c=2,nspot.r=20,nspot.c=19))
getSpacing("topbottom",list(ngrid.r=2,ngrid.c=2,nspot.r=20,nspot.c=19))

Extract the Print Layout of an Array from the GAL File

Description

From the Block, Row and Column information in a genelist, determine the number of grid rows and columns on the array and the number of spot rows and columns within each grid.

Usage

getLayout(gal, guessdups=FALSE)
getLayout2(galfile)
getDupSpacing(ID)

Arguments

Details

A GenePix Array List (GAL) file is a list of genes and associated information produced by an Axon microarray scanner. The function getLayout determines the print layout from a data frame created from a GAL file or gene list. The data.frame must contain columns Block , Column and Row . (The number of tip columns is assumed to be either one or four.)

On some arrays, each probe may be duplicated a number of times ( ndups ) at regular intervals ( spacing ) in the GAL file. getDupSpacing determines valid values for ndups and spacing from a vector of IDs. If guessdups=TRUE , then getLayout calls getDupSpacing .

The function getLayout2 attempts to determine the print layout from the header information of an actual GAL file.

Value

A printlayout object, which is a list with the following components. The last two components are present only if guessdups=TRUE .

*

Seealso

An overview of LIMMA functions for reading data is given in 03.ReadingData .

Author

Gordon Smyth and James Wettenhall

Examples

# gal <- readGAL()
# layout <- getLayout(gal)

Fit Linear Model to Microarray Data by Generalized Least Squares

Description

Fit a linear model genewise to expression data from a series of microarrays. The fit is by generalized least squares allowing for correlation between duplicate spots or related arrays. This is a utility function for lmFit .

Usage

gls.series(M,design=NULL,ndups=2,spacing=1,block=NULL,correlation=NULL,weights=NULL,list())

Arguments

Details

This is a utility function used by the higher level function lmFit . Most users should not use this function directly but should use lmFit instead.

This function is for fitting gene-wise linear models when some of the expression values are correlated. The correlated groups may arise from replicate spots on the same array (duplicate spots) or from a biological or technical replicate grouping of the arrays. This function is normally called by lmFit and is not normally called directly by users.

Note that the correlation is assumed to be constant across genes. If correlation=NULL then a call is made to duplicateCorrelation to estimated the correlation.

Value

A list with components

*

Seealso

duplicateCorrelation .

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Gordon Smyth

Gene Ontology or KEGG Pathway Analysis

Description

Test for over-representation of gene ontology (GO) terms or KEGG pathways in one or more sets of genes, optionally adjusting for abundance or gene length bias.

Usage

list(list("goana"), list("MArrayLM"))(de, coef = ncol(de), geneid = rownames(de), FDR = 0.05, trend = FALSE, list())
list(list("kegga"), list("MArrayLM"))(de, coef = ncol(de), geneid = rownames(de), FDR = 0.05, trend = FALSE, list())
list(list("goana"), list("default"))(de, universe = NULL, species = "Hs", prior.prob = NULL, covariate=NULL,
      plot=FALSE, list())
list(list("kegga"), list("default"))(de, universe = NULL, restrict.universe = FALSE, species = "Hs", species.KEGG = NULL,
      convert = FALSE, gene.pathway = NULL, pathway.names = NULL,
      prior.prob = NULL, covariate=NULL, plot=FALSE, list())
getGeneKEGGLinks(species.KEGG = "hsa", convert = FALSE)
getKEGGPathwayNames(species.KEGG = NULL, remove.qualifier = FALSE)

Arguments

Details

These functions perform over-representation analyses for Gene Ontology terms or KEGG pathways. The default methods accept a gene set as a vector of Entrez Gene IDs or multiple gene sets as a list of such vectors. An over-represention analysis is then done for each set. The MArrayLM method extracts the gene sets automatically from a linear model fit object.

The p-values returned by goana and kegga are unadjusted for multiple testing. The authors have chosen not to correct automatically for multiple testing because GO terms and KEGG pathways are often overlapping, so standard methods of p-value adjustment may be very conservative. Users should be aware though that p-values are unadjusted, meaning that only very small p-values should be used for published results.

goana uses annotation from the appropriate Bioconductor organism package. The species can be any character string XX for which an organism package org.XX.eg.db is installed. Examples are "Hs" for human for "Mm" for mouse. See alias2Symbol for other possible values for species .

kegga reads KEGG pathway annotation from the KEGG website. For kegga , the species name can be provided in either Bioconductor or KEGG format. Examples of KEGG format are "hsa" for human, "mmu" for mouse of "dme" for fly. kegga can be used for any species supported by KEGG, of which there are more than 14,000 possibilities. By default, kegga obtains the KEGG annotation for the specified species from the http://rest.kegg.jp website. Alternatively one can supply the required pathway annotation to kegga in the form of two data.frames. If this is done, then an internet connection is not required.

The gene ID system used by kegga for each species is determined by KEGG. For human and mouse, the default (and only choice) is Entrez Gene ID. For Drosophila, the default is FlyBase CG annotation symbol. The format of the IDs can be seen by typing head(getGeneKEGGLinks(species)) , for example head(getGeneKEGGLinks("hsa")) or head(getGeneKEGGLinks("dme")) . Entrez Gene IDs can always be used. If Entrez Gene IDs are not the default, then conversion can be done by specifying "convert=TRUE" .

Another possibility is to use KEGG orthology IDs as the gene IDs, and these can be used for any species. In that case, set species.KEGG="ko" .

The ability to supply data.frame annotation to kegga means that kegga can in principle be used in conjunction with any user-supplied set of annotation terms.

The default goana and kegga methods accept a vector prior.prob giving the prior probability that each gene in the universe appears in a gene set. This vector can be used to correct for unwanted trends in the differential expression analysis associated with gene length, gene abundance or any other covariate (Young et al, 2010). The MArrayLM object computes the prior.prob vector automatically when trend is non- NULL .

If prior.prob=NULL , the function computes one-sided hypergeometric tests equivalent to Fisher's exact test. If prior probabilities are specified, then a test based on the Wallenius' noncentral hypergeometric distribution is used to adjust for the relative probability that each gene will appear in a gene set, following the approach of Young et al (2010).

The MArrayLM methods performs over-representation analyses for the up and down differentially expressed genes from a linear model analysis. In this case, the universe is all the genes found in the fit object.

trend=FALSE is equivalent to prior.prob=NULL . If trend=TRUE or a covariate is supplied, then a trend is fitted to the differential expression results and this is used to set prior.prob .

The statistical approach provided here is the same as that provided by the goseq package, with one methodological difference and a few restrictions. Unlike the goseq package, the gene identifiers here must be Entrez Gene IDs and the user is assumed to be able to supply gene lengths if necessary. The goseq package has additional functionality to convert gene identifiers and to provide gene lengths. The only methodological difference is that goana and kegga computes gene length or abundance bias using tricubeMovingAverage instead of monotonic regression. While tricubeMovingAverage does not enforce monotonicity, it has the advantage of numerical stability when de contains only a small number of genes.

Value

The goana default method produces a data frame with a row for each GO term and the following columns:

• The last two column names above assume one gene set with the name DE . In general, there will be a pair of such columns for each gene set and the name of the set will appear in place of "DE" .

The goana method for MArrayLM objects produces a data frame with a row for each GO term and the following columns:

The row names of the data frame give the GO term IDs.

The output from kegga is the same except that row names become KEGG pathway IDs, Term becomes Pathway and there is no Ont column.

Seealso

topGO , topKEGG

The goseq package provides an alternative implementation of methods from Young et al (2010). Unlike the limma functions documented here, goseq will work with a variety of gene identifiers and includes a database of gene length information for various species.

The gostats package also does GO analyses without adjustment for bias but with some other options.

See 10.GeneSetTests for a description of other functions used for gene set testing.

Note

kegga requires an internet connection unless gene.pathway and pathway.names are both supplied.

The default for kegga with species="Dm" changed from convert=TRUE to convert=FALSE in limma 3.27.8. Users wanting to use Entrez Gene IDs for Drosophila should set convert=TRUE , otherwise fly-base CG annotation symbol IDs are assumed (for example "Dme1_CG4637").

The default for restrict.universe=TRUE in kegga changed from TRUE to FALSE in limma 3.33.4.

Bug fix: results from kegga with trend=TRUE or with non-NULL covariate were incorrect prior to limma 3.32.3. The results were biased towards significant Down p-values and against significant Up p-values.

Author

Gordon Smyth and Yifang Hu

References

Young, M. D., Wakefield, M. J., Smyth, G. K., Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology 11, R14. http://genomebiology.com/2010/11/2/R14

Examples

## Linear model usage:

fit <- lmFit(y, design)
fit <- eBayes(fit)

# Standard GO analysis

go.fisher <- goana(fit, species="Hs")
topGO(go.fisher, sort = "up")
topGO(go.fisher, sort = "down")

# GO analysis adjusting for gene abundance

go.abund <- goana(fit, geneid = "GeneID", trend = TRUE)
topGO(go.abund, sort = "up")
topGO(go.abund, sort = "down")

# GO analysis adjusting for gene length bias
# (assuming that y$genes$Length contains gene lengths)

go.len <- goana(fit, geneid = "GeneID", trend = "Length")
topGO(go.len, sort = "up")
topGO(go.len, sort = "down")

## Default usage with a list of gene sets:

go.de <- goana(list(DE1 = EG.DE1, DE2 = EG.DE2, DE3 = EG.DE3))
topGO(go.de, sort = "DE1")
topGO(go.de, sort = "DE2")
topGO(go.de, ontology = "BP", sort = "DE3")
topGO(go.de, ontology = "CC", sort = "DE3")
topGO(go.de, ontology = "MF", sort = "DE3")

## Standard KEGG analysis

k <- kegga(fit, species="Hs")
k <- kegga(fit, species.KEGG="hsa") # equivalent to previous
topKEGG(k, sort = "up")
topKEGG(k, sort = "down")

Row and Column Positions on Microarray

Description

Grid and spot row and column positions.

Usage

gridr(layout)
gridc(layout)
spotr(layout)
spotc(layout)

Arguments

Value

Vector of length prod(unlist(layout)) giving the grid rows ( gridr ), grid columns ( gridc ), spot rows ( spotr ) or spot columns ( spotc ).

Author

Gordon Smyth

Stemmed Heat Diagram

Description

Creates a heat diagram showing the co-regulation of genes under one condition with a range of other conditions.

Usage

heatDiagram(results, coef, primary=1, names=NULL, treatments=colnames(coef), limit=NULL,
            orientation="landscape", low="green", high="red", cex=1, mar=NULL,
            ncolors=123, list())
heatdiagram(stat, coef, primary=1, names=NULL, treatments=colnames(stat),
            critical.primary=4, critical.other=3, limit=NULL, orientation="landscape",
            low="green", high="red", cex=1, mar=NULL, ncolors=123, list())

Arguments

Details

Users are encouraged to use heatDiagram rather than heatdiagram as the later function may be removed in future versions of limma.

This function plots an image of gene expression profiles in which rows (or columns for portrait orientation) correspond to treatment conditions and columns (or rows) correspond to genes. Only genes which are significantly differentially expressed in the primary condition are included. Genes are sorted by differential expression under the primary condition.

Note: the plot produced by this function is unique to the limma package. It should not be confused with "heatmaps" often used to display results from cluster analyses.

Value

An image is created on the current graphics device. A matrix with named rows containing the coefficients used in the plot is also invisibly returned.

Seealso

image .

Author

Gordon Smyth

Examples

MA <- normalizeWithinArrays(RG)
design <- cbind(c(1,1,1,0,0,0),c(0,0,0,1,1,1))
fit <- lmFit(MA,design=design)
contrasts.mouse <- cbind(Control=c(1,0),Mutant=c(0,1),Difference=c(-1,1))
fit <- eBayes(contrasts.fit(fit,contrasts=contrasts.mouse))
results <- decideTests(fit,method="global",p=0.1)
heatDiagram(results,fit$coef,primary="Difference")

Prompt for Method Help Topics

Description

For any S4 generic function, find all methods defined in currently loaded packages. Prompt the user to choose one of these to display the help document.

Usage

helpMethods(genericFunction)

Arguments

Seealso

showMethods

Author

Gordon Smyth

Examples

helpMethods(show)

Convert Gene Identifiers to Indices for Gene Sets

Description

Make a list of gene identifiers into a list of indices for gene sets.

Usage

ids2indices(gene.sets, identifiers, remove.empty=TRUE)

Arguments

Details

This function used to create input for romer , mroast and camera function. Typically, identifiers is the vector of Entrez Gene IDs, and gene.sets is obtained constructed from a database of gene sets, for example a representation of the Molecular Signatures Database (MSigDB) downloaded from http://bioinf.wehi.edu.au/software/MSigDB .

Value

list of integer vectors, each vector containing the indices of a gene set in the vector identifiers .

Seealso

romer , mroast , camera

There is a topic page on 10.GeneSetTests .

Author

Gordon Smyth and Yifang Hu

Examples

download.file("http://bioinf.wehi.edu.au/software/MSigDB/human_c2_v5p2.rdata",
"human_c2_v5p2.rdata", mode = "wb")

load("human_c2_v5p2.rdata")
c2.indices <- ids2indices(Hs.c2, y$genes$GeneID)
camera(y, c2.indices, design)

Image Plot of Microarray Statistics

Description

Creates an image of colors or shades of gray that represent the values of a statistic for each spot on a spotted microarray. This function can be used to explore any spatial effects across the microarray.

Usage

imageplot(z, layout, low = NULL, high = NULL, ncolors = 123, zerocenter = NULL, 
zlim = NULL, mar=c(2,1,1,1), legend=TRUE, list())

Arguments

Details

This function may be used to plot the values of any spot-specific statistic, such as the log intensity ratio, background intensity or a quality measure such as spot size or shape. The image follows the layout of an actual microarray slide with the bottom left corner representing the spot (1,1,1,1). The color range is used to represent the range of values for the statistic. When this function is used to plot the red/green log-ratios, it is intended to be an in silico version of the classic false-colored red-yellow-green image of a scanned two-color microarray.

This function is related to the earlier plot.spatial function in the sma package and to the later maImage function in the marray package. It differs from plot.spatial most noticeably in that all the spots are plotted and the image is plotted from bottom left rather than from top left. It is intended to display spatial patterns and artefacts rather than to highlight only the extreme values as does plot.spatial . It differs from maImage in that any statistic may be plotted and in its use of a red-yellow-green color scheme for log-ratios, similar to the classic false-colored jpeg image, rather than the red-black-green color scheme associated with heat maps.

Value

An plot is created on the current graphics device.

Seealso

maImage in the marray package, image in the graphics package.

An overview of diagnostic functions available in LIMMA is given in 09.Diagnostics .

Author

Gordon Smyth

Examples

M <- rnorm(8*4*16*16)
imageplot(M,layout=list(ngrid.r=8,ngrid.c=4,nspot.r=16,nspot.c=16))

Write Imageplots to Files

Description

Write imageplots to files in PNG format, six plots to a file in a 3 by 2 grid arrangement.

Usage

imageplot3by2(RG, z="Gb", prefix=paste("image",z,sep="-"), path=NULL,
              zlim=NULL, common.lim=TRUE, list())

Arguments

Details

At the time of writing, this function writes plots in PNG format in an arrangement optimized for A4-sized paper.

Value

No value is returned, but one or more files are written to the working directory. The number of files is determined by the number of columns of RG .

Seealso

An overview of diagnostic functions available in LIMMA is given in 09.Diagnostics .

Author

Gordon Smyth

Intra-Spot Correlation for Two Color Data

Description

Estimate the within-block correlation associated with spots for spotted two color microarray data.

Usage

intraspotCorrelation(object, design, trim=0.15)

Arguments

Details

This function estimates the correlation between two channels observed on each spot. The correlation is estimated by fitting a heteroscedastic regression model to the M and A-values of each gene. The function also returns a consensus correlation, which is a robust average of the individual correlations, which can be used as input for functions lmscFit .

The function may take long time to execute.

Value

A list with components

*

Seealso

This function uses remlscore from the statmod package.

An overview of methods for single channel analysis in limma is given by 07.SingleChannel .

Author

Gordon Smyth

References

Smyth, G. K. (2005). Individual channel analysis of two-colour microarray data. Proceedings of the 55th Session of the International Statistics Institute , 5-12 April 2005, Sydney, Australia, Paper 116. http://www.statsci.org/smyth/pubs/ISI2005-116.pdf

Examples

#  See lmscFit
corfit <- intraspotCorrelation(MA, design)
all.correlations <- tanh(corfit$atanh.correlations)
boxplot(all.correlations)

Check for Full Column Rank

Description

Test whether a numeric matrix has full column rank.

Usage

is.fullrank(x)
nonEstimable(x)

Arguments

Details

is.fullrank is used to check the integrity of design matrices in limma, for example after subsetting operations.

nonEstimable is used by lmFit to report which coefficients in a linear model cannot be estimated.

Value

is.fullrank returns TRUE or FALSE .

nonEstimable returns a character vector of names for the columns of x which are linearly dependent on previous columns. If x has full column rank, then the value is NULL .

Author

Gordon Smyth

Examples

# TRUE
is.fullrank(1)
is.fullrank(cbind(1,0:1))

# FALSE
is.fullrank(0)
is.fullrank(matrix(1,2,2))
nonEstimable(matrix(1,2,2))

Test for Numeric Argument

Description

Test whether argument is numeric or a data.frame with numeric columns.

Usage

isNumeric(x)

Arguments

Details

This function is used to check the validity of arguments for numeric functions. It is an attempt to emulate the behavior of internal generic math functions.

isNumeric differs from is.numeric in that data.frames with all columns numeric are accepted as numeric.

Value

TRUE or FALSE

Seealso

is.numeric , Math

Author

Gordon Smyth

Examples

isNumeric(3)
isNumeric("a")
x <- data.frame(a=c(1,1),b=c(0,1))
isNumeric(x)   # TRUE
is.numeric(x)  # FALSE

Kooperberg Model-Based Background Correction for GenePix data

Description

This function uses a Bayesian model to background correct GenePix microarray data.

Usage

kooperberg(RG, a = TRUE, layout = RG$printer, verbose = TRUE)

Arguments

Details

This function is for use with GenePix data and is designed to cope with the problem of large numbers of negative intensities and hence missing values on the log-intensity scale. It avoids missing values in most cases and at the same time dampens down the variability of log-ratios for low intensity spots. See Kooperberg et al (2002) for more details.

kooperberg uses the foreground and background intensities, standard deviations and number of pixels to compute empirical estimates of the model parameters as described in equation 2 of Kooperberg et al (2002).

Value

An RGList containing the components

*

Seealso

04.Background gives an overview of background correction functions defined in the LIMMA package.

Author

Matthew Ritchie

References

Kooperberg, C., Fazzio, T. G., Delrow, J. J., and Tsukiyama, T. (2002) Improved background correction for spotted DNA microarrays. Journal of Computational Biology 9 , 55-66.

Ritchie, M. E., Silver, J., Oshlack, A., Silver, J., Holmes, M., Diyagama, D., Holloway, A., and Smyth, G. K. (2007). A comparison of background correction methods for two-colour microarrays. Bioinformatics 23, 2700-2707. https://www.ncbi.nlm.nih.gov/pubmed/17720982

Examples

#  This is example code for reading and background correcting GenePix data
#  given GenePix Results (gpr) files in the working directory (data not
#  provided).
# get the names of the GenePix image analysis output files in the current directory
genepixFiles <- dir(pattern="*\.gpr$")
RG <- read.maimages(genepixFiles, source="genepix", other.columns=c("F635 SD","B635 SD",
"F532 SD","B532 SD","B532 Mean","B635 Mean","F Pixels","B Pixels"))
RGmodel <- kooperberg(RG)
MA <- normalizeWithinArrays(RGmodel)

View Limma User's Guide

Description

Finds the location of the Limma User's Guide and optionally opens it.

Usage

limmaUsersGuide(view=TRUE)

Arguments

Details

The function vignette("limma") will find the short limma Vignette which describes how to obtain the Limma User's Guide. The User's Guide is not itself a true vignette because it is not automatically generated using Sweave during the package build process. This means that it cannot be found using vignette , hence the need for this special function.

If the operating system is other than Windows, then the PDF viewer used is that given by Sys.getenv("R_PDFVIEWER") . The PDF viewer can be changed using Sys.putenv(R_PDFVIEWER=) .

This function is used by drop-down Vignettes menu when the Rgui interface for Windows is used.

Value

Character string giving the file location.

Seealso

vignette , openPDF , openVignette , Sys.getenv , Sys.putenv

Author

Gordon Smyth

Examples

limmaUsersGuide(view=FALSE)

Linear Model for Series of Arrays

Description

Fit linear model for each gene given a series of arrays

Usage

lmFit(object, design=NULL, ndups=1, spacing=1, block=NULL, correlation, weights=NULL,
      method="ls", list())

Arguments

Details

This function fits multiple linear models by weighted or generalized least squares. It accepts data from a experiment involving a series of microarrays with the same set of probes. A linear model is fitted to the expression data for each probe. The expression data should be log-ratios for two-color array platforms or log-expression values for one-channel platforms. (To fit linear models to the individual channels of two-color array data, see lmscFit .) The coefficients of the fitted models describe the differences between the RNA sources hybridized to the arrays. The probe-wise fitted model results are stored in a compact form suitable for further processing by other functions in the limma package.

The function allows for missing values and accepts quantitative precision weights through the weights argument. It also supports two different correlation structures. If block is not NULL then different arrays are assumed to be correlated. If block is NULL and ndups is greater than one then replicate spots on the same array are assumed to be correlated. It is not possible at this time to fit models with both a block structure and a duplicate-spot correlation structure simultaneously.

If object is a matrix then it should contain log-ratios or log-expression data with rows corresponding to probes and columns to arrays. (A numeric vector is treated the same as a matrix with one column.) For objects of other classes, a matrix of expression values is taken from the appropriate component or slot of the object. If object is of class MAList or marrayNorm , then the matrix of log-ratios (M-values) is extracted. If object is of class ExpressionSet , then the expression matrix is extracted. (This may contain log-expression or log-ratio values, depending on the platform.) If object is of class PLMset then the matrix of chip coefficients chip.coefs is extracted.

The arguments design , ndups , spacing and weights will be extracted from the data object if available and do not normally need to set explicitly in the call. On the other hand, if any of these are set in the function call then they will over-ride the slots or components in the data object . If object is an PLMset , then weights are computed as 1/pmax(object@se.chip.coefs, 1e-05)^2 . If object is an ExpressionSet object, then weights are not computed.

If the argument block is used, then it is assumed that ndups=1 .

The correlation argument has a default value of 0.75 , but in normal use this default value should not be relied on and the correlation value should be estimated using the function duplicateCorrelation . The default value is likely to be too high in particular if used with the block argument.

The actual linear model computations are done by passing the data to one the lower-level functions lm.series , gls.series or mrlm . The function mrlm is used if method="robust" . If method="ls" , then gls.series is used if a correlation structure has been specified, i.e., if ndups>1 or block is non-null and correlation is different from zero. If method="ls" and there is no correlation structure, lm.series is used.

Value

An MArrayLM object containing the result of the fits.

The rownames of object are preserved in the fit object and can be retrieved by rownames(fit) where fit is output from lmFit . The column names of design are preserved as column names and can be retrieved by colnames(fit) .

Seealso

lmFit uses getEAWP to extract expression values, gene annotation and so from the data object .

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Gordon Smyth

Examples

# Simulate gene expression data for 100 probes and 6 microarrays
# Microarray are in two groups
# First two probes are differentially expressed in second group
# Std deviations vary between genes with prior df=4
sd <- 0.3*sqrt(4/rchisq(100,df=4))
y <- matrix(rnorm(100*6,sd=sd),100,6)
rownames(y) <- paste("Gene",1:100)
y[1:2,4:6] <- y[1:2,4:6] + 2
design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
options(digits=3)

# Ordinary fit
fit <- lmFit(y,design)
fit <- eBayes(fit)
topTable(fit,coef=2)
dim(fit)
colnames(fit)
rownames(fit)[1:10]
names(fit)

# Fold-change thresholding
fit2 <- treat(fit,lfc=0.1)
topTreat(fit2,coef=2)

# Volcano plot
volcanoplot(fit,coef=2,highlight=2)

# Mean-difference plot
plotMD(fit,column=2)

# Q-Q plot of moderated t-statistics
qqt(fit$t[,2],df=fit$df.residual+fit$df.prior)
abline(0,1)

# Various ways of writing results to file
write.fit(fit,file="exampleresults.txt")
write.table(fit,file="exampleresults2.txt")

# Fit with correlated arrays
# Suppose each pair of arrays is a block
block <- c(1,1,2,2,3,3)
dupcor <- duplicateCorrelation(y,design,block=block)
dupcor$consensus.correlation
fit3 <- lmFit(y,design,block=block,correlation=dupcor$consensus)

# Fit with duplicate probes
# Suppose two side-by-side duplicates of each gene
rownames(y) <- paste("Gene",rep(1:50,each=2))
dupcor <- duplicateCorrelation(y,design,ndups=2)
dupcor$consensus.correlation
fit4 <- lmFit(y,design,ndups=2,correlation=dupcor$consensus)
dim(fit4)
fit4 <- eBayes(fit4)
topTable(fit4,coef=2)

Fit Linear Model to Individual Channels of Two-Color Data

Description

Fit a linear model to the individual log-intensities for each gene given a series of two-color arrays

Usage

lmscFit(object, design, correlation)

Arguments

Details

For two color arrays, the channels measured on the same set of arrays are correlated. The M and A however are uncorrelated for each gene. This function fits a linear model to the set of M and A-values for each gene after re-scaling the M and A-values to have equal variances. The input correlation determines the scaling required. The input correlation is usually estimated using intraspotCorrelation before using lmscFit .

Missing values in M or A are not allowed.

Value

An object of class MArrayLM

Seealso

lm.fit .

An overview of methods for single channel analysis in limma is given by 07.SingleChannel .

Author

Gordon Smyth

References

Smyth, GK (2005). Individual channel analysis of two-colour microarray data. Proceedings of the 55th Session of the International Statistics Institute , 5-12 April 2005, Sydney, Australia; Internatational Statistics Institute; Paper 116. http://www.statsci.org/smyth/pubs/ISI2005-116.pdf

Smyth, GK, and Altman, NS (2013). Separate-channel analysis of two-channel microarrays: recovering inter-spot information. BMC Bioinformatics 14, 165. http://www.biomedcentral.com/1471-2105/14/165

Examples

# Subset of data from ApoAI case study in Limma User's Guide
# Avoid non-positive intensities
RG <- backgroundCorrect(RG,method="normexp")
MA <- normalizeWithinArrays(RG)
MA <- normalizeBetweenArrays(MA,method="Aq")
targets <- data.frame(Cy3=I(rep("Pool",6)),Cy5=I(c("WT","WT","WT","KO","KO","KO")))
targets.sc <- targetsA2C(targets)
targets.sc$Target <- factor(targets.sc$Target,levels=c("Pool","WT","KO"))
design <- model.matrix(~Target,data=targets.sc)
corfit <- intraspotCorrelation(MA,design)
fit <- lmscFit(MA,design,correlation=corfit$consensus)
cont.matrix <- cbind(KOvsWT=c(0,-1,1))
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
topTable(fit2,adjust="fdr")

Fit Linear Model to Microrray Data by Ordinary Least Squares

Description

Fit a linear model genewise to expression data from a series of arrays. This function uses ordinary least squares and is a utility function for lmFit .

Usage

lm.series(M,design=NULL,ndups=1,spacing=1,weights=NULL)

Arguments

Details

This is a utility function used by the higher level function lmFit . Most users should not use this function directly but should use lmFit instead.

The linear model is fit for each gene by calling the function lm.fit or lm.wfit from the base library.

Value

A list with components

*

Seealso

lm.fit .

An overview of linear model functions in limma is given by 06.LinearModels .

Author

Gordon Smyth

Examples

# See lmFit for examples

Univariate Lowess With Prior Weights

Description

Univariate locally weighted linear regression allowing for prior weights. Returns fitted values and residuals.

Usage

loessFit(y, x, weights=NULL, span=0.3, iterations=4L, min.weight=1e-5, max.weight=1e5,
         equal.weights.as.null=TRUE, method="weightedLowess")

Arguments

Details

This function is essentially a wrapper function for lowess and weightedLowess with added error checking. The idea is to provide the classic univariate lowess algorithm of Cleveland (1979) but allowing for prior weights and missing values.

The venerable lowess code is fast, uses little memory and has an accurate interpolation scheme, so it is an advantage to use it when prior weights are not needed. This functions calls lowess when weights=NULL , but returns values in original rather than sorted order and allows missing values. The treatment of missing values is analogous to na.exclude .

By default, weights that are all equal (even all zero) are treated as if they were NULL , so lowess is called in this case also.

When unequal weights are provided, this function calls weightedLowess by default, although two other possibilities are also provided. weightedLowess implements a similar algorithm to lowess except that it uses the prior weights both in the local regressions and in determining which other observations to include in the local neighbourhood of each observation.

Two alternative algorithms for weighted lowess curve fitting are provided as options. If method="loess" , then a call is made to loess(y~x,weights=weights,span=span,degree=1,family="symmetric", . This method differs from weightedLowess in that the prior weights are ignored when determining the neighbourhood of each observation.

If method="locfit" , then repeated calls are made to locfit:::locfit.raw with deg=1 . In principle, this is similar to "loess" , but "locfit" makes some approximations and is very much faster and uses much less memory than "loess" for long data vectors.

The arguments span and iterations here have the same meaning as for weightedLowess and loess . span is equivalent to the argument f of lowess while iterations is equivalent to iter+1 for lowess . It gives the total number of fits rather than the number of robustifying fits.

When there are insufficient observations to estimate the loess curve, loessFit returns a linear regression fit. This mimics the behavior of lowess but not that of loess or locfit.raw .

Value

A list with components

*

Seealso

If weights=NULL , this function calls lowess . Otherwise it calls weightedLowess , locfit.raw or loess . See the help pages of those functions for references and credits.

Compare with loess in the stats package.

See 05.Normalization for an outline of the limma package normalization functions.

Note

With unequal weights, "loess" was the default method prior to limma version 3.17.25. The default was changed to "locfit" in limma 3.17.25, and then to "weightedLowess" in limma 3.19.16. "weightedLowess" will potentially give somewhat different results to the older algorithms because the local neighbourhood of each observation is determined differently (more carefully).

Author

Gordon Smyth

References

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association 74, 829-836.

Examples

x <- (1:100)/101
y <- sin(2*pi*x)+rnorm(100,sd=0.4)
out <- loessFit(y,x)
plot(x,y)
lines(x,out$fitted,col="red")

# Example using weights

y <- x-0.5
w <- rep(c(0,1),50)
y[w==0] <- rnorm(50,sd=0.1)
pch <- ifelse(w>0,16,1)
plot(x,y,pch=pch)
out <- loessFit(y,x,weights=w)
lines(x,out$fitted,col="red")

Logarithm of cosh

Description

Compute log(cosh(x)) without floating overflow or underflow

Usage

logcosh(x)

Arguments

Details

The computation uses asymptotic expressions for very large or very small arguments. For intermediate arguments, log(cosh(x)) is returned.

Value

Numeric vector or matrix of same dimensions as x .

Seealso

logsumexp

Author

Gordon K Smyth

Examples

x <- c(1e-8,1e-7,1e-6,1e-5,1e-4,1,3,50,800)
logcosh(x)
log(cosh(x))

Log Sum of Exponentials

Description

Compute log( exp(x)+exp(y) ) without floating overflow or underflow

Usage

logsumexp(x, y)

Arguments

Details

The computation uses logcosh() .

Value

Numeric vector or matrix of same dimensions as x .

Seealso

logcosh