BACKGROUND
high-throughput genomics assays have revolutionized our understanding of the molecular etiology of human disease, molecular biology of cell lineage and genetic regulation of gene expression.
transcriptome profiling in particular has been widely applied to detect variation in transcript levels attributable to differences in disease state, cell type or regulatory genetics.
as transcriptome profiling studies have expanded in size and scope, they have grown increasingly complex and consider multiple sources of biological and technical variation.
recent studies have simultaneously considered multiple dimensions of variation to understand the impact of cell type  <cit> , tissue type  <cit> , brain region  <cit> , experimental stimuli  <cit> , time duration following stimulus  <cit>  or ancestry  <cit>  on the genetic regulation of gene expression.
more studies are including a disease axis, for example to characterize the role of regulatory genetics on late onset alzheimer’s disease in multiple brain regions  <cit> .
the fundamental challenge in the analysis of complex datasets is to quantify and interpret the contribution of multiple sources of variation.
indeed the most pressing questions concern the relationship between these sources of variation.
how does cell or tissue type affect the genetic regulation of gene expression, and does it vary by ancestry  <cit> ?
what is the relative contribution of experimental stimulus versus regulatory genetics to variation in gene expression  <cit> ?
is technical variability of rna-seq low enough to study regulatory genetics and disease biology, and what are the major drivers of this technical variability  <cit> ?
a rich understanding of complex datasets requires answering these questions with both a genome-wide summary and gene-level resolution.
standard computational workflows employ principal components analysis  <cit>  and hierarchical clustering  <cit>  to summarize genome-wide expression patterns, and differential expression  to perform gene-level analyses.
recently, statistical methods that decompose variation in gene expression into the variance attributable to multiple variables in the experimental design have yielded valuable insight into the biological and technical components driving expression variation .
moreover, linear mixed models have been widely used in the analysis and interpretation of genome-wide association studies .
the linear mixed model is uniquely suited to interpreting drivers of variation in complex gene expression studies.
yet the lack of a convenient workflow and scalable implementation for analysis and visualization have prevented wider application of this rich statistical framework.
applying this analysis framework to gene expression data currently requires particular expertise in regression modeling, computational statistics, the r programming language and data visualization.
even then, the time required to implement the analysis is often prohibitive.
as gene expression datasets become more complex, the analysis and interpretation of the data is becoming the rate-limiting step.
we have developed the variancepartition software and workflow to facilitate rapid analysis and improve interpretation of complex gene expression datasets.
the software and workflow enables any analyst to perform a sophisticated analysis and visualize the results in hours using a few lines of r code.
variancepartition leverages the power of the linear mixed model  to jointly quantify the contribution of multiple sources of variation in high-throughput genomics studies.
in applications to transcriptome profiling, variancepartition fits a linear mixed model for each gene and partitions the total variance into the fraction attributable to each aspect of the study design, plus the residual variation.
because it is built on the first principles of the linear mixed model, variancepartition has well characterized theoretical properties  and accurately estimates the variance fractions even for complex experimental designs where the standard anova method is either inaccurate or not applicable.
moreover, variancepartition gives strong interpretations about the drivers of expression variation, and we demonstrate that these findings are reproducible across multiple datasets.
here we apply variancepartition to four well-characterized gene expression studies to demonstrate how the workflow facilitates interpretation of drivers of expression variation in complex study designs with multiple dimensions of variation.
we illustrate how variancepartition enables rapid interpretation of the drivers of expression variation in these complex datasets.
implementation
overview of the software
the variancepartition r package implements a computational workflow  that is complementary to standard analyses and provides particular insight into datasets with multiple dimensions of variation.
variancepartition provides a user-friendly, parallelized interface for genome-wide analysis and publication quality visualizations to examine the results.
because the variance fractions are simple to describe and interpret, variancepartition can give particular insight into how each dimension of variation contributes to transcriptional variability.
a typical variancepartition analysis comprises: 1) fitting a linear mixed model to quantify the contribution of each dimension of variation to each gene, 2) visualizing the results, and 3) integrating additional data about each gene to interpret drivers of this variation.
the variancepartition workflow requires only a few lines of r code for pre-processing, analysis and visualization and this enables rapid interpretation of complex datasets.
fig.
<dig> analysis workflow of gene expression data and meta-data.
standard analysis consists of interpreting gene expression data with respect to variables in the metadata using genome-wide analysis such as a principal components analysis and b hierarchical clustering, and gene-level analysis such as c differential expression.
the variancepartition workflow uses a rich statistical framework in the form of a linear mixed model and produces gene-level results and a genome-wide summary to simultaneously interpret gene expression data in the context of multiple variables in the metadata.
the workflow produces d gene-level results quantifying the contribution of each metadata variable to the variation in expression of each gene, and e a violin plot to summarize the genome-wide trend and rank the total contribution of each variable.
f the gene-level results can be used to identify genes that show high expression variation across individuals  or tissue .
furthermore, variancepartition facilitates examination of specific genes, and integrating external data enables further interpretation of the drivers of expression variation
the variancepartition software is implemented in r and is optimized for genome-wide analysis of large-scale transcriptome profiling datasets.
variancepartition uses the packages lme <dig>  <cit>  foreach  <cit> , iterators  <cit>  and doparallel  <cit>  to efficiently fit a linear mixed model for each gene in parallel on a multicore machine with a small memory footprint.
the precision weights from limma/voom  <cit>  are seamlessly incorporated into the analysis workflow.
built-in publication quality visualizations are implemented in ggplot <dig>  <cit> .
the variancepartition software including extensive documentation is available from http://bioconductor.org/packages/variancepartition and is compatible with bioconductor ≥v <dig>  for r ≥v <dig> .
linear mixed model framework
variancepartition summarizes the contribution of each variable in terms of the fraction of variation explained .
while the concept of fve is widely applied to univariate regression by reporting the r
<dig> value from a simple linear model, variancepartition extends fve to applications with complex study designs with multiple variables of interest.
the linear mixed model framework of variancepartition allows multiple dimensions of variation to be considered jointly in a single model and accommodates discrete variables with a large number of categories.
this analysis has a similar motivation as the standard anova method.
yet the linear mixed model framework has several statistical and practical advantages that make it more accurate and generally applicable to complex study designs with multiple dimensions of variation .
each gene is analyzed separately using the linear mixed model
<dig> y=∑jxjβj+∑kzkαk+ε
<dig> αk∼n
<dig> ε∼n
where y is the expression of a single gene across all samples, x
j is the matrix of j
th fixed effect with coefficients β
j,z
k is the matrix corresponding to the k
th random effect with coefficients α
k drawn from a normal distribution with variance σαk <dig>  the noise term, ε, is drawn from a normal distribution with variance σε <dig>  all parameters are estimated with maximum likelihood  <cit>  as simulations under a range of experimental designs indicate that this approach gives the most accurate fve estimates .
variance terms for the fixed effects are computed using the post hoc calculation
<dig> σ^βj2=var.
the total variance is
<dig> σ^total2=∑jσ^βj2+∑kσ^αk2+σ^ε <dig>
so that the fraction of variance explained by the j
th fixed effect is
<dig> σ^βj2/σ^total <dig>
by the k
th random effect is
<dig> σ^αk2/σ^total <dig>
and the residual variance is
<dig> σ^ε2/σ^total <dig>
in the standard application of variancepartition, these fractions sum to  <dig> and are always positive by definition.
moreover, the fraction of variation is also interpretable in terms of intra-class correlation, a metric used to assess biological and technical reproducibility  <cit> .
each gene is processed separately so that only visualization and reporting of genome-wide summary statistics use the results from multiple genes.
parameter estimation
the formulation of the linear mixed model is very general and includes as special cases models where only fixed effects or only random effects are used.
when only fixed effects are used, this model corresponds to a fixed effects analysis of variance  where parameters can be estimated with ordinary least squares.
when random effects are specified, the variance terms can be estimated with maximum likelihood or restricted maximum likelihood   <cit> .
since reml does not directly estimate parameters for fixed effects, these coefficients are estimated after the fact by plugging in estimates for the variance components  <cit> .
we focus on the most general case  that includes both fixed and random effects.
in this case parameters can be estimated with maximum likelihood.
maximum likelihood estimates are used exclusively in the main text and are the default in the variancepartition software when random effects are specified because this method performs best in simulations.
relationship to existing methods
the fixed effects anova model has been widely applied for decades to decompose variance into multiple components of variation  <cit> .
yet this approach is often inadequate to address the questions that are posed by complex gene expression datasets.
the linear mixed model used by variancepartition has three distinct advantages compared to anova.
first, by placing a gaussian prior on variables modeled as random effects, the linear mixed model more accurately estimates the fraction of variance explained.
even as the number of categories in a discrete variable increases, the linear mixed model still produces accurate estimates because the prior shrinks the estimate for each category towards the zero.
conversely, the fixed effects anova is fit with a linear regression model using ordinary least squares.
this method is known to suffer from overfitting and over-estimates the variance fractions for variables with many categories.
these properties are well established  <cit>  and are consistent with our simulation study .
second, the linear mixed model can decompose variance into multiple components in situations where the fixed effects anova cannot be applied because the design matrix is degenerate .
this situation is very common for the types of question of relevant to complex gene expression studies.
for example, sex and ancestry are invariant properties of an individual, so jointly analyzing variation across these  <dig> dimensions of variation involves a degenerate design matrix.
in cases like these, the linear mixed model can accurately estimate the desired variance fractions , while anova will fail to estimate any of these values because the parameters are not identifiable.
thus anova is inadequate for the type of analysis we performed here with variancepartition using linear mixed model.
finally, the linear mixed model can quantify how variation attributable to one aspect of the study design depends on another, such as the case of cross-individual expression variation depending on tissue/cell type.
anova does not have this capability.
interpretation of percent variance explained
the percent variance explained can be interpreted as the intra-class correlation .
consider the simplest example of the i
th sample from the k
th individual
<dig> yi,k=μ+αk+εi,k
<dig> αk∼n
<dig> εi,k∼n
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:
<dig> icc=cor
<dig> =cor
<dig> =covvarvar
<dig> =covσα2+σε <dig>
<dig> =σα2σα2+σε <dig>
the correlation between samples from different individuals is:
<dig> =cor
<dig> =cor
<dig> =covσα2+σε <dig>
<dig> =0σα2+σε <dig>
<dig> = <dig>
this interpretation in terms of fraction of variation explained  naturally generalizes to multiple variance components  <cit> .
see additional file  <dig> for more details.
variation across individual within 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  <cit> .
this analysis examines the expression variation across individuals within multiple cell types, or another subset of the data.
a given sample is only from one cell type, so this analysis asks a question about a subset of the data.
the data is implicitly divided into subsets based on cell type and variation explained by individual is evaluated within each subset.
this subsetting means that the variance fractions no longer sum to  <dig>  but the model still allows ranking of dimensions of variation based on genome-wide contribution to variance and enables analysis of gene-level results.
see the additional file  <dig> for more details.
modeling measurement error in rna-seq data
uncertainty in the measurement of rna-seq data can be modeled with observation-level precision weights that model the relationship between expression magnitude and sampling variance  <cit> .
variancepartition naturally incorporates these precision weights to create a heteroskedastic linear mixed model  <cit>  that can explicitly account from the measurement uncertainty due to the finite count nature of rna-seq data.
let the precision w
i denote the inverse of the variance of the observation y
i for the i
th observation.
the precisions can be used to re-weight the samples in a regression to account for the variation in the uncertainty about each observation.
weighting by the precision upweights samples with low measurement error and down weights samples with high measurement error.
denoting the vector of precision weights for a single gene across all samples as w, the model is fit by weighting the residual variance from equation
<dig> ε∼nσε2).
these weights are estimated using limma/voom  <cit>  in a preprocessing step and are then incorporated into the variancepartition analysis.
RESULTS
analysis of geuvadis rna-seq dataset
consider  <dig> rna-seq experiments from the geuvadis study  <cit>  of lymphoblastoid cell lines from  <dig> individuals of  <dig> ancestries and  <dig> sexes sequenced across  <dig> laboratories.
for a single gene, the total variance can be partitioned into the contributions of these components of variation plus residual variance:
<dig> σtotal2=σindividual2+σlab2+σancestry2+σsex2+σε <dig>
the contribution of each driver of variation can be interpreted based on the fraction of total variation it explains.
thus the fraction of variance due to variation across individuals is
<dig> σindividual2/σtotal <dig>
and the fractions from all components of variation sum to  <dig>
applying variancepartition to the geuvadis  <cit>  dataset illustrates how the method can decouple biological and technical variation, and further decompose biological variation into multiple components.
expression variation across individuals, ancestries and sexes is biological, variation across the labs where the samples were sequenced comprise technical artifacts, while the residual variation remains uncharacterized.
results from representative genes illustrate how variancepartition identifies genes where the majority of variation in expression is explained by a single variable such as individual or sex, while variation in other genes is driven by multiple variables .
since the variance fractions sum to  <dig> for each gene, it is simple to compare results across genes and across sources of variation.
visualizing these results genome-wide illustrates that variation across individuals is the major source of expression variation and explains a median of  <dig> % of variance genome-wide .
the median variance explained by laboratory , ancestry  and sex  is substantially smaller.
we note that the variance explained by individual increases to  <dig> % when ancestry is removed from the analysis since ancestry is a biological property of each individual .
fig.
<dig> analysis of geuvadis dataset identifies drivers of expression variation.
a total variance for each gene is partitioned into the fraction attributable to each dimension of variation in the study.
b violin and box plots of percent variation in gene expression explained by each variable.
three representative genes and their major sources of variation are indicated.
c boxplot of uty expression stratified by sex.
d boxplot of ccdc85b expression stratified by lab.
inset shows scatter plot of percent gc content versus percent variance explained by lab.
red line indicates linear regression line with coefficient of determination and p-value shown.
e boxplot of znf <dig> expression stratified by individual for a subset of individuals with at least  <dig> technical replicate.
inset illustrates a cis-eqtl for znf <dig> where expression is stratified by genotype at rs <dig>  f probability of each gene having a cis-eqtl plotted against the percent variance explained by individual.
dashed lines indicate the genome-wide average probability , and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual.
points indicate a sliding window average of the probability of genes in each window having a cis-eqtl.
window size is  <dig> genes with an overlap of  <dig> genes between windows.
the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis
yet particular genes show substantial deviation from the genome-wide trend.
this is particularly noticeable for sex, where of the  <dig> genes for which sex explains more that 10% variance  <dig> are on the x or y chromosomes.
for example, variation across sex explains 98% variance in uty on the y chromosome .
while differential expression measures the differences in mean expression between the sexes, variancepartition measures the variance within and between each sex.
this analysis indicates that variation across sexes explains very little variation genome-wide, but has a strong effect on a small number of genes.
integrating additional data with gene-level results from variancepartition can give a clear interpretation of the drivers of variation.
for example,  <dig> % of variation in ccdc85b is explained by variation across laboratory.
this gene has a very high gc content of  <dig> % and is consistent with the genome-wide pattern where the degree of variation across laboratories is positively correlated with gc content .
while technical variation in rna-seq is known to depend on gc content  <cit> , variancepartition gives a clear illustration of how the effect of technical artifacts varies substantially across genes.
moreover, this analysis can be used to identify other correlates underlying technical issues in expression variation.
in addition, variancepartition gives a strong interpretation to genes whose expression varies across individuals by relating the gene-level results to cis-regulatory variation.
for example, the fact that  <dig> % of variation in znf <dig> is explained by individual suggests that this variation is driven by genetics, and, in fact, znf <dig> has a cis-eqtl .
this observation is also seen genome-wide, as genes with a greater fraction of variation across individuals have a significantly higher probability of having a cis-eqtl detected in this study .
this analysis explicitly demonstrates how expression variation across individuals is driven by cis-regulatory variation.
analysis of seqc rna-seq dataset
the sequencing quality control  project  <cit>  evaluated the technical reproducibility of rna-seq data by sequencing the same  <dig> rna samples at  <dig> laboratories, using  <dig> total library constructions and up to  <dig> lanes on each of  <dig> illumina hiseq  <dig> flowcells for a total of  <dig> rna-seq experiments.
the goal of the study was to determine the degree to which these technical factors explain variation in gene expression measurements.
this complex dataset has multiple levels of variation and variancepartition provides a rigorous statistical framework to quantify and interpret these sources of variation in a single analysis.
as expected, variation across the  <dig> rna samples is the major axis of variation, explaining a median of  <dig> % of variation in expression .
but the real interest is in the sources of technical variability.
the fact that the technical variables laboratory , library , flowcell , and lane  explain a small fraction of the total variation indicate that these rna-seq experiments were highly reproducible genome-wide.
interpreting these values in terms of the intra-class correlation indicates that two experiments from the sample rna sample but which differ in all other aspects of the study design are highly correlated .
conversely, two experiments from the same lane, but different rna samples, etc, show negligible correlation as is expected when technical variation is low.
fig.
<dig> analysis of sequencing quality control  dataset decouples sources of technical variation.
a violin and box plots of percent variation in gene expression explained by each variable.
b boxplot of percent variance explained by rna sample for human genes and external rna controls consortium  spike-in controls.
p-value is from one-sided mann-whitney test.
c scatter plot of percent gc content and percent variance explained by laboratory.
red line indicates linear regression line with regression coefficient, coefficient of determination and p-value shown
analysis and visualization with variancepartition succinctly illustrates that while variation across laboratories and library constructions is not negligible, it is small compared with the magnitude of biological variation for the large majority of genes.
moreover, variation across flowcells and lanes is very small in this dataset.
thus variancepartition illustrates that rna-seq data is highly reproducible genome-wide with a small subset of genes showing large technical artifacts.
however, there are notable deviations from these genome-wide trends.
first, there are a set of transcripts that show little variation between the  <dig> rna samples and, in fact, these correspond to spike-in synthetic rna added to each sample at a standardized concentration to act as controls having equal abundance in all experiments  <cit> .
as expected, spike-in transcripts show significantly less variation across the  <dig> rna samples than human genes .
second, although technical effects are low for most genes, a small number of genes show high variation across laboratories and library constructions.
in fact, the fraction of variation across laboratories correlates with the gc content of each gene , and recapitulates the known role of gc content with reproducibility of rna-seq data .
analysis of immvar microarray dataset
the immune variation  project assayed gene expression in cd <dig> +cd <dig> − monocytes and cd <dig> + t-cells on the affymetrix human gene  <dig>  st array platform in order to characterize the role of cell type in genetic regulation of gene expression  <cit> .
analysis of  <dig> individuals with data from both cell types reveals that multiple variables contribute to expression variation in this dataset .
since variancepartition reports the contribution of each variable while simultaneously correcting for all other values, it is apparent that the variation across cell types is the strongest biological driver of variation  followed by variation across individuals .
although cell type has a smaller median effect than batch, it is notable that cell type explains >50% of the variation for  <dig>  genes.
the observation that batch and cell type are the strongest drivers of variation is largely consistent with results from principal components analysis  .
we note that the relationship between variancepartition and pca depends on both the fraction of expression variation explained by a particular variable across all genes as well as the dimension of the variable.
while variation across the  <dig> cell types explains less expression variation than variation across the  <dig> batches, the first principal component separates samples by cell type because this variable spans a lower-dimensional space.
fig.
<dig> analysis of immvar dataset interprets multiple dimensions of expression variation.
a violin and box plots of percent variation in gene expression explained by each variable.
b principal components analysis of gene expression with experiments colored by batch.
c total variance for each gene is partitioned into the fraction attributable to each dimension of variation in the study design.
d expression of uty stratified by sex.
e expression of tlr <dig> stratified by cell type.
f expression of gstm <dig> stratified by individual.
g scatter plot of percent gc content and percent variance explained by batch.
red line indicates linear regression line with regression coefficient, coefficient of determination and p-value shown.
h results from variancepartition analysis allowing the contribution of individual to vary in each cell type.
i probability of each gene having a cis-eqtl plotted against the percent variance explained by individual within each cell type.
dashed lines indicate the genome-wide average probability, and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual within each cell type.
points indicate a sliding window average of the probability of genes in each window having a cis-eqtl.
window size is  <dig> genes with an overlap of  <dig> genes between windows.
the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis
meanwhile, sex drives expression variation in a small number of genes, while the age of each individual has a negligible effect.
we note that despite the large batch effect observed in this dataset, the biological variation across cell type, individual and sex are still large enough to make meaningful conclusions about cell-specific regulatory genetics when this technical effect is accounted for  <cit> .
moreover, variancepartition identifies genes that vary along different aspects of the study design , and visualization of a subset of these genes illustrates the strong expression differences when stratified by sex, cell type and individual .
variancepartition enables further interpretation of the batch effect because it gives results at a gene-level resolution.
the samples were processed in  <dig> technical batches and this axis of variation explains a median of  <dig> % of total variation, indicating a large technical effect.
consistent with other analyses, the fraction of variation explained by batch at the gene-level is positively correlated with gc content .
by leveraging the flexibility of the linear mixed model, variancepartition can quantify the variation across individuals within each cell type.
since the variance is analyzed within multiple subsets of the data and each sample is only in a single subset, the total variation explained no longer sums to  <dig> as it does for standard application of variancepartition.
yet the results allow ranking of dimensions of variation based on genome-wide contribution to variance and enables analysis of gene-level results .
this analysis uses the fact that  <dig> individuals within monocytes have at least  <dig> technical replicate, while  <dig> individuals within t-cells have at least  <dig> technical replicate.
the variation across individuals within t-cells  and monocytes  is substantially larger than when the two cell types were combined .
the fact that the contribution of individual varies between cell types is consistent with cell-specific regulatory genetics  <cit> .
finally, the fraction of variation explained by individual within each cell type at the gene-level is directly related to the probability of each gene having cis-eqtl within the corresponding cell type .
analysis of gtex rna-seq dataset
application of variancepartition to post mortem rna-seq data of multiple tissues tissues from the gtex consortium  <cit>  decouples the influence of multiple biological and technical drivers of expression variation.
we analyzed  <dig> experiments from  <dig> individuals in  <dig> tissues  in order to restrict the analysis to tissues with rna-seq data for most individuals .
variation across tissues is the major source of variation  while the technical variables expression batch , ischemic time , rna isolation batch , and rin  have a moderate effect on expression variation genome-wide .
variation across expression batches is correlated with gc content but to a lesser degree that other datasets .
cumulatively, these technical variables explain only  <dig> % of the total expression variation.
concerns about reliability of rna-seq data from post mortem samples has been raised due to the potential effects of rna degradation following cell death  <cit> .
variancepartition analysis indicates that variation in ischemic time has as relatively small effect genome-wide and the fraction of variance it explains is comparable to technical effects, yet the effect varies substantially across genes.
fig.
<dig> analysis of gtex dataset identifies drivers of expression variation at multiple levels.
a violin and box plots of percent variation in gene expression explained by each variable.
b results from variancepartition analysis allowing the contribution of individual to vary in each tissue.
c probability of each gene having a cis-eqtl plotted against the percent variance explained by individual within each tissue.
dashed lines indicate the genome-wide average probability, and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual within each tissue.
the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis.
d fraction of variation in glmp explained by each source of variation.
e glmp has a cis-eqtl active in blood but not skin
the flexibility of the linear mixed model framework allows variancepartition to analyze cross-individual variation within each tissue.
we note again that since the variance is analyzed within multiple subsets of the data, the total variation explained no longer sums to  <dig> here.
while variation across individuals explains only a median of  <dig> % of variation when all tissues types are considered together, there is substantial variation across individuals within each tissue separately .
cross-individual variation is highest in blood , while skin , blood vessel , and adipose tissue  exhibit lower cross-individual variation.
the fraction of variation explained by individual within each tissue is directly related to the probability of each gene having a cis-eqtl within the corresponding tissue .
this association is not as strong as in other datasets likely due to the smaller number of individuals and to the relatively small fraction of variation across individuals in adipose tissue.
at the gene-level, variancepartition can prioritize genes based on multiple criteria.
for examples, glmp exhibits higher variation across individuals within blood but low variation in skin .
this is consistent of a tissue-specific regulatory variation, and, in fact, the cis-eqtl rs <dig> influences gene expression in blood but not in skin .
discussion
as the scope of gene expression studies continues to expand, the need to quantify and interpret multiple drivers of expression variation is becoming essential.
here we present variancepartition, a publicly available software package that leverages the power of the linear mixed model to quantify the contribution of multiple sources of variation in complex gene expression datasets.
for each gene, this analysis partitions the total expression variance into the fraction attributable to each aspect of the study design.
a variancepartition analysis gives a genome-wide summary of the drivers of variation, but also produces gene-level results to identify genes that deviate from the genome-wide trend.
the fraction of expression variation is easily interpretable across genes, drivers of variation and datasets.
thus variancepartition produces a more detailed and quantitative genome-wide overview than the standard principal components analysis   <cit>  and hierarchical clustering   <cit>  approaches.
pca and hc focus on the major axis of variation, and they overlook the secondary drivers of variation that can be well characterized with variancepartition.
moreover, the gene-level results from variancepartition indicate genes that deviate from the genome-wide trend and integration with additional data can enable a further interpretation.
while pca and hc do not give gene-level results, differential expression  analysis reports gene-level fold change and corresponding p-value for each aspect of the study design.
yet de analysis does not produce a clear genome-wide summary, and the fold change and p-values are not easily comparable across multiple drivers of variation.
analysis of publicly available gene expression studies demonstrate that variancepartition recovers striking patterns of biological and technical variation that are reproducible across multiple datasets.
at a genome-wide level, expression variation across individuals and cell types is large enough to overcome the technical variation of transcriptome profiling.
yet at the gene-level there is substantial deviation from the genome-wide trend due to a range of biological and technical factors.
by quantifying the variance attributable to each aspect of the study design, variancepartition facilitates the interpretation of these gene-level effects in the context of additional information.
we demonstrate reproducible findings that cross-individual variation is driven by cis-eqtl’s and technical variation across laboratories associated with gc content.
moreover, variation across individuals and the relationship to cis-eqtl’s depend on the cell or tissue type.
CONCLUSIONS
the variancepartition workflow and implementation makes the rich linear mixed model framework easily applicable for interpreting drivers of variation in complex gene expression data.
variancepartition provides a general statistical and visualization framework for studying drivers of variation in rna-seq datasets in many types of high-throughput genomic assays including rna-seq , protein quantification, metabolite quantification, metagenomic assays, methylation arrays and epigenomic sequencing assays.
although we have focused here on large-scale studies, variancepartition analysis has given valuable insight into rna-seq datasets with as few as  <dig> samples.
the variancepartition software is an open source r package and is freely available on bioconductor.
the software can easily be applied to rna-seq quantifications from featurecounts  <cit> , htseq  <cit> , kallisto  <cit> , sailfish  <cit> , salmon  <cit> , rsem  <cit>  and cufflinks  <cit>  which have been processed in r with limma/voom  <cit> , deseq <dig>  <cit> , tximport  <cit>  and ballgown  <cit> .
the software provides a user-friendly interface for analysis and visualization with extensive documentation, and will enable routine application to a range of genomics datasets.
availability and requirements
project name: variancepartition
project home page:
http://bioconductor.org/packages/variancepartition
operating systems: linux, mac os x, windows
programming language: r ≥ v <dig>
other requirements: bioconductor ≥ v <dig>
license: gpl-2
additional file
additional file  <dig> supplementary note.
additional details describing statistical methods and software.
includes results from simulation study and additional figures from data analysis.
abbreviations
anovaanalysis of variance
dedifferential expression
eqtlexpression quantitative trait locus
fvefraction of variance explained
hchierarchical clustering
iccintra-class correlation
iqrinter-quartile range
pcaprincipal components analysis
remlrestricted maximum likelihood

 high-throughput genomics assays have revolutionized our understanding of the molecular etiology of human disease, molecular biology of cell lineage and genetic regulation of gene expression. transcriptome profiling in particular has been widely applied to detect variation in transcript levels attributable to differences in disease state, cell type or regulatory genetics. as transcriptome profiling studies have expanded in size and scope, they have grown increasingly complex and consider multiple sources of biological and technical variation. recent studies have simultaneously considered multiple dimensions of variation to understand the impact of cell type, tissue type, brain region, experimental stimuli, time duration following stimulus or ancestry on the genetic regulation of gene expression. more studies are including a disease axis, for example to characterize the role of regulatory genetics on late onset alzheimer’s disease in multiple brain regions. the fundamental challenge in the analysis of complex datasets is to quantify and interpret the contribution of multiple sources of variation. indeed the most pressing questions concern the relationship between these sources of variation. how does cell or tissue type affect the genetic regulation of gene expression, and does it vary by ancestry? what is the relative contribution of experimental stimulus versus regulatory genetics to variation in gene expression? is technical variability of rna-seq low enough to study regulatory genetics and disease biology, and what are the major drivers of this technical variability? a rich understanding of complex datasets requires answering these questions with both a genome-wide summary and gene-level resolution. standard computational workflows employ principal components analysis and hierarchical clustering to summarize genome-wide expression patterns, and differential expression to perform gene-level analyses. recently, statistical methods that decompose variation in gene expression into the variance attributable to multiple variables in the experimental design have yielded valuable insight into the biological and technical components driving expression variation. moreover, linear mixed models have been widely used in the analysis and interpretation of genome-wide association studies. the linear mixed model is uniquely suited to interpreting drivers of variation in complex gene expression studies. yet the lack of a convenient workflow and scalable implementation for analysis and visualization have prevented wider application of this rich statistical framework. applying this analysis framework to gene expression data currently requires particular expertise in regression modeling, computational statistics, the r programming language and data visualization. even then, the time required to implement the analysis is often prohibitive. as gene expression datasets become more complex, the analysis and interpretation of the data is becoming the rate-limiting step. we have developed the variancepartition software and workflow to facilitate rapid analysis and improve interpretation of complex gene expression datasets. the software and workflow enables any analyst to perform a sophisticated analysis and visualize the results in hours using a few lines of r code. variancepartition leverages the power of the linear mixed model to jointly quantify the contribution of multiple sources of variation in high-throughput genomics studies. in applications to transcriptome profiling, variancepartition fits a linear mixed model for each gene and partitions the total variance into the fraction attributable to each aspect of the study design, plus the residual variation. because it is built on the first principles of the linear mixed model, variancepartition has well characterized theoretical properties and accurately estimates the variance fractions even for complex experimental designs where the standard anova method is either inaccurate or not applicable. moreover, variancepartition gives strong interpretations about the drivers of expression variation, and we demonstrate that these findings are reproducible across multiple datasets. here we apply variancepartition to four well-characterized gene expression studies to demonstrate how the workflow facilitates interpretation of drivers of expression variation in complex study designs with multiple dimensions of variation. we illustrate how variancepartition enables rapid interpretation of the drivers of expression variation in these complex datasets. implementation overview of the software the variancepartition r package implements a computational workflow that is complementary to standard analyses and provides particular insight into datasets with multiple dimensions of variation. variancepartition provides a user-friendly, parallelized interface for genome-wide analysis and publication quality visualizations to examine the results. because the variance fractions are simple to describe and interpret, variancepartition can give particular insight into how each dimension of variation contributes to transcriptional variability. a typical variancepartition analysis comprises: 1) fitting a linear mixed model to quantify the contribution of each dimension of variation to each gene, 2) visualizing the results, and 3) integrating additional data about each gene to interpret drivers of this variation. the variancepartition workflow requires only a few lines of r code for pre-processing, analysis and visualization and this enables rapid interpretation of complex datasets. fig. analysis workflow of gene expression data and meta-data. standard analysis consists of interpreting gene expression data with respect to variables in the metadata using genome-wide analysis such as a principal components analysis and b hierarchical clustering, and gene-level analysis such as c differential expression. the variancepartition workflow uses a rich statistical framework in the form of a linear mixed model and produces gene-level results and a genome-wide summary to simultaneously interpret gene expression data in the context of multiple variables in the metadata. the workflow produces d gene-level results quantifying the contribution of each metadata variable to the variation in expression of each gene, and e a violin plot to summarize the genome-wide trend and rank the total contribution of each variable. f the gene-level results can be used to identify genes that show high expression variation across individuals or tissue. furthermore, variancepartition facilitates examination of specific genes, and integrating external data enables further interpretation of the drivers of expression variation the variancepartition software is implemented in r and is optimized for genome-wide analysis of large-scale transcriptome profiling datasets. variancepartition uses the packages lme foreach, iterators and doparallel to efficiently fit a linear mixed model for each gene in parallel on a multicore machine with a small memory footprint. the precision weights from limma/voom are seamlessly incorporated into the analysis workflow. built-in publication quality visualizations are implemented in ggplot. the variancepartition software including extensive documentation is available from http://bioconductor.org/packages/variancepartition and is compatible with bioconductor ≥v for r ≥v. linear mixed model framework variancepartition summarizes the contribution of each variable in terms of the fraction of variation explained. while the concept of fve is widely applied to univariate regression by reporting the r value from a simple linear model, variancepartition extends fve to applications with complex study designs with multiple variables of interest. the linear mixed model framework of variancepartition allows multiple dimensions of variation to be considered jointly in a single model and accommodates discrete variables with a large number of categories. this analysis has a similar motivation as the standard anova method. yet the linear mixed model framework has several statistical and practical advantages that make it more accurate and generally applicable to complex study designs with multiple dimensions of variation. each gene is analyzed separately using the linear mixed model y=∑jxjβj+∑kzkαk+ε αk∼n ε∼n where y is the expression of a single gene across all samples, x j is the matrix of j th fixed effect with coefficients β j,z k is the matrix corresponding to the k th random effect with coefficients α k drawn from a normal distribution with variance σαk the noise term, ε, is drawn from a normal distribution with variance σε all parameters are estimated with maximum likelihood as simulations under a range of experimental designs indicate that this approach gives the most accurate fve estimates. variance terms for the fixed effects are computed using the post hoc calculation σ^βj2=var. the total variance is σ^total2=∑jσ^βj2+∑kσ^αk2+σ^ε so that the fraction of variance explained by the j th fixed effect is σ^βj2/σ^total by the k th random effect is σ^αk2/σ^total and the residual variance is σ^ε2/σ^total in the standard application of variancepartition, these fractions sum to and are always positive by definition. moreover, the fraction of variation is also interpretable in terms of intra-class correlation, a metric used to assess biological and technical reproducibility. each gene is processed separately so that only visualization and reporting of genome-wide summary statistics use the results from multiple genes. parameter estimation the formulation of the linear mixed model is very general and includes as special cases models where only fixed effects or only random effects are used. when only fixed effects are used, this model corresponds to a fixed effects analysis of variance where parameters can be estimated with ordinary least squares. when random effects are specified, the variance terms can be estimated with maximum likelihood or restricted maximum likelihood. since reml does not directly estimate parameters for fixed effects, these coefficients are estimated after the fact by plugging in estimates for the variance components. we focus on the most general case that includes both fixed and random effects. in this case parameters can be estimated with maximum likelihood. maximum likelihood estimates are used exclusively in the main text and are the default in the variancepartition software when random effects are specified because this method performs best in simulations. relationship to existing methods the fixed effects anova model has been widely applied for decades to decompose variance into multiple components of variation. yet this approach is often inadequate to address the questions that are posed by complex gene expression datasets. the linear mixed model used by variancepartition has three distinct advantages compared to anova. first, by placing a gaussian prior on variables modeled as random effects, the linear mixed model more accurately estimates the fraction of variance explained. even as the number of categories in a discrete variable increases, the linear mixed model still produces accurate estimates because the prior shrinks the estimate for each category towards the zero. conversely, the fixed effects anova is fit with a linear regression model using ordinary least squares. this method is known to suffer from overfitting and over-estimates the variance fractions for variables with many categories. these properties are well established and are consistent with our simulation study. second, the linear mixed model can decompose variance into multiple components in situations where the fixed effects anova cannot be applied because the design matrix is degenerate. this situation is very common for the types of question of relevant to complex gene expression studies. for example, sex and ancestry are invariant properties of an individual, so jointly analyzing variation across these dimensions of variation involves a degenerate design matrix. in cases like these, the linear mixed model can accurately estimate the desired variance fractions, while anova will fail to estimate any of these values because the parameters are not identifiable. thus anova is inadequate for the type of analysis we performed here with variancepartition using linear mixed model. finally, the linear mixed model can quantify how variation attributable to one aspect of the study design depends on another, such as the case of cross-individual expression variation depending on tissue/cell type. anova does not have this capability. interpretation of percent variance explained the percent variance explained can be interpreted as the intra-class correlation. consider the simplest example of the i th sample from the k th individual yi,k=μ+αk+εi,k αk∼n εi,k∼n 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: icc=cor =cor =covvarvar =covσα2+σε =σα2σα2+σε the correlation between samples from different individuals is: =cor =cor =covσα2+σε =0σα2+σε = this interpretation in terms of fraction of variation explained naturally generalizes to multiple variance components. see additional file for more details. variation across individual within 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. this analysis examines the expression variation across individuals within multiple cell types, or another subset of the data. a given sample is only from one cell type, so this analysis asks a question about a subset of the data. the data is implicitly divided into subsets based on cell type and variation explained by individual is evaluated within each subset. this subsetting means that the variance fractions no longer sum to but the model still allows ranking of dimensions of variation based on genome-wide contribution to variance and enables analysis of gene-level results. see the additional file for more details. modeling measurement error in rna-seq data uncertainty in the measurement of rna-seq data can be modeled with observation-level precision weights that model the relationship between expression magnitude and sampling variance. variancepartition naturally incorporates these precision weights to create a heteroskedastic linear mixed model that can explicitly account from the measurement uncertainty due to the finite count nature of rna-seq data. let the precision w i denote the inverse of the variance of the observation y i for the i th observation. the precisions can be used to re-weight the samples in a regression to account for the variation in the uncertainty about each observation. weighting by the precision upweights samples with low measurement error and down weights samples with high measurement error. denoting the vector of precision weights for a single gene across all samples as w, the model is fit by weighting the residual variance from equation ε∼nσε2). these weights are estimated using limma/voom in a preprocessing step and are then incorporated into the variancepartition analysis.  analysis of geuvadis rna-seq dataset consider rna-seq experiments from the geuvadis study of lymphoblastoid cell lines from individuals of ancestries and sexes sequenced across laboratories. for a single gene, the total variance can be partitioned into the contributions of these components of variation plus residual variance: σtotal2=σindividual2+σlab2+σancestry2+σsex2+σε the contribution of each driver of variation can be interpreted based on the fraction of total variation it explains. thus the fraction of variance due to variation across individuals is σindividual2/σtotal and the fractions from all components of variation sum to applying variancepartition to the geuvadis dataset illustrates how the method can decouple biological and technical variation, and further decompose biological variation into multiple components. expression variation across individuals, ancestries and sexes is biological, variation across the labs where the samples were sequenced comprise technical artifacts, while the residual variation remains uncharacterized. results from representative genes illustrate how variancepartition identifies genes where the majority of variation in expression is explained by a single variable such as individual or sex, while variation in other genes is driven by multiple variables. since the variance fractions sum to for each gene, it is simple to compare results across genes and across sources of variation. visualizing these results genome-wide illustrates that variation across individuals is the major source of expression variation and explains a median of % of variance genome-wide. the median variance explained by laboratory, ancestry and sex is substantially smaller. we note that the variance explained by individual increases to % when ancestry is removed from the analysis since ancestry is a biological property of each individual. fig. analysis of geuvadis dataset identifies drivers of expression variation. a total variance for each gene is partitioned into the fraction attributable to each dimension of variation in the study. b violin and box plots of percent variation in gene expression explained by each variable. three representative genes and their major sources of variation are indicated. c boxplot of uty expression stratified by sex. d boxplot of ccdc85b expression stratified by lab. inset shows scatter plot of percent gc content versus percent variance explained by lab. red line indicates linear regression line with coefficient of determination and p-value shown. e boxplot of znf expression stratified by individual for a subset of individuals with at least technical replicate. inset illustrates a cis-eqtl for znf where expression is stratified by genotype at rs f probability of each gene having a cis-eqtl plotted against the percent variance explained by individual. dashed lines indicate the genome-wide average probability, and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual. points indicate a sliding window average of the probability of genes in each window having a cis-eqtl. window size is genes with an overlap of genes between windows. the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis yet particular genes show substantial deviation from the genome-wide trend. this is particularly noticeable for sex, where of the genes for which sex explains more that 10% variance are on the x or y chromosomes. for example, variation across sex explains 98% variance in uty on the y chromosome. while differential expression measures the differences in mean expression between the sexes, variancepartition measures the variance within and between each sex. this analysis indicates that variation across sexes explains very little variation genome-wide, but has a strong effect on a small number of genes. integrating additional data with gene-level results from variancepartition can give a clear interpretation of the drivers of variation. for example, % of variation in ccdc85b is explained by variation across laboratory. this gene has a very high gc content of % and is consistent with the genome-wide pattern where the degree of variation across laboratories is positively correlated with gc content. while technical variation in rna-seq is known to depend on gc content, variancepartition gives a clear illustration of how the effect of technical artifacts varies substantially across genes. moreover, this analysis can be used to identify other correlates underlying technical issues in expression variation. in addition, variancepartition gives a strong interpretation to genes whose expression varies across individuals by relating the gene-level results to cis-regulatory variation. for example, the fact that % of variation in znf is explained by individual suggests that this variation is driven by genetics, and, in fact, znf has a cis-eqtl. this observation is also seen genome-wide, as genes with a greater fraction of variation across individuals have a significantly higher probability of having a cis-eqtl detected in this study. this analysis explicitly demonstrates how expression variation across individuals is driven by cis-regulatory variation. analysis of seqc rna-seq dataset the sequencing quality control project evaluated the technical reproducibility of rna-seq data by sequencing the same rna samples at laboratories, using total library constructions and up to lanes on each of illumina hiseq flowcells for a total of rna-seq experiments. the goal of the study was to determine the degree to which these technical factors explain variation in gene expression measurements. this complex dataset has multiple levels of variation and variancepartition provides a rigorous statistical framework to quantify and interpret these sources of variation in a single analysis. as expected, variation across the rna samples is the major axis of variation, explaining a median of % of variation in expression. but the real interest is in the sources of technical variability. the fact that the technical variables laboratory, library, flowcell, and lane explain a small fraction of the total variation indicate that these rna-seq experiments were highly reproducible genome-wide. interpreting these values in terms of the intra-class correlation indicates that two experiments from the sample rna sample but which differ in all other aspects of the study design are highly correlated. conversely, two experiments from the same lane, but different rna samples, etc, show negligible correlation as is expected when technical variation is low. fig. analysis of sequencing quality control dataset decouples sources of technical variation. a violin and box plots of percent variation in gene expression explained by each variable. b boxplot of percent variance explained by rna sample for human genes and external rna controls consortium spike-in controls. p-value is from one-sided mann-whitney test. c scatter plot of percent gc content and percent variance explained by laboratory. red line indicates linear regression line with regression coefficient, coefficient of determination and p-value shown analysis and visualization with variancepartition succinctly illustrates that while variation across laboratories and library constructions is not negligible, it is small compared with the magnitude of biological variation for the large majority of genes. moreover, variation across flowcells and lanes is very small in this dataset. thus variancepartition illustrates that rna-seq data is highly reproducible genome-wide with a small subset of genes showing large technical artifacts. however, there are notable deviations from these genome-wide trends. first, there are a set of transcripts that show little variation between the rna samples and, in fact, these correspond to spike-in synthetic rna added to each sample at a standardized concentration to act as controls having equal abundance in all experiments. as expected, spike-in transcripts show significantly less variation across the rna samples than human genes. second, although technical effects are low for most genes, a small number of genes show high variation across laboratories and library constructions. in fact, the fraction of variation across laboratories correlates with the gc content of each gene, and recapitulates the known role of gc content with reproducibility of rna-seq data. analysis of immvar microarray dataset the immune variation project assayed gene expression in cd +cd − monocytes and cd + t-cells on the affymetrix human gene st array platform in order to characterize the role of cell type in genetic regulation of gene expression. analysis of individuals with data from both cell types reveals that multiple variables contribute to expression variation in this dataset. since variancepartition reports the contribution of each variable while simultaneously correcting for all other values, it is apparent that the variation across cell types is the strongest biological driver of variation followed by variation across individuals. although cell type has a smaller median effect than batch, it is notable that cell type explains >50% of the variation for genes. the observation that batch and cell type are the strongest drivers of variation is largely consistent with results from principal components analysis. we note that the relationship between variancepartition and pca depends on both the fraction of expression variation explained by a particular variable across all genes as well as the dimension of the variable. while variation across the cell types explains less expression variation than variation across the batches, the first principal component separates samples by cell type because this variable spans a lower-dimensional space. fig. analysis of immvar dataset interprets multiple dimensions of expression variation. a violin and box plots of percent variation in gene expression explained by each variable. b principal components analysis of gene expression with experiments colored by batch. c total variance for each gene is partitioned into the fraction attributable to each dimension of variation in the study design. d expression of uty stratified by sex. e expression of tlr stratified by cell type. f expression of gstm stratified by individual. g scatter plot of percent gc content and percent variance explained by batch. red line indicates linear regression line with regression coefficient, coefficient of determination and p-value shown. h results from variancepartition analysis allowing the contribution of individual to vary in each cell type. i probability of each gene having a cis-eqtl plotted against the percent variance explained by individual within each cell type. dashed lines indicate the genome-wide average probability, and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual within each cell type. points indicate a sliding window average of the probability of genes in each window having a cis-eqtl. window size is genes with an overlap of genes between windows. the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis meanwhile, sex drives expression variation in a small number of genes, while the age of each individual has a negligible effect. we note that despite the large batch effect observed in this dataset, the biological variation across cell type, individual and sex are still large enough to make meaningful conclusions about cell-specific regulatory genetics when this technical effect is accounted for. moreover, variancepartition identifies genes that vary along different aspects of the study design, and visualization of a subset of these genes illustrates the strong expression differences when stratified by sex, cell type and individual. variancepartition enables further interpretation of the batch effect because it gives results at a gene-level resolution. the samples were processed in technical batches and this axis of variation explains a median of % of total variation, indicating a large technical effect. consistent with other analyses, the fraction of variation explained by batch at the gene-level is positively correlated with gc content. by leveraging the flexibility of the linear mixed model, variancepartition can quantify the variation across individuals within each cell type. since the variance is analyzed within multiple subsets of the data and each sample is only in a single subset, the total variation explained no longer sums to as it does for standard application of variancepartition. yet the results allow ranking of dimensions of variation based on genome-wide contribution to variance and enables analysis of gene-level results. this analysis uses the fact that individuals within monocytes have at least technical replicate, while individuals within t-cells have at least technical replicate. the variation across individuals within t-cells and monocytes is substantially larger than when the two cell types were combined. the fact that the contribution of individual varies between cell types is consistent with cell-specific regulatory genetics. finally, the fraction of variation explained by individual within each cell type at the gene-level is directly related to the probability of each gene having cis-eqtl within the corresponding cell type. analysis of gtex rna-seq dataset application of variancepartition to post mortem rna-seq data of multiple tissues tissues from the gtex consortium decouples the influence of multiple biological and technical drivers of expression variation. we analyzed experiments from individuals in tissues in order to restrict the analysis to tissues with rna-seq data for most individuals. variation across tissues is the major source of variation while the technical variables expression batch, ischemic time, rna isolation batch, and rin have a moderate effect on expression variation genome-wide. variation across expression batches is correlated with gc content but to a lesser degree that other datasets. cumulatively, these technical variables explain only % of the total expression variation. concerns about reliability of rna-seq data from post mortem samples has been raised due to the potential effects of rna degradation following cell death. variancepartition analysis indicates that variation in ischemic time has as relatively small effect genome-wide and the fraction of variance it explains is comparable to technical effects, yet the effect varies substantially across genes. fig. analysis of gtex dataset identifies drivers of expression variation at multiple levels. a violin and box plots of percent variation in gene expression explained by each variable. b results from variancepartition analysis allowing the contribution of individual to vary in each tissue. c probability of each gene having a cis-eqtl plotted against the percent variance explained by individual within each tissue. dashed lines indicate the genome-wide average probability, and curves indicate logistic regression smoothed probabilities as a function of the percent variance explained by individual within each tissue. the p-value indicates the probability that a more extreme coefficient relating the eqtl probability to percent variation explained by individual is observed under the null hypothesis. d fraction of variation in glmp explained by each source of variation. e glmp has a cis-eqtl active in blood but not skin the flexibility of the linear mixed model framework allows variancepartition to analyze cross-individual variation within each tissue. we note again that since the variance is analyzed within multiple subsets of the data, the total variation explained no longer sums to here. while variation across individuals explains only a median of % of variation when all tissues types are considered together, there is substantial variation across individuals within each tissue separately. cross-individual variation is highest in blood, while skin, blood vessel, and adipose tissue exhibit lower cross-individual variation. the fraction of variation explained by individual within each tissue is directly related to the probability of each gene having a cis-eqtl within the corresponding tissue. this association is not as strong as in other datasets likely due to the smaller number of individuals and to the relatively small fraction of variation across individuals in adipose tissue. at the gene-level, variancepartition can prioritize genes based on multiple criteria. for examples, glmp exhibits higher variation across individuals within blood but low variation in skin. this is consistent of a tissue-specific regulatory variation, and, in fact, the cis-eqtl rs influences gene expression in blood but not in skin. discussion as the scope of gene expression studies continues to expand, the need to quantify and interpret multiple drivers of expression variation is becoming essential. here we present variancepartition, a publicly available software package that leverages the power of the linear mixed model to quantify the contribution of multiple sources of variation in complex gene expression datasets. for each gene, this analysis partitions the total expression variance into the fraction attributable to each aspect of the study design. a variancepartition analysis gives a genome-wide summary of the drivers of variation, but also produces gene-level results to identify genes that deviate from the genome-wide trend. the fraction of expression variation is easily interpretable across genes, drivers of variation and datasets. thus variancepartition produces a more detailed and quantitative genome-wide overview than the standard principal components analysis and hierarchical clustering approaches. pca and hc focus on the major axis of variation, and they overlook the secondary drivers of variation that can be well characterized with variancepartition. moreover, the gene-level results from variancepartition indicate genes that deviate from the genome-wide trend and integration with additional data can enable a further interpretation. while pca and hc do not give gene-level results, differential expression analysis reports gene-level fold change and corresponding p-value for each aspect of the study design. yet de analysis does not produce a clear genome-wide summary, and the fold change and p-values are not easily comparable across multiple drivers of variation. analysis of publicly available gene expression studies demonstrate that variancepartition recovers striking patterns of biological and technical variation that are reproducible across multiple datasets. at a genome-wide level, expression variation across individuals and cell types is large enough to overcome the technical variation of transcriptome profiling. yet at the gene-level there is substantial deviation from the genome-wide trend due to a range of biological and technical factors. by quantifying the variance attributable to each aspect of the study design, variancepartition facilitates the interpretation of these gene-level effects in the context of additional information. we demonstrate reproducible findings that cross-individual variation is driven by cis-eqtl’s and technical variation across laboratories associated with gc content. moreover, variation across individuals and the relationship to cis-eqtl’s depend on the cell or tissue type.  the variancepartition workflow and implementation makes the rich linear mixed model framework easily applicable for interpreting drivers of variation in complex gene expression data. variancepartition provides a general statistical and visualization framework for studying drivers of variation in rna-seq datasets in many types of high-throughput genomic assays including rna-seq, protein quantification, metabolite quantification, metagenomic assays, methylation arrays and epigenomic sequencing assays. although we have focused here on large-scale studies, variancepartition analysis has given valuable insight into rna-seq datasets with as few as samples. the variancepartition software is an open source r package and is freely available on bioconductor. the software can easily be applied to rna-seq quantifications from featurecounts, htseq, kallisto, sailfish, salmon, rsem and cufflinks which have been processed in r with limma/voom, deseq, tximport and ballgown. the software provides a user-friendly interface for analysis and visualization with extensive documentation, and will enable routine application to a range of genomics datasets. availability and requirements project name: variancepartition project home page: http://bioconductor.org/packages/variancepartition operating systems: linux, mac os x, windows programming language: r ≥ v other requirements: bioconductor ≥ v license: gpl-2 additional file additional file supplementary note. additional details describing statistical methods and software. includes results from simulation study and additional figures from data analysis. abbreviations anovaanalysis of variance dedifferential expression eqtlexpression quantitative trait locus fvefraction of variance explained hchierarchical clustering iccintra-class correlation iqrinter-quartile range pcaprincipal components analysis remlrestricted maximum likelihood 
