Saving plot to file

In order to save the plot to a file, use the ggsave() function:

fig <- plotVarPart(vp)
ggsave(file, fig)

Plot expression stratified by other variables

variancePartition also includes plotting functions to visualize the variation across a variable of interest. plots the expression of a gene stratified by the specified variable. In the example dataset, users can plot a gene expression trait stratified by Tissue (Figure 2a) or Individual (Figure 2b).

# get gene with the highest variation across Tissues
# create data.frame with expression of gene i and Tissue
#       type for each sample
i <- which.max(varPart$Tissue)
GE <- data.frame(Expression = geneExpr[i, ], Tissue = info$Tissue)

# Figure 2a
# plot expression stratified by Tissue
plotStratify(Expression ~ Tissue, GE, main = rownames(geneExpr)[i])

# get gene with the highest variation across Individuals
# create data.frame with expression of gene i and Tissue
#       type for each sample
i <- which.max(varPart$Individual)
GE <- data.frame(
  Expression = geneExpr[i, ],
  Individual = info$Individual
)

# Figure 2b
# plot expression stratified by Tissue
label <- paste("Individual:", format(varPart$Individual[i] * 100,
  digits = 3
), "%")
main <- rownames(geneExpr)[i]
plotStratify(Expression ~ Individual, GE,
  colorBy = NULL,
  text = label, main = main
)

For gene141, variation across tissues explains 52.9% of variance in gene expression. For gene43, variation across Individuals explains 91.4% of variance in gene expression.

Intuition about the backend

At the heart of variancePartition, a regression model is fit for each gene separately and summary statistics are extracted and reported to the user for visualization and downstream analysis. For a single model fit, calcVarPart() computes the fraction of variance explained by each variable. calcVarPart() is defined by this package, and computes these statistics from either a fixed effects model fit with lm() or a linear mixed model fit with lme4::lmer(). fitExtractVarPartModel() loops over each gene, fits the regression model and returns the variance fractions reported by calcVarPart().

Fitting the regression model and extracting variance statistics can also be done directly:

library("lme4")

# fit regression model for the first gene
form_test <- geneExpr[1, ] ~ Age + (1 | Individual) + (1 | Tissue)
fit <- lmer(form_test, info, REML = FALSE)

# extract variance statistics
calcVarPart(fit)

## Individual     Tissue        Age  Residuals 
##   8.90e-01   2.47e-02   4.35e-05   8.50e-02

Should a variable be modeled as fixed or random effect?

Categorical variables should (almost) always be modeled as a random effect. The difference between modeling a categorical variable as a fixed versus random effect is minimal when the sample size is large compared to the number of categories (i.e. levels). So variables like disease status, sex or time point will not be sensitive to modeling as a fixed versus random effect. However, variables with many categories like Individual must be modeled as a random effect in order to obtain statistically valid results. So to be on the safe side, categorical variable should be modeled as a random effect.

% R and variancePartition handle catagorical variables stored as a very naturally. If categorical variables are stored as an or , they must be converted to a before being used with variancePartition

variancePartition fits two types of models:

1. linear mixed model where all categorical variables are modeled as random effects and all continuous variables are fixed effects. The function lme4::lmer() is used to fit this model.

2. fixed effected model, where all variables are modeled as fixed effects. The function lm() is used to fit this model.

Which variables should be included?

In my experience, it is useful to include all variables in the first analysis and then drop variables that have minimal effect. However, like all multiple regression methods, variancePartition will divide the contribution over multiple variables that are strongly correlated. So, for example, including both sex and height in the model will show sex having a smaller contribution to variation gene expression than if height were omitted, since there variables are strongly correlated. This is a simple example, but should give some intuition about a common issue that arises in analyses with variancePartition.

variancePartition can naturally assess the contribution of both individual and sex in a dataset. As expected, genes for which sex explains a large fraction of variation are located on chrX and chrY. If the goal is to interpret the impact of sex, then there is no issue. But recall the issue with correlated variables and note that individual is correlated with sex, because each individual is only one sex regardless of how many samples are taken from a individual. It follows that including sex in the model reduces the apparent contribution of individual to gene expression. In other words, the ICC for individual will be different if sex is included in the model.

In general, including variables in the model that do not vary within individual will reduce the apparent contribution of individual as estimated by variancePartition. For example, sex and ethnicity never vary between multiple samples from the same individual and will always reduce the apparent contribution of individual. However, disease state and age may or may not vary depending on the study design.

