Transcript Presentation Title Goes Here

Significance analysis,
Microarray ANOVA, and
False Discovery Rates
European Institute of Statistical
Genetics
Liege, Belgium
September 3, 2007
Greg Gibson
Analysis Gene by Gene
Transcriptome Volcano plots
Significance (-log 10 P-value)
-15
A
A
-10
B
-5
D
C
C
0
Fold change (log 2 male/ female)
Variance estimators
• Permutation approach may be more appropriate where you have
many treatments with low replication
• Gene-specific approach means that the power for each gene
varies, but should be conservative
Gary Churchill
Variance components
Variance Component Contributions
Effect
Contribution of
Genotype to Variance
1.0
% positive
Mean Proportion
Sex + SG
74
0.60
Genotype + SG
51
0.33
Age
58
0.21
SexGenotype
32
0.31
Sex
61
0.56
Genotype
28
0.25
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1.0
Contribution of Sex to Variance
Variables in a Microarray Study
Arrays
Dyes
Always random; Technical replicates
Fixed; Usually 2 levels
Treatments
Genes
The classes of RNA samples
May be duplicated on a slide
Nested Factors: Individual within Variety
Spot within Gene
Two Approaches to MA-ANOVA
• Kerr and Churchill:
- Fit a single ANOVA model and assess significance of the
Gene×Treatment interaction against a common error term
- Use permutation testing to establish the appropriate error
• Wolfinger:
- Fit a global normalization model, then perform separate MMANOVA for each gene
- Use the error term from each ANOVA to assess significance
Kerr and Churchill ANOVA
Spots!
Statistical Linear Model:
yijkgr = µ + Ai + Dj + ADij + Gg + (VG)kg + (AG)igr + (DG)kg + eijkgr
Variety-by-Gene
effects
We assume that the error eijkgr is randomly distributed
with mean 0, and is not correlated among varieties.
Quantities of interest are expression levels of gene
specifically attributable to different varieties:
(VG)kg- (VG)k’g
Two-step Analysis of Variance
1. Express transcript level relative to sample mean
log(fluorescence) =  + Array + Dye + Residual
2. For each gene, assess significance of treatment effects on
the Residual (ie. relative expression level)
Residual =  + Spot + Sex + Geno + Treat + Interact + Error
(note: Fitting ‘Spot’ eliminates need for a Reference sample)
Mixed-Model ANOVA I
Stage 1: Fit Normalization Model
log2Y =  + array + dye + array*dye + error
Notes:
– Random effects are in red; others are fixed.
– Some preprocessing may be desirable (e.g. loess).
– This model averages across all genes.
Mixed-Model ANOVA
II
Stage 2: Analyze Results from Normalization Step
•
•
•
•
Quality control using parameter estimates and residuals
Inference on effects and variance components of interest
Construct data to be used for clustering
Use dynamic visualization for exploration of data and results
Mixed-Model ANOVA
III
Stage 3: Fit Gene-Specific Models
For each of the 1700 genes, construct residuals (Relative
Fluorescence Intensities, RFA) from Stage 1 and use them as
input to the following gene-specific model:
RFA =  + sex + drug + sex*drug +
animal + array(animal) + dye(array animal) + error
Advantages of MM-ANOVA
• Assessment of SIGNIFICANCE of effects, particularly in
multifactorial designs
• Comparison of variance contributions
• Does away with ‘wasteful’ reference samples
• Combined analysis of independent experiments,
including different platforms
• Robust statistical inference
Issues with ANOVA models
• Little agreement on normalization step
• Model fitting is not always straight-forward
• Appropriate assignment of error terms difficult
• How do we interpret significance levels?
– A popular misconception is that one must assume normality to
use ANOVA. However, non-parametric inferential methods
and/or Empirical Bayes methods can be combined with ANOVA
Mixed Model For Affy ProbeLevel Data (after Normalization)
log2 (P Mijk )  Ti  Pj  Ak  e ijk ,
i  1,...,4
j  0,...,10
k  1,...,31
Ak ~ N (0,  a2 ),e ijk ~ N (0,  2 )
Hypothesis Testing
True (but unknown) situation
Reported by the m Tests:
H0 is True
m0
H0 is False
m1
H0 was NOT rejected
m-S
True negative
correct
m0 – F
False negative
type II error
m1 - T
H0 WAS rejected
S
False positive
type I error
F
True positive
correct
T
Multiple Hypothesis Testing Example
True (but unknown) situation
Reported by 100 Tests:
H0 is True
90
H0 is False
10
H0 was NOT rejected
85
True negative
correct
80
False negative
type II error
5
H0 WAS rejected
15
False positive
type I error
10
True positive
correct
5
False Discovery Rates (FDR)
• The FDR measures the proportion of false
positives among all genes called significant:
# false positives
# called significant
F
F+T
10
15
67%
•
Typically, you want this ratio to be small, say
10%, implying that most positives are TRUE
•
The FDR gives the rate at which further
biological verification will result in dead ends
<
False Positive Rates (FPR)
• The FPR measures the rate at which truly null
genes are called significant:
# false positives
# truly null
F
m0
10
90
11%
•
The p-value is a measure of significance in
terms of the FPR (Type 1 error rate)
•
If 5 genes are called significant, but none are,
then the FPR would be 5/100 = 5%
q values
•
The q-value is a measure of significance in
terms of the FDR (Type 1 error rate)
•
A simple formulation is to assume that all of the
genes are truly null (Benjamini and Hochberg)
•
Storey recommends estimating the proportion of
truly null genes (Est π0 = Est m0/m) from the
observed distribution of p-values
•
Typically there are some truly positive genes, so
m0 < m, say π0 ~ 70%, and so the number of
truly null genes is reduced, allowing a more
accurate estimate of the FDR.
Estimating π0 at 
For example
•
Suppose you have 10,000 probes on an array,
and 100 are called significant at p <0.001.
•
You would expect 0.1%, or 10 genes to be
significant at this level by chance.
•
According to BH formulation, the FDR is 10 false
positives out of 100 called significant = 10%
•
But if 2,000 genes really are differentially
expressed then only 8,000 are true nulls, and
the corresponding number of false positives may
be 8, so the FDR would be 8/100 = 8%
Q-plots