using a linear mixed model, variancepartition quantifies variation in each expression trait attributable to differences in disease status, sex, cell or tissue type, ancestry, genetic background, experimental stimulus, or technical variables. analysis of four large-scale transcriptome profiling datasets illustrates that variancepartition recovers striking patterns of biological and technical variation that are reproducible across multiple datasets. we describe a statistical and visualization framework, variancepartition, to prioritize drivers of variation based on a genome-wide summary, and identify genes that deviate from the genome-wide trend. our open source software, variancepartition, enables rapid interpretation of complex gene expression studies as well as other high-throughput genomics assays. as large-scale studies of gene expression with multiple sources of biological and technical variation become widely adopted, characterizing these drivers of variation becomes essential to understanding disease biology and regulatory genetics.
 as large-scale studies of gene expression with multiple sources of biological and technical variation become widely adopted, characterizing these drivers of variation becomes essential to understanding disease biology and regulatory genetics.  we describe a statistical and visualization framework, variancepartition, to prioritize drivers of variation based on a genome-wide summary, and identify genes that deviate from the genome-wide trend. using a linear mixed model, variancepartition quantifies variation in each expression trait attributable to differences in disease status, sex, cell or tissue type, ancestry, genetic background, experimental stimulus, or technical variables. analysis of four large-scale transcriptome profiling datasets illustrates that variancepartition recovers striking patterns of biological and technical variation that are reproducible across multiple datasets.  our open source software, variancepartition, enables rapid interpretation of complex gene expression studies as well as other high-throughput genomics assays. variancepartition is available from bioconductor: http://bioconductor.org/packages/variancepartition. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords transcriptome profilingrna-seqlinear mixed modelhttp://dx.doi.org/ /100000050national heart, lung, and blood instituteu hl107388-04schadt eric e. http://dx.doi.org/ /100000049national institute on agingu ag046170-02schadt eric e. issue-copyright-statement© the author 2016 