In biological datasets technical variability (i.e. batch effects) can often reduce the apparent biological signal. In RNA-seq analysis, it is common for the the impact of this technical variability to be removed before downstream analysis. Instead of including these batch variable in the variancePartition analysis, it is simple to complete the expression residuals with the batch effects removed and then feeds these residuals to variancePartition. This will increase the fraction of variation explained by biological variables since technical variability is reduced.

form <- ~ Individual + Tissue + Batch + Age + Height

# Compute Canonical Correlation Analysis (CCA)
# between all pairs of variables
# returns absolute correlation value
C <- canCorPairs(form, info)

# Plot correlation matrix
# between all pairs of variables
plotCorrMatrix(C)

Extracting additional information from model fits

Advanced users may want to perform the model fit and extract results in separate steps in order to examine the fit of the model for each gene. Thus the work of can be divided into two steps: 1) fit the regression model, and 2) extracting variance statistics.

form <- ~ Age + (1 | Individual) + (1 | Tissue) + (1 | Batch)

# Fit model
results <- fitVarPartModel(geneExpr, form, info)

# Extract results
varPart <- extractVarPart(results)

Note that storing the model fits can use a lot of memory (~10Gb with 20K genes and 1000 experiments). I do not recommend unless you have a specific need for storing the entire model fit.

Instead, fitVarPartModel() can extract any desired information using any function that accepts the model fit from lm() or lmer(). The results are stored in a list and can be used for downstream analysis.

# Fit model and run summary() function on each model fit
vpSummaries <- fitVarPartModel(geneExpr, form, info, fxn = summary)

# Show results of summary() for the first gene
vpSummaries[[1]]

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: y.local ~ Age + (1 | Individual) + (1 | Tissue) + (1 | Batch)
##    Data: data
## Weights: data$w.local
## Control: control
## 
##      AIC      BIC   logLik deviance df.resid 
##      397      413     -193      385       94 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0580 -0.5865  0.0147  0.6603  1.9708 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Individual (Intercept) 10.82275 3.2898  
##  Batch      (Intercept)  0.00192 0.0438  
##  Tissue     (Intercept)  0.30010 0.5478  
##  Residual                1.02997 1.0149  
## Number of obs: 100, groups:  Individual, 25; Batch, 4; Tissue, 3
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) -10.60243    1.09493   -9.68
## Age           0.00318    0.01610    0.20
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age -0.739

Removing batch effects before fitting model

Gene expression studies often have substantial batch effects, and variancePartition can be used to understand the magnitude of the effects. However, we often want to focus on biological variables (i.e. individual, tissue, disease, sex) after removing the effect of technical variables. Depending on the size of the batch effect, I have found it useful to correct for the batch effect first and then perform a variancePartition analysis afterward. Subtracting this batch effect can reduce the total variation in the data, so that the contribution of other variables become clearer.

Standard analysis:

form <- ~ (1 | Tissue) + (1 | Individual) + (1 | Batch) + Age
varPart <- fitExtractVarPartModel(geneExpr, form, info)

Analysis on residuals:

library("limma")
# subtract out effect of Batch
fit <- lmFit(geneExpr, model.matrix(~Batch, info))
res <- residuals(fit, geneExpr)

# fit model on residuals
form <- ~ (1 | Tissue) + (1 | Individual) + Age

varPartResid <- fitExtractVarPartModel(res, form, info)

Remove batch effect with linear mixed model

# subtract out effect of Batch with linear mixed model
modelFit <- fitVarPartModel(geneExpr, ~ (1 | Batch), info)
res <- residuals(modelFit)

# fit model on residuals
form <- ~ (1 | Tissue) + (1 | Individual) + Age

varPartResid <- fitExtractVarPartModel(res, form, info)

If the two-step process requires too much memory, the residuals can be computed more efficiently. Here, run the function inside the call to fitVarPartModel() to avoid storing the large intermediate results.

# extract residuals directly without storing intermediate results
residList <- fitVarPartModel(geneExpr, ~ (1 | Batch), info,
  fxn = residuals
)

# convert list to matrix
residMatrix <- do.call(rbind, residList)

Variation within multiple subsets of the data

So far, we have focused on interpreting one variable at a time. But the linear mixed model behind variancePartition is a very powerful framework for analyzing variation at multiple levels. We can easily extend the previous analysis of the contribution of individual and tissue on variation in gene expression to examine the contribution of individual within each tissue. This analysis is as easy as specifying a new formula and rerunning variancePartition. Note that is analysis will only work when there are replicates for at least some individuals within each tissue in order to assess cross-individual variance with in a tissue.

