# An Empirical Bayesian Method for Detecting Differentially Expressed Genes Using EST Data

Authors: Na You; Junmei Liu; Chang Xuan Mao
Journal: International Journal of Plant Genomics
 (2008)

---

## Abstract

Detection of differentially expressed genes from  expressed sequence tags (ESTs) data has received much attention. An empirical Bayesian method is introduced in which gene expression patterns are estimated and used to define detection statistics. Significantly differentially expressed genes can be declared given detection statistics. Simulation is done to evaluate the performance of proposed method. Two real applications are studied.


---

## Body

## 1. Introduction

It is important
to detect differentially expressed genes, for example, exploring the key genes
related to certain diseases. As the EST sequencing technology develops, a large
number of EST databases from a variety of tissues are available. Enormous EST
collections provide opportunities to quantify gene expression levels [1].
Efficient statistical methods are in great demand.Several methods have been proposed to detect
significantly differentially expressed (SDE) genes from EST data [2]. Fisher's
exact test was used by the Cancer Genome Anatomy Project [3]. Audic and
Claverie [4] developed a Bayesian method. GT statistic [5] and R statistic [6]
were proposed for multilibrary comparison. In each method, gene-specific
detection statistics quantify differences of gene expression levels and SDE
genes are declared by their rankings.An empirical Bayesian method is proposed to detect SDE
genes. The relative gene expression abundances are estimated in each library,
and a new detection statistic is derived for each gene. In Section2,
simulation experiments suggest that the proposed method outperforms those
existing methods. Real applications are also studied in Section 2. Statistical
methods are described in Section 3. The possibility of extending the method for
multiple libraries is indicated in Section 4.


## 2. Results

Let(π11,π12,…,π1c) and (π21,π22,…,π2c) be the gene
expression patterns in two libraries, where πji is the relative
abundance of gene i in library j. The absolute difference between relative abundances
is Di=|π1i−π2i|. Given a sample of ESTs from library j, an empirical Bayes estimator π^ji for πji is defined in Section 3. Given gene i seen in both
samples, define D^i=|π^1i−π^2i|. Given gene i seen in only
one sample, for example, sample 2, define D^i=|π^1i−π^2i| if π^1i<π^2i and D^i=0 otherwise,
which is conservative in the sense that D^i possibly
underestimates Di. Gene i is declared to
be SDE if D^i is relatively
large.


### 2.1. Simulation

In a simulation experiment, EST frequencies are
generated from a multinomial distribution with sample sizesj and probability
vector (πj1,πj2,…,πjc), where c=1000, πji=λji/∑k=1cλjk, (λ11,λ12,…,λ1c) from G1, (λ21,λ22,…,λ2c) from G2, and G1 and G2 are two
distributions over (0,∞). The proposed methods, Fihser's exact test, χ2 test, AC
statistic, and R statistic, are studied. Given a cutoff point τ, the efficiency of a statistical method is measured
by pτ, the expected percentage of the true first τ SDE genes being
correctly declared as the first τ SDE genes. The
average of estimated pτ is calculated
from 500 replications.In the first four experiments,s1=s2=2000 and the results
are presented in Figure 1. Note that G1=U(0,10), Beta(2,5), 0.2δ(2)+0.4δ(5)+0.2δ(10), Gamma(3,0.1) and G2=Beta(2,1), Beta(2,5), Beta(2,2), Beta(2,2), respectively,
where U(a,b) is the uniform
distribution on (a,b), δ(a) is degenerate
at a, Beta(a,b) is transformed
from the beta distribution with shape parameters a and b by λ=p/(1−p) for p∈(0,1), and Gamma(a,b) is the gamma
distribution with shape a and scale b. For each cutoff point τ=10,20,…,100,pτ are calculated.
Clearly the proposed method has better performance than others.Simulation results of Fisher's exact test (◦), χ2 test (Δ), AC statistic
(+), R statistic
(×), and the
proposed statistic (•) in detecting
SDE genes using two EST samples of the same size.
(a)(b)(c)(d)In the second four experiments,(s1,s2)=(2000,4000), (4000,2000), (2000,4000), and (4000,2000), respectively,
and the results are presented in Figure 2. Note that G1=Gamma(3,0.1) and G2=Beta(2,2) in Figures 2(a)
and 2(b) and G1=U(0,10) and G2=Beta(2,1) in Figures 2(c)
and 2(d). The proposed method is usually the best one among all methods
studied.Simulation results of Fisher's exact test (◦), χ2 test (Δ), AC statistic
(+), R statistic
(×), and the
proposed statistic (•) in detecting
SDE genes using two EST samples of different sizes.
(a)(b)(c)(d)


### 2.2. Real Applications

