Transcript 2000.stat - University of Wisconsin–Madison

Mining
for Low-abundance Transcripts
in Microarray Data
Yi Lin1, Samuel T. Nadler2,
Alan D. Attie2, Brian S. Yandell1,3
1Statistics, 2Biochemistry, 3Horticulture,
University of Wisconsin-Madison
September 2000
Basic Idea
• DNA microarrays present tremendous
opportunities to understand complex processes
• Transcription factors and receptors expressed at
low levels, missed by most other methods
• Analyze low expression genes with normal scores
• Robustly adapt to changing variability across
average expression levels
• Basis for clustering and other exploratory methods
• Data-driven p-value sensitive to low abundance &
changes in variability with gene expression
DNA Microarrays
• 100-100,000 items per “chip”
– cDNA, oligonucleotides, proteins
– dozens of technologies, changing rapidly
• expose live tissue from organism
– mRNA  cDNA  hybridize with chip
– read expression “signal” as intensity (fluorescence)
• design considerations
– compare conditions across or within chips
– worry about image capture of signal
Low Abundance Genes
• background adjustment
– remove local “geography”
– comparing within and between chips
• negative values after adjustment
– low abundance genes
• virtually absent in one condition
• could be important genes: transcription factors, receptors
– large measurement variability
• early technology (bleeding edge)
• prevalence across genes on a chip
– 0-20% per chip
– 10-50% across multiple conditions
Log Transformation?
• tremendous scale range in mean intensities
– 100-1000 fold common
• concentrations of chemicals (pH)
– fold changes have intuitive appeal
• looks pretty good in practice
• want transformed data to be roughly normal
– easy to test if no difference across conditions
– looking for genes that are “outliers”
• beyond edge of bell shaped curve
• provide formal or informal thresholds
Exploratory Methods
• Clustering methods (Eisen 1998, Golub 1999)
• Self-organizing maps (Tamayo 1999)
– search for genes with similar changes across conditions
– do not determine significance of changes in expression
– require extensive pre-filtering to eliminate
• low intensity
• modest fold changes
– may detect patterns unrelated to fold change
• comparison of discrimination methods (Dudoit 2000)
Confirmatory Methods
• ratio-based decisions for 2 conditions (Chen 1997)
– constant variance of ratio on log scale, use normality
• anova (Kerr 2000, Dudoit 2000)
– handles multiple conditions in anova model
– constant variance on log scale, use normality
• Bayesian inference (Newton 2000, Tsodikov 2000)
– Gamma-Gamma model
– variance proportional to squared intensity
• error model (Roberts 2000, Hughes 2000)
– variance proportional to squared intensity
– transform to log scale, use normality
0. acquire data
Q, B
7. standardize
Z=Y – center
spread
2. rank order genes
R=rank(A)/(n+1)
3. normal scores
N=qnorm(R)
4. contrast
conditions
Y=N1 – N2
Y = contrast
1. adjust for
background
A=Q – B
6. center & spread
X = mean
5. mean intensity
X=mean(N)
Normal Scores Procedure
adjusted expression
rank order
normal scores
A=Q–B
R = rank(A) / (n+1)
N = qnorm( R )
average intensity
difference
variance
standardization
X = (N1+N2)/2
Y = N1 – N2
Var(Y | X) 2(X)
S = [Y – (X)]/(X)
Motivate Normal Scores
• natural transformation to normality is log(A)
• background intensity B = bd +B
– measured with error 
– attenuation d may depend on condition
• gene measurement Q = [exp(g+h+)+b]d+Q
–
–
–
–
gene signal g
degree of hybridization h:
intrinsic noise  (variance may depend on g)
attenuation  (depends thickness of sample, etc.)
• subtract background: A = Q – B
– adjusted measurements A = dexp(g+h+)+
– symmetric measurement error =B – Q
Motivate Normal Scores (cont.)
• adjusted measurements: A = Q – B = dG+
• log expression level: log(G) = g+h+
• gene signal g confounded with hybridization h
– unless hybridization h independent of condition
• G is observed if
– no measurement error 
– no dye or array effect (no attenuation d)
– no background intensity
• natural under this model to consider N=log(G)
– normal scores almost as good
Robust Center & Spread
• genes sorted based on X
• partitioned into many (about 400) slices
– containing roughly the same number of genes
• slices summarized by median and MAD
– MAD = median absolute deviation
– robust to outliers (e.g. changing genes)
• MAD ~ same distribution across X up to scale
– MADi = i Zi, Zi ~ Z, i = 1,…,400
– log(MADi ) = log(i) + log( Zi)
• median ~ same idea
Robust Center & Spread
• MAD ~ same distribution across X up to scale
– log(MADi ) = log(i) + log( Zi), I = 1,…,400
• regress log(MADi) on Xi with smoothing splines
– smoothing parameter tuned automatically
• generalized cross validation (Wahba 1990)
• globally rescale anti-log of smooth curve
– Var(Y|X)  2(X)
• can force 2(X) to be decreasing
• similar idea for median
– E(Y|X)  (X)
Motivation for Spread
• log expression level: log(G) = g+h+
– hybridization h negligible or same across conditions
– intrinsic noise  may depend on gene signal g
• compare two conditions 1 and 2
– Y = N1 – N2  log(G1) – log(G2) = g1 – g2+ 1 -2
• no differential expression: g1=g2=g
– Var(Y|g) = 2(g)
– g  X suggests condition on X instead, but:
• Var(Y | X) not exactly 2(X)
• cannot be determined without further assumptions
Simulation of Spread Recovery
•
•
•
•
10,000 genes: log expression Normal(4,2)
5% altered genes: add Normal(0,2)
no measurement error, attenuation
estimate robust spread
Robust Spread Estimation
Q-Q Plot
5
0
-5
0.2
spread
0.4
0.6
Sample Quantiles
0.8
Estimate of Spread
-4
-2
0
2
mean expression
4
-4
-2
0
2
Theoretical Quantiles
4
Bonferroni-corrected p-values
• standardized normal scores
– S = [Y – (X)]/(X) ~ Normal(0,1) ?
– genes with differential expression more dispersed
• Zidak version of Bonferroni correction
– p = 1 – (1 – p1)n
– 13,000 genes with an overall level p = 0.05
• each gene should be tested at level 1.95*10-6
• differential expression if S > 4.62
– differential expression if |Y – (X)| > 4.62(X)
• too conservative? weight by X?
– Dudoit (2000)
Simulation Study
• simulations with two conditions
– 10,000 genes
– g1 ,g2 ~ Normal(4,2) for nonchanging genes
• 5% with differential expression
– gc ~ Normal(4,2) + Normal(3-rank(X)/(n+1),1/2)
– up- or down-regulated with probability 1/2
• intrinsic noise  ~ Normal(0,0.5)
• attenuation  = 1
• measurement error variance = 0, 1, 2, 5, 10, 20
500
100
200
300
400
0
12
5
10
20
0
Number of Changed Genes Captured
Success in Capture
0
100
200
300
400
(a) Rank by Normal Scores
500
Comparison of Methods
• differential expression for two conditions
• Newton (2000) J Comp Biol
– Gamma-Gamma-Bernouli model
– Bayesian odds of differential expression
• Chen (1997)
– constant ratio of expressions
– underlying log-normality
• normal scores (Lin 2000)
– some (unknown) transformation to normality
– robust, smooth estimate of spread & center
– Bonferroni-style p-values
Capturing Changed Genes: No Noise
Fold Change
0.5 1.0 2.0
5.0
0
0.1 0.2
0
0.02
0.10 0.50 2.00 10.00 50.00
No Noise: Average Intensity
Capturing Changed Genes: Noise
Fold Change
0.20
1.00
5.00
0
0.05
0
0.05
0.20
1.00
5.00 20.00
Noise 10: Average Intensity
Comparison with E. coli Data
• 4,000+ genes (whole genome)
• Newton (2000) J Comp Biol
– Bayesian odds of differential expression
• IPTG-b known to affect only a few genes
– ~150 genes at low abundance
– including key genes
E. coli with IPTG-b
Normal Scores IPTG-b
Cy3/Cy5 ratio
5.0
50.0
Cy3/Cy5 ratio
1 e-03 1 e-01
1 e+03
0
500.0
Raw IPTG-b
0
0.5
0
0.1
0
1 e-02
1 e+00
Mean Intensity
1 e+02
0.05
0.20
1.00
5.00
Mean Intensity
20.00
Diabetes & Obesity Study
• 13,000+ gene fragments (11,000+ genes)
– oligonuleotides, Affymetrix gene chips
– mean(PM) - mean(NM) adjusted expression levels
• six conditions in 2x3 factorial
– lean vs. obese
– B6, F1, BTBR mouse genotype
• adipose tissue
– influence whole-body fuel partitioning
– might be aberrant in obese and/or diabetic subjects
• Nadler (2000) PNAS
Low Abundance Obesity Genes
• low mean expression on at least 1 of 6 conditions
– negative adjusted values
– ignored by clustering routines
• transcription factors
– I-kB modulates transcription - inflammatory processes
– RXR nuclear hormone receptor - forms heterodimers
with several nuclear hormone receptors
• regulation proteins
– protein kinase A
– glycogen synthase kinase-3
• roughly 100 genes
– 90 new since Nadler (2000) PNAS
Table 1. Genes Differentially Expressed with Obesity as Identified by
the Normal Scores Procedure
Accesion #
Description
aa608251
Mus musculus Dp1l1 mRNA for polyposis locus protein 1-like 1 (TB2 protein-like
1)
M.musculus Gal beta-1,3-GalNAc-specific GalNAc alpha-2,6-sialyltransferase
gene
No significant homologies
Mus musculus cytochrome B561 (mcyt)
M.musculus gene for beta-3-adrenergic receptor
Mouse mRNA for transition protein 1 TP1
No significant homologies
Mus musculus protamine
Mus musculus transition protein 1 (Tnp1)
M.musculus tyrosinase-related protein-2
Mus musculus serum response factor (Srf)
Mus musculus hepsin (Hpn)
No significant homologies
Mus musculus mRNA for SRp25 nuclear protein
Mouse vas deferens androgen related protein (MVDP)
No significant homologies
Mouse cAMP-dependent protein kinase type II regulatory subunit
No significant homologies
Mouse 2',3'-cyclic-nucleotide 3'-phosphodiesterase
No significant homologies
No significant homologies
Mus musculus Balb/c mammary-derived growth inhibitor (MDGI)
No significant homologies
Mouse C33/R2/IA4
Mus musculus adipsin (Adn)
Mus musculus protamine 2
x93999
w45955
U16297
X72862*
X12521
aa544897
AA062269
AA144629
X63349
AA116710
AA105229
aa267014
aa408365
j05663
w49331
J02935
aa190005
m31810
AA033381
aa616759
u02883
AA198324
D14883
W36455
m27501
Obese/Lean
Lean/Obese
Average
p-value
0
940.01
0.535
0.0000
0.01
74.78
0.513
0.0000
0.01
0.01
0.01
0.01
0.01
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.04
0.04
0.04
80.52
95.95
104.54
105.3
138.69
41.58
42.71
43.42
57.02
59.34
62.72
65.53
30.75
31.1
32.02
34.86
35.21
36.63
36.67
38.1
39.31
25.01
25.63
27.47
0.356
0.506
1.575
0.751
0.346
0.877
0.71
1.599
0.847
0.318
1.154
0.38
0.682
0.849
0.399
1.056
0.512
0.506
0.522
1.002
0.428
1.134
12.089
0.99
0.0001
0.0000
0.0000
0.0000
0.0000
0.0009
0.0008
0.0000
0.0001
0.0047
0.0000
0.0001
0.0105
0.0093
0.0148
0.0015
0.0037
0.0028
0.0027
0.0010
0.0016
0.0119
0.0000
0.0143
X87096
AA048974
ET63455
X65026
AA138292
ET62360
w64688
AA035870
AA111277
w19022
w48392
AA168767
X99143
AA137962
W10926
AA103862
J04179
X03505
U09507
AA023099
u77460
x66224
AA154294
W85163
M.musculus brevican
No significant homologies
Mus musculus serum amyloid A-4 protein (Saa4)
M.musculus mRNA for GTP-binding protein
Similar to Rat S6 kinase
CALCIUM-DEPENDENT PHOSPHOLIPASE A2 PRECURSOR
No significant homologies
No significant homologies
Mus musculus hippocalcin-like 1 (Hpcal1)
No significant homologies
similar X75018 Id4 helix-loop-helix protein
Sequence withdrawn
M.musculus hair keratin, mHb6
similar to Human ras-related protein rab-14
Mus musculus mRNA for ubiquitin like protein
Similar to Human hqp0256
Mouse chromatin nonhistone high mobility group protein (HGM-I(Y))
Mouse gene exon 2 for serum amyloid A (SAA) 3 protein
Mus musculus p21 (Waf1)
Mus musculus dUTPase
Mus musculus anaphylatoxin C3a receptor
M.musculus retinoid X receptor-beta
Mus musculus protein tyrosine phosphatase (Ptpn16)
Mus musculus migration inhibitory factor
51.31
51.32
51.57
54.05
54.25
56.31
57.06
57.27
61.46
62.54
66.48
70.1
71.31
81.78
82.62
84.63
88.29
88.76
92.52
122.37
130.36
197.77
233.62
272.49
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0
0
0.34
0.565
0.31
0.81
0.419
0.417
0.943
1.091
0.348
0.694
0.573
0.09
0.203
1.076
0.762
0.268
1.138
0.693
0.724
0.221
1.191
0.256
0.306
0.168
0.0038
0.0001
0.0108
0.0001
0.0001
0.0001
0.0000
0.0000
0.0008
0.0000
0.0000
0.0072
0.0064
0.0000
0.0000
0.0013
0.0000
0.0000
0.0000
0.0002
0.0000
0.0000
0.0000
0.0000
2.00
0.50
0.20
0.05
Obese vs. Lean
5.00
Low Abundance Genes for Obesity
0.02
0.05
0.20
0.50
2.00
5.00
Average Intensity for Obesity
20.00
Microarray ANOVAs
• Kerr (2000)
• gene by condition interaction
– Nijk = genei + conditionj + gene*conditionij + rep errorijk
• conditions organized in factorial design
– experimental units may be whole or part of array
• genes are random effects
– focus on outliers (BLUPs), not variance components
– gene*conditionij = differential expression
– allow variance to depend on genei main effect
• replication to improve precision, catch gross errors
Microarray Random Effects
• variance component for non-changing genes
– robust estimate of MS(G*C) using smoothed MAD
– rescale normal score response N by spread (X)
– look for differential expression
• or use clustering methods
• variance component for replication
– robust estimate of MSE using smoothed MAD
– look for outliers = gross errors
-5
Obesity Dominance
0
5
10
Obesity Genotype Main Effects
-10
-5
0
Obesity Additive
5
Microarray QTLs
• condition may be genotype
– whole organism or pattern of genes
– genotype may be inferred rather than known exactly
• Nijk = genei + QTLj + gene*QTLij + individualijk
• QTL genotype depends on flanking markers
– mixture model across possible QTL genotypes
– single vs. multiple QTL
• single QTL may influence numerous genes
– epistasis = inter-genic interaction
– modification of biochemical pathway(s)