# specify formula to model within/between individual variance
# separately for each tissue
# Note that including +0 ensures each tissue is modeled explicitly
# Otherwise, the first tissue would be used as baseline
form <- ~ (Tissue + 0 | Individual) + Age + (1 | Tissue) + (1 | Batch)

# fit model and extract variance percents
varPart <- fitExtractVarPartModel(geneExpr, form, info, showWarnings = FALSE)

# violin plot
plotVarPart(sortCols(varPart), label.angle = 60)

This analysis corresponds to a varying coefficient model, where the correlation between individuals varies for each tissue ]. Since the variation across individuals is modeled within each tissue, the total variation explained does not sum to 1 as it does for standard application of variancePartition. So interpretation as intra-class does not strictly apply and use of plotPercentBars() is no longer applicable. Yet the variables in the study design are still ranked in terms of their genome-wide contribution to expression variation, and results can still be analyzed at the gene level. See Variation with multiple subsets of the data for statistical details.

Detecting problems caused by collinearity of variables

Including variables that are highly correlated can produce misleading results and overestimate the contribution of variables modeled as fixed effects. This is usually not an issue, but can arise when statistically redundant variables are included in the model. In this case, the model is "degenerate" or "computationally singular" and parameter estimates from this model are not meaningful. Dropping one or more of the covariates will fix this problem.

A check of collinearity is built into fitVarPartModel() and fitExtractVarPartModel(), so the user will be warned if this is an issue.

Alternatively, the user can use the colinearityScore() function to evaluate whether this is an issue for a single model fit:

form <- ~ (1 | Individual) + (1 | Tissue) + Age + Height

# fit model
res <- fitVarPartModel(geneExpr[1:4, ], form, info)

# evaluate the collinearity score on the first model fit
# this reports the correlation matrix between coefficient estimates
# for fixed effects
# the collinearity score is the maximum absolute correlation value
# If the collinearity score > .99 then the variance partition
# estimates may be problematic
# In that case, a least one variable should be omitted
colinearityScore(res[[1]])

## [1] 0.777
## attr(,"vcor")
##             (Intercept)     Age  Height
## (Intercept)       1.000 -0.4191 -0.7774
## Age              -0.419  1.0000 -0.0575
## Height           -0.777 -0.0575  1.0000

Including weights computed separately

variancePartition automatically used precision weights computed by limma::voom(), but the user can also specify custom weights using the weightsMatrix argument.

form <- ~ (1 | Individual) + (1 | Tissue) + Age + Height

# Specify custom weights
# In this example the weights are simulated from a
# uniform distribution and are not meaningful.
weights <- matrix(runif(length(geneExpr)), nrow = nrow(geneExpr))

# Specify custom weights
res <- fitExtractVarPartModel(geneExpr[1:4, ], form, info,
  weightsMatrix = weights[1:4, ]
)

In addition, setting the useWeights=FALSE will suppress usage of the weights in all cases, i.e. when the weights are specified manually or implicitly with the results of limma::voom().

Including interaction terms

Typical analysis assumes that the effect of each variable on gene expression does not depend on other variables in the model. Sometimes this assumption is too strict, and we want to model an interaction effect whereby the effect of Batch depends on Tissue. This can be done easly by specifying an interaction term, (1|Batch:Tissue). Since Batch has 4 categories and Tissue has 3, this interaction term implicity models a new 3*4 = 12 category variable in the analysis. This new interaction term will absorb some of the variance from the Batch and Tissue term, so an interaction model should always include the two constituent variables.

Here we fit an interaction model, but we observe that interaction between Batch and Tissue does not explain much expression variation.

form <- ~ (1 | Individual) + Age + Height + (1 | Tissue) + (1 | Batch) +
  (1 | Batch:Tissue)

# fit model
vpInteraction <- fitExtractVarPartModel(geneExpr, form, info)

plotVarPart(sortCols(vpInteraction))

limma::voom()

Read RNA-seq counts into R, normalize for library size within and between experiments with TMM (Robinson and Oshlack 2010), estimate precision weights with limma::voom().

library("limma")
library("edgeR")

# identify genes that pass expression cutoff
isexpr <- rowSums(cpm(geneCounts) > 1) >= 0.5 * ncol(geneCounts)

# create data structure with only expressed genes
gExpr <- DGEList(counts = geneCounts[isexpr, ])

# Perform TMM normalization
gExpr <- calcNormFactors(gExpr)

# Specify variables to be included in the voom() estimates of
# uncertainty.
# Recommend including variables with a small number of categories
# that explain a substantial amount of variation
design <- model.matrix(~Batch, info)