One example
concerns Chinese spring wheat drought stressed leaf cDNA library (7235) and
root cDNA library (#ASP), available at TIGR gene indexes database (downloaded
athttp://www.tigr.org/tdb/tgi, 01/06/2006). In each EST sample, there are
totally 790 and 1306 sequenced ESTs, respectively. After removing the
unannotated 103 and 194 ESTs, the annotated ESTs are clustered into 465 and 804
groups with each group associated with a unique gene. Only those well-annotated
ESTs are used. The first 20 SDE genes by the proposed method are listed in
Table 1, among which 7, 7, 7, and 7 genes are in the set of first 20 SDE genes
by Fisher's exact test, χ2 test, AC
statistic, and R statistic, respectively.Table 1
The first 20 SDE genes in wheat leaf and root
libraries by the proposed method (x1i-the EST number
of gene i from leaf
library, x2i-that from root
library, 0/1-absence/presence in the set of the first 20 SDE genes).
Genex1ix2i1000DiFisherχ2ACRTC2495319027.101111TC23443827.881111TC23215184.870000TC26419184.870000TC26431184.870000TC24980503.401111TC23786062.620000TC26436062.620000TC24819062.620000TC264557121.850000TC23314151.590000TC24981151.590000TC24795051.570000TC24804051.570000TC26553051.570000TC26356411.370000TC23560401.371111TC24669401.371111TC24679401.371111TC26379401.371111Another example concerns pinus gene expression level
comparison in root gravitropism April 2003 test library (#FH3) and root control
2 (late) library (#FH4), also from TIGR, in which 2513 and 1132 ESTs associated
with 1211 and 605 genes are well annotated and clustered. Table2 lists the
first 20 SDE genes by the proposed method, among which 4, 4, 5, and 3 genes are
in the set of the first 20 SDE genes by Fisher's exact test, χ2 test, AC
statistic, and R statistic, respectively.Table 2
The first 20 SDE genes in #FH3 and #FH4 by the
proposed method (x1i-the EST number
of gene i from #FH3, x2i-that from
#FH4, 0/1-absence/presence in the set of the first 20 SDE genes).
Genex1ix2i1000DiFisherχ2ACRTC40351497.621111TC403556106.251110TC517791925.120010TC40566774.030000TC51682774.030000TC462901423.790000TC463721553.700000TC407681333.400000TC409121343.360000TC519951232.940000TC404201252.560000TC5170818122.440000TC404051142.430000TC40361062.341111TC46426062.341111TC4627619122.120000TC40388911.820000TC40647921.820000TC40350931.810000TC40731821.630000


## 2.1. Simulation

In a simulation experiment, EST frequencies are
generated from a multinomial distribution with sample sizesj and probability
vector (πj1,πj2,…,πjc), where c=1000, πji=λji/∑k=1cλjk, (λ11,λ12,…,λ1c) from G1, (λ21,λ22,…,λ2c) from G2, and G1 and G2 are two
distributions over (0,∞). The proposed methods, Fihser's exact test, χ2 test, AC
statistic, and R statistic, are studied. Given a cutoff point τ, the efficiency of a statistical method is measured
by pτ, the expected percentage of the true first τ SDE genes being
correctly declared as the first τ SDE genes. The
average of estimated pτ is calculated
from 500 replications.In the first four experiments,s1=s2=2000 and the results
are presented in Figure 1. Note that G1=U(0,10), Beta(2,5), 0.2δ(2)+0.4δ(5)+0.2δ(10), Gamma(3,0.1) and G2=Beta(2,1), Beta(2,5), Beta(2,2), Beta(2,2), respectively,
where U(a,b) is the uniform
distribution on (a,b), δ(a) is degenerate
at a, Beta(a,b) is transformed
from the beta distribution with shape parameters a and b by λ=p/(1−p) for p∈(0,1), and Gamma(a,b) is the gamma
distribution with shape a and scale b. For each cutoff point τ=10,20,…,100,pτ are calculated.
Clearly the proposed method has better performance than others.Simulation results of Fisher's exact test (◦), χ2 test (Δ), AC statistic
(+), R statistic
(×), and the
proposed statistic (•) in detecting
SDE genes using two EST samples of the same size.
(a)(b)(c)(d)In the second four experiments,(s1,s2)=(2000,4000), (4000,2000), (2000,4000), and (4000,2000), respectively,
and the results are presented in Figure 2. Note that G1=Gamma(3,0.1) and G2=Beta(2,2) in Figures 2(a)
and 2(b) and G1=U(0,10) and G2=Beta(2,1) in Figures 2(c)
and 2(d). The proposed method is usually the best one among all methods
studied.Simulation results of Fisher's exact test (◦), χ2 test (Δ), AC statistic
(+), R statistic
(×), and the
proposed statistic (•) in detecting
SDE genes using two EST samples of different sizes.
(a)(b)(c)(d)


## 2.2. Real Applications

One example
concerns Chinese spring wheat drought stressed leaf cDNA library (7235) and
root cDNA library (#ASP), available at TIGR gene indexes database (downloaded
athttp://www.tigr.org/tdb/tgi, 01/06/2006). In each EST sample, there are
totally 790 and 1306 sequenced ESTs, respectively. After removing the
unannotated 103 and 194 ESTs, the annotated ESTs are clustered into 465 and 804
groups with each group associated with a unique gene. Only those well-annotated
ESTs are used. The first 20 SDE genes by the proposed method are listed in
Table 1, among which 7, 7, 7, and 7 genes are in the set of first 20 SDE genes
by Fisher's exact test, χ2 test, AC
statistic, and R statistic, respectively.Table 1
The first 20 SDE genes in wheat leaf and root
libraries by the proposed method (x1i-the EST number
of gene i from leaf
library, x2i-that from root
library, 0/1-absence/presence in the set of the first 20 SDE genes).
Genex1ix2i1000DiFisherχ2ACRTC2495319027.101111TC23443827.881111TC23215184.870000TC26419184.870000TC26431184.870000TC24980503.401111TC23786062.620000TC26436062.620000TC24819062.620000TC264557121.850000TC23314151.590000TC24981151.590000TC24795051.570000TC24804051.570000TC26553051.570000TC26356411.370000TC23560401.371111TC24669401.371111TC24679401.371111TC26379401.371111Another example concerns pinus gene expression level
comparison in root gravitropism April 2003 test library (#FH3) and root control
2 (late) library (#FH4), also from TIGR, in which 2513 and 1132 ESTs associated
with 1211 and 605 genes are well annotated and clustered. Table2 lists the
first 20 SDE genes by the proposed method, among which 4, 4, 5, and 3 genes are
in the set of the first 20 SDE genes by Fisher's exact test, χ2 test, AC
statistic, and R statistic, respectively.Table 2
The first 20 SDE genes in #FH3 and #FH4 by the
proposed method (x1i-the EST number
of gene i from #FH3, x2i-that from
#FH4, 0/1-absence/presence in the set of the first 20 SDE genes).
Genex1ix2i1000DiFisherχ2ACRTC40351497.621111TC403556106.251110TC517791925.120010TC40566774.030000TC51682774.030000TC462901423.790000TC463721553.700000TC407681333.400000TC409121343.360000TC519951232.940000TC404201252.560000TC5170818122.440000TC404051142.430000TC40361062.341111TC46426062.341111TC4627619122.120000TC40388911.820000TC40647921.820000TC40350931.810000TC40731821.630000


## 3. Methods

Suppose that
there arec genes in a
library. Let xi be the number
of ESTs from gene i, a Poisson variable with mean λi. Given a prior distribution G on the λi, the posterior  mean of λi is E(λi∣xi)=(xi+1)hG(xi+1)/hG(xi), where hG(x)=∫λx/x!e−λdG(λ) is a Poisson
mixture. A gene is observed if and only if xi≥1. Conditioning on xi≥1, xi follows a
zero-truncated Poisson mixture hG(x)/(1−hG(0)) or a mixture fQ(x) of truncated
Poisson, where (1)fQ(x)=hG(x)1−hG(0)=∫λxx!(eλ−1)dQ(λ),dQ(λ)=(1−e−λ)dG(λ)∫(1−e−η)dG(η).
Let θ(Q)=hG(0)/(1−hG(0)) be the odds
that a gene is unseen. Write E(λi∣xi)=fQ(1)/θ(Q) if xi=0 and E(λi∣xi)=(xi+1)fQ(xi+1)/fQ(xi) otherwise.Letnx denote the
number of genes with exactly x ESTs in the
sample. The nonparametric maximum likelihood estimator Q^ for Q is
(2)Q^=argmax∑x≥1nxlog⁡fQ(x),
whose calculation is discussed in [7]. It is difficult
to estimate θ(Q) well [8]. There
are lower bound estimators, for example, θ(Q)^=n1(n1−1)/{2n(n2+1)} [9], where n=∑x≥1nx is the number of observed expressed genes. An empirical Bayes
estimator for λi is(3)λ^i=E(λi∣xi)^={fQ^(1)θ(Q)^,xi=0,(xi+1)fQ^(xi+1)fQ^(xi),xi≥1.As the relative
abundance πi satisfies πi=λi/∑k=1cλk, let π^i=λ^i/s^, where(4)s^=∑k=1c^λ^k=nfQ^(1)+∑x≥1nx(x+1)fQ^(x+1)fQ^(x),c^=n{1+θ(Q)^}.


## 4. Discussion

A new statistical method is proposed to compare the
gene expression patterns in two cDNA libraries. It can be extended to
multilibrary comparison, for example, considering all pairwise comparisons
among multiple libraries [3].



---

*Source: 817210-2008-03-13.xml*