# Estimate precision weights for each gene and sample
# This models uncertainty in expression measurements
vobjGenes <- voom(gExpr, design)

# Define formula
form <- ~ (1 | Individual) + (1 | Tissue) + (1 | Batch) + Age

# variancePartition seamlessly deals with the result of voom()
# by default, it seamlessly models the precision weights
# This can be turned off with useWeights=FALSE
varPart <- fitExtractVarPartModel(vobjGenes, form, info)

DESeq2

Process and normalize the gene-level counts before running variancePartition analysis.

library("DESeq2")

# create DESeq2 object from gene-level counts and metadata
dds <- DESeqDataSetFromMatrix(
  countData = geneCounts,
  colData = info,
  design = ~1
)

# Estimate library size correction scaling factors
dds <- estimateSizeFactors(dds)

# identify genes that pass expression cutoff
isexpr <- rowSums(fpm(dds) > 1) >= 0.5 * ncol(dds)

# compute log2 Fragments Per Million
# Alternatively, fpkm(), vst() or rlog() could be used
quantLog <- log2(fpm(dds)[isexpr, ] + 1)

# Define formula
form <- ~ (1 | Individual) + (1 | Tissue) + (1 | Batch) + Age

# Run variancePartition analysis
varPart <- fitExtractVarPartModel(quantLog, form, info)

Note that DESeq2 does not compute precision weights like limma::voom(), so they are not used in this version of the analysis.

Isoform quantification

Other software performs isoform-level quantification using reads that map to multiple transcripts. These include kallisto (Bray et al. 2016), sailfish (Patro, Mount, and Kingsford 2014), salmon Patro, Duggal, and Kingsford (2015), RSEM (Li and Dewey 2011)and cufflinks (Trapnell et al. 2010.)

library("tximportData")
library("tximport")
library("readr")

# Get data from folder where tximportData is installed
dir <- system.file("extdata", package = "tximportData")
samples <- read.table(file.path(dir, "samples.txt"), header = TRUE)
files <- file.path(dir, "kallisto", samples$run, "abundance.tsv")
names(files) <- paste0("sample", 1:6)

tx2gene <- read.csv(file.path(dir, "tx2gene.csv"))

# reads results from kallisto
txi <- tximport(files,
  type = "kallisto", tx2gene = tx2gene,
  countsFromAbundance = "lengthScaledTPM"
)

# define metadata (usually read from external source)
info_tximport <- data.frame(
  Sample = sprintf("sample%d", 1:6),
  Disease = c("case", "control")[c(rep(1, 3), rep(2, 3))]
)

# Extract counts from kallisto
y <- DGEList(txi$counts)

# compute library size normalization
y <- calcNormFactors(y)

# apply voom to estimate precision weights
design <- model.matrix(~Disease, data = info_tximport)
vobj <- voom(y, design)

# define formula
form <- ~ (1 | Disease)

# Run variancePartition analysis (on only 10 genes)
varPart_tx <- fitExtractVarPartModel(
  vobj[1:10, ], form,
  info_tximport
)

ballgown

Quantifications from Cufflinks/Tablemaker and RSEM can be processed and read into R with the Bioconductor package ballgown.

library("ballgown")

# Get data from folder where ballgown is installed
data_directory <- system.file("extdata", package = "ballgown")

# Load results of Cufflinks/Tablemaker
bg <- ballgown(
  dataDir = data_directory, samplePattern = "sample",
  meas = "all"
)

# extract gene-level FPKM quantification
# Expression can be convert to log2-scale if desired
gene_expression <- gexpr(bg)

# extract transcript-level FPKM quantification
# Expression can be convert to log2-scale if desired
transcript_fpkm <- texpr(bg, "FPKM")

# define metadata (usually read from external source)
info_ballgown <- data.frame(
  Sample = sprintf("sample%02d", 1:20),
  Batch = rep(letters[1:4], 5),
  Disease = c("case", "control")[c(rep(1, 10), rep(2, 10))]
)

# define formula
form <- ~ (1 | Batch) + (1 | Disease)

# Run variancePartition analysis
# Gene-level analysis
varPart_gene <- fitExtractVarPartModel(
  gene_expression, form,
  info_ballgown
)

## Warning in filterInputData(exprObj, formula, data, useWeights = useWeights): Sample names of responses (i.e. columns of exprObj) do not match
## sample names of metadata (i.e. rows of data).  Recommend consistent
## names so downstream results are labeled consistently.

## Warning in .fitExtractVarPartModel(exprObj, formula, data, REML = REML, : Model failed for 1 responses.
##   See errors with attr(., 'errors')

# Transcript-level analysis
varPart_transcript <- fitExtractVarPartModel(
  transcript_fpkm, form,
  info_ballgown
)

## Warning in filterInputData(exprObj, formula, data, useWeights = useWeights): Sample names of responses (i.e. columns of exprObj) do not match
## sample names of metadata (i.e. rows of data).  Recommend consistent
## names so downstream results are labeled consistently.

## Warning in .fitExtractVarPartModel(exprObj, formula, data, REML = REML, : Model failed for 2 responses.
##   See errors with attr(., 'errors')

Note that ballgownrsem can be used for a similar analysis of RSEM results.

See tutorial for more details.

Tissue is major source of variation

Individual is major source of variation

Implementation in R

An R formula is used to define the terms in the fixed and random effects, and fitVarPartModel() fits the specified model for each gene separately. If random effects are specified, lme4::lmer() is used behind the scenes to fit the model, while lm() is used if there are only fixed effects. fitVarPartModel() returns a list of the model fits, and returns the variance partition statistics for each model in the list. fitExtractVarPartModel() combines the actions of fitVarPartModel() and into one function call. calcVarPart() is called behind the scenes to compute variance fractions for both fixed and mixed effects models, but the user can also call this function directly on a model fit with lm() or lmer().

Interpretation of percent variance explained

The percent variance explained can be interpreted as the intra-class correlation (ICC) when a special case of Equation 1 is used. Consider the simplest example of the i^th sample from the k^th individual with only an intercept term, one random effect corresponding to individual, and an error term. In this case ICC corresponds to the correlation between two samples from the same individual. This value is equal to the fraction of variance explained by individual. For example, consider the correlation between samples from the same individual: The correlation between samples from different individuals is: This interpretation in terms of fraction of variation explained (FVE) naturally generalizes to multiple variance components. Consider two sources of variation, individual and cell type with variances σ_id² and σ_cell², respectively. Applying a generalization of the the previous derivation, two samples are correlated according to:

Notice that the correlation between samples from the same individual and same cell type corresponds to the sum of the fraction explained by individual + fraction explained by cell type. This defines ICC for individual and tissue, as well as the combined ICC and relates these values to FVE.

In order to illustrate how this FVE and ICC relate to the correlation between samples in multilevel datasets, consider a simple example of 5 samples from 2 individuals and 2 tissues:

Modeling the separate effects of individual and tissue gives the following covariance structure between samples when a linear mixed model is used:

The covariance matrix is symmetric so that blank entries take the value on the opposite side of the diagonal. The covariance can be converted to correlation by dividing by σ_id² + σ_cell² + σ_ε², and this gives the results from above. This example generalizes to any number of variance components (Pinheiro and Bates 2000).

Variation with multiple subsets of the data

The linear mixed model underlying variancePartition allows the effect of one variable to depend on the value of another variable. Statistically, this is called a varying coefficient model (Pinheiro and Bates 2000; Galecki and Burzykowski 2013). This model arises in variancePartition analysis when the variation explained by individual depends on tissue or cell type.

A given sample is only from one cell type, so this analysis asks a question about a subset of the data. The the data is implicitly divided into subsets base on cell type and variation explained by individual is evaluated within each subset. The data is not actually divided onto subset, but the statistical model essentially examples samples with each cell type. This subsetting means that the variance fractions do not sum to 1.

Consider a concrete example with variation from across individual and cell types (T-cells and monocytes) with data from the i^th sample from the k^th individual, sex of s and cell type c. Modeling the variation across individuals within cell type corresponds to

with corresponding variance components:

Since the dataset is now divided into multiple subsets, direct interpretation of the fraction of variation explained (FVE) as intra-class correlation does not apply. Instead, we compute a “pseudo-FVE” by approximating the total variance attributable to cell type by using a weighted average of the within cell type variances weighted by the sample size within each cell type. Thus the values of pseudo-FVE do not have the simple interpretation as in the standard application of variancePartition, but allows ranking of variables based on genome-wide contribution to variance and analysis of gene-level results.

Modelling error in gene expression measurements

Uncertainty in the measurement of gene expression can be modeled with precision weights and tests of differentially expression using limma::voom() model this uncertainty directly with a heteroskedastic linear regression (Law et al. 2014). variancePartition can use these precision weights in a heteroskedastic linear mixed model implemented in lme4 (Bates et al. 2015). These precision weights are used seamlessly by calling fitVarPartModel() or fitExtractVarPartModel() on the output of limma::voom(). Otherwise the user can specify the weights with the weightsMatrix parameter.