Bayesian Random-Effects Meta-Analysis Models in Gene Expression Studies

Abstract

Random-effects meta-analysis models are commonly applied in combining effect sizes from individual gene expression studies. However, study heterogeneity is unknown and may arise from a variation of sample quality and experimental conditions. High heterogeneity of effect sizes can reduce the statistical power of the models. In addition, classical random-effects meta-analysis models are based on a normal approximation, which may be limited to small samples and its results may be biased toward the null value. A Bayesian approach was used to avoid the approximation and the biases. We applied a sample-quality weight to adjust the study heterogeneity in the Bayesian random-effects meta-analysis model with weighted between-study variance on a sample quality indicator and illustrated the application of this approach in Alzheimer’s gene expression studies.

1. Introduction

Advances in the development of high-throughput technologies have enabled researchers to identify and quantify a large range of gene expression differences in many diseases. The number of gene expression studies has been increasing over the past decades as a result of advanced technologies. Several of them examine and address the same biological questions, even they have been implemented under different experimental conditions. Meta-analysis of gene expression data, therefore, has become a widely applied approach in combining results from multiple studies to identify common expression patterns that sometimes cannot be detected in individual studies. The meta-analytic approach has been shown to increase statistical power, accuracy, and generalizability of results [1, 2, 3, 4]. The use of meta-analysis techniques depends on the type of response and study objectives and most analyses in microarray studies have emphasized identifying differentially expressed (DE) genes or genes that distinguish the group of samples.

Random-effects (RE) meta-analysis models are commonly applied in combining effect sizes from individual gene expression studies. However, study heterogeneity is unknown and may arise from the variation of sample quality and experimental conditions and the study heterogeneity can decrease the statistical power of the models. To maintain power, we can increase the number of studies [5] or apply an appropriate estimation method for incorporating study heterogeneity into the models. Typically, the classical RE models assume studies are independently and uniformly sampled from a population of studies. However, studies are possibly designed based on the results of previous studies. The independence assumption and an infinite population of studies may not exist. Bayesian random-effects (BRE) models have been applied to handle the uncertainty of parameters in the models. The uncertainty is incorporated through a prior distribution. A summary of evidence after the data have been observed is described by the likelihood of the models. Multiplying the prior distribution and the likelihood function will provide a posterior distribution of the parameters [6, 7].

Sample quality has a substantial impact on results of gene expression studies [8, 9]. Low heterogeneity can be found in meta-analyses containing good-quality samples, while poor-quality samples can induce high heterogeneity of effect sizes. We recently evaluated the relationships between DE and heterogeneous genes in meta-analyses of Alzheimer’s gene expression data. Some overlapped DE and heterogeneous genes were detected in meta-analyses containing borderline- or poor-quality samples, while no heterogeneous genes were identified in meta-analyses containing good-quality samples [10]. The data obtained from borderline- or poor-quality samples can increase study heterogeneity and decrease the efficiency of meta-analyses [11, 12].

Small samples in gene expression studies may limit the application of classical RE models and its results may be biased toward the null or the observed value is closer to the null hypothesis than the true value. The BRE model can be used to avoid the approximation and the biases. We introduced a meta-analytic approach that included a sample-quality weight to adjust study heterogeneity in the BRE model [13]. The gene expression data therefore would include both up-weighted good-quality samples and down-weighted borderline-quality samples. Therefore, we first reviewed the classical RE models, the BRE model, and the weighted BRE model in the methods section and then illustrated an application of the methods in Alzheimer’s gene expression studies. Our results are then compiled in the results section and followed by discussion and conclusions.

2. Methods

2.1 Standard random-effects model

Choi et al. [14] introduced an unbiased standardized mean difference in expression between groups for each gene g [14, 15]. The meta-analytic model for combining the effect sizes is based on a two-level hierarchical model as follows:

$$y_{\mathit{ig}}={\theta }_{\mathit{ig}}+{\epsilon }_{\mathit{ig}},{\epsilon }_{\mathit{ig}}\sim N\left(0,{\sigma }_{\mathit{ig}}^2\right)$$

$${\theta }_{\mathit{ig}}={\beta }_g+{\delta }_{\mathit{ig}},{\delta }_{\mathit{ig}}\sim N\left(0,{\tau }_g^2\right)$$E1

where $y_{\mathit{ig}}$ denotes the expression for gene g in study i = 1, …, k; θ_ig denotes the true difference in mean expression; ${\sigma }_{\mathit{ig}}^2$ denotes the within study variability displaying sampling errors conditional on $i\mathrm{th}$ study; ${\beta }_g$ denotes the common effects or the average measure of differential expression across individual datasets for each gene or the parameter of interest; ${\delta }_{\mathit{ig}}$ denotes the random effects; and ${\tau }_g^2$ denotes the between-study variability displaying the variability between studies. We estimate the parameter of the common effect using a weighted least squares estimation. We minimize the sum of squares error by differentiating with respect to ${\widehat{\beta }}_g$ for each gene in each study, which yields.

$${\widehat{\beta }}_g=\frac{{\sum }_{i=1}^k{w}_{\mathit{ig}}{y}_{\mathit{ig}}}{{\sum }_{i=1}^k{w}_{\mathit{ig}}},w_{\mathit{ig}}={\sigma }_{\mathit{ig}}^{-2}.$$E2

Generally, an unbiased estimator for ${\tau }_g^2$ can be derived from the method of moments and the estimator can attain negative values, therefore its truncated version, or the DerSimonian-Laird (DSL) estimator, is considered instead:

$${\widehat{\tau }}_{\mathit{DSL}\left(g\right)}^2={\widehat{\tau }}_g^2=max\left\{0,\frac{Q_g-\left(k_g-1\right)}{S_{1g}-\left(S_{2g}/S_{1g}\right)}\right\}$$E3

where $Q_g={\sum }_{i=1}^k{w}_{\mathit{ig}}{\left(y_{\mathit{ig}}-{\widehat{\beta }}_g\right)}^2$; $w_{\mathit{ig}}={\sigma }_{\mathit{ig}}^{-2}$; and $S_{\mathit{rg}}={\sum }_{i=1}^k{w}_{\mathit{ig}}^r$. There may be a small bias of the DSL estimator but the bias occurs where ${\widehat{\tau }}_g^2$ is close to zero or homogeneity [16], and the bias of the DSL estimator could not be traded off by variance reductions [17, 18]. Therefore, the DSL estimator is commonly applied when fitting random-effects models for a meta-analysis [19, 20]. In this study, we estimated ${\widehat{\beta }}_g\left({\widehat{\tau }}_g^2\right)$ for each gene from the microarray data with the weight $w_{\mathit{ig}}={\left({\widehat{\sigma }}_{\mathit{ig}}^2+{\widehat{\tau }}_g^2\right)}^{-1}$ by the generalized least squares method [14], providing the minimum variance unbiased estimator for ${\beta }_g$, to obtain the statistic for each gene $\left(z_g\right)$,

$${\widehat{\beta }}_{\mathit{DSL}\left(g\right)}\left({\widehat{\tau }}_{\mathit{DSL}\left(g\right)}^2\right)={\widehat{\beta }}_g\left({\widehat{\tau }}_g^2\right)=\frac{{\sum }_{i=1}^k{\left({\widehat{\sigma }}_{\mathit{ig}}^2+{\widehat{\tau }}_g^2\right)}^{-1}{y}_{\mathit{ig}}}{{\sum }_{i=1}^k{\left({\widehat{\sigma }}_{\mathit{ig}}^2+{\widehat{\tau }}_g^2\right)}^{-1}},$$E4

$$\mathit{Var}\left[{\widehat{\beta }}_{\mathit{DSL}\left(g\right)}\left({\widehat{\tau }}_{\mathit{DSL}\left(g\right)}^2\right)\right]=\mathit{Var}\left[{\widehat{\beta }}_g\left({\widehat{\tau }}_g^2\right)\right]=\frac{1}{{\sum }_{i=1}^k{\left({\widehat{\sigma }}_{\mathit{ig}}^2+{\widehat{\tau }}_g^2\right)}^{-1}},$$E5

$$\text{such that}{z}_g=\frac{{\widehat{\beta }}_{\mathit{DSL}\left(g\right)}\left({\widehat{\tau }}_{\mathit{DSL}\left(g\right)}^2\right)}{\sqrt{\mathit{Var}\left({\widehat{\beta }}_{\mathit{DSL}\left(g\right)}\left({\widehat{\tau }}_{\mathit{DSL}\left(g\right)}^2\right)\right)}}\sim N\left(0,1\right).$$E6

The standard random-effects model currently estimates the between-study variance $\left({\tau }_g^2\right)$ or the study heterogeneity using the DSL estimator.

3. Bayesian random-effects model (BRE)

In contrast to the classical RE model, the data and model parameters in the BRE model are considered to be random quantities [21]. We applied the BRE model to allow for the uncertainty of the between-study variance in this study. The model for gene g is written as

$$y_{\mathit{ig}}\left|{\theta }_{\mathit{ig}}\right.\sim N\left({\theta }_{\mathit{ig}},{\sigma }_{\mathit{ig}}^2\right),$$

$${\theta }_{\mathit{ig}}\left|{\beta }_g,{\tau }_g\right.\sim N\left({\beta }_g,{\tau }_g^2\right),$$

$${\beta }_g\sim N\left(0,1000\right),\text{and}$$

$${\tau }_g\sim \text{uniform}\left(0,1\right).$$E7

The kernel of the posterior distribution can be written as

$$\begin{array}{l}p\left({\beta }_g,{\theta }_{1g},\dots ,{\theta }_{\mathit{kg}},{\tau }_g^2\right)\propto p\left({\mathit{\theta }}_g|{\mathbf{y}}_g,{\mathbf{\sigma }}_g^2\right)p\left({\beta }_g,{\tau }_g^2|{\mathit{\theta }}_g\right)\\ \\ \propto {\prod }_{i=1}^{k}p\left({\theta }_{\mathit{ig}}|y_{\mathit{ig}},{\sigma }_{\mathit{ig}}^2\right)p\left({\theta }_{\mathit{ig}}|{\beta }_g,{\tau }_g^2\right)\pi \left({\beta }_g\right)\pi \left({\tau }_g^2\right),\end{array}$$E8

where ${\mathbf{y}}_g=\left(y_{1g},\dots ,y_{\mathit{kg}}\right)$, ${\mathbf{\sigma }}_g^2=\left({\sigma }_{1g}^2,\dots ,{\sigma }_{\mathit{kg}}^2\right)$, and ${\mathit{\theta }}_g=\left({\theta }_{1g},\dots ,{\theta }_{\mathit{kg}}\right)$ for gene g in the ith study; i = 1,…,k. The $\pi \left({\beta }_g\right)$ and $\pi \left({\tau }_g^2\right)$ are non-informative priors given as ${\beta }_g\sim N\left(0,1000\right)$, and ${\tau }_g\sim$ uniform (0,1).

The choice of prior distributions for scale parameters can affect analysis results, particularly in small samples. With scale parameters, the distributional form and the location of the prior distributions are obtained [22]. Uniform distributions are appropriate non-informative priors for ${\tau }_g^2$ [6, 13].

4. Sample-quality weights

The quality control (QC) criteria for identifying poor-quality samples in this study were the 3′:5′ GAPDH ratio greater than 3 and/or percent of present calls less than 30% for Affymetrix arrays; and detection rate less than 30% for Illumina Bead Arrays, in addition to data visualizations [8, 23]. Poor-quality samples were excluded before data preprocessing. Furthermore, the inverse of the within-study variance is considered an optimal weight for meta-analysis. The variance of weighted mean (${\widehat{\beta }}_g$) is minimized when the weights are taken from the variance of the samples $y_{\mathit{ig}}$. A high variance gives low weights in meta-analysis [24, 25]. In our recent study, the weight intermingled with the QC indicators called as “zero-to-one weight” was most appropriate for detecting DE genes [13]. The QC indicators adjusted the within-study variance in the weighted function as:

$$w_{P6}={\left({\sigma }_{\mathit{ig}}^{2\left(w_{P1}\right)}+{\widehat{\tau }}_g^2\right)}^{-1},$$E9

where $w_{P1}\in \left\{2^{-S_{\mathit{ij}}},0.01{\tilde{P}}_{\mathit{ij}}\right\}$${\tilde{P}}_{\mathit{ij}}$ denotes the percent of present calls of the jth sample in the ith study. A high value of the $P_{\mathit{ij}}$ weight indicates good-quality samples, providing high values of zero-to-one weights $\left(w_{P,\mathit{ij}}\right)$to give more weight to the expression data.

5. Weighted between-study variance model

We adjusted the between-study variance in the BRE model (Eq. (9)) by multiplying with an average weight over the total sample in the ith study for gene g$\left({\overline{w}}_{\mathit{ig}}={\sum }_{j=1}^{n_{\mathit{ig}\left(a\right)}+n_{\mathit{ig}\left(\mathrm{c}\right)}}{w}_{\mathit{ijg}}/\left(n_{\mathit{ig}\left(a\right)}+n_{\mathit{ig}\left(\mathrm{c}\right)}\right)\right)$. The BRE weighted between-study variance model for gene g is given by

$$y_{\mathit{ig}}\left|{\theta }_{\mathit{ig}}\right.\sim N\left({\theta }_{\mathit{ig}},{\sigma }_{\mathit{ig}}^2\right),$$

$${\theta }_{\mathit{ig}}\left|{\beta }_g,{\tau }_g{\overline{w}}_{\mathit{ig}}\right.\sim N\left({\beta }_g,{\tau }_g^2{\overline{w}}_{\mathit{ig}}\right),$$

$${\beta }_g\sim N\left(0,1000\right),\text{and}$$

$${\tau }_g\sim \text{uniform}\left(0,1\right).$$E10

We implemented two chains each with 20,000 iterations, a 15,000 burn-in period, and a thinning of 3 in the Bayesian model, and assessed the convergence of the model using the Gelman and Rubin diagnostic [26]. The posterior mean was standardized by posterior standard deviation as the posterior distribution was symmetric and normal. We then applied a Benjamini and Hochberg (BH) procedure to control the false discovery rate (FDR) for multiple gene testing. The performance of several BRE models for unweighted and weighted data, Gibbs and Metropolis-Hastings sampling algorithms, weighted common effect, and weighted between-study variance and classical RE models for unweighted and weighted data were evaluated in simulation studies [10, 13]. The classical RE and BRE meta-analysis models were implemented using programs from MAMA, R2jags, and metaDE packages in the R programming environment [27, 28, 29].

6. Results

We reviewed publicly available gene expression data from the NCBI GEO database. Twelve series of RNA expression profiling in the GEO database were selected for initial review. Eligible criteria for data acquisition were as follows: (1) the datasets were publicly accessible, (2) the samples were from human brain regions, (3) the series included samples from healthy controls, (4) the datasets included phenotypic data and published manuscripts describing the data, (5) the datasets without redundant samples, and (6) the raw or normalized intensity data were defined as gene expression levels. For each study we reviewed the minimum information about a microarray experiment (MIAME) from the GEO website, including research methods and results described in the manuscripts, and data summaries of the phenotypic data. This presented study included four publicly available Alzheimer’s disease (AD) gene expression datasets of post-mortem brain samples. The Gene Expression Omnibus accession numbers: GSE1297 [30], GSE5281 [31], GSE29378 [32], and GSE48350 [33] containing the gene expression and phenotypic data were included in this presented study. Some of these accession numbers (GSE5281, GSE13214, and GSE48350) include samples from multiple brain regions; we restricted our attention to only those samples acquired from hippocampus and to AD and control subjects in each dataset. The QC criteria for identifying poor-quality samples were having a 3′/5′ glyceraldehyde-3-phosphate dehydrogenase (GAPDH) ratio greater than three and/or percent of present calls less than 30% [23]. We then conducted within study data preprocessing, quantile normalization, and data aggregating. Our meta-analysis was therefore performed on 12,037 target genes in 131 subjects (68 AD cases and 63 controls) from the four studies using the Affymetrix and Illumina platforms (Figure 1). We then considered five ways of metadata sets and primarily examined the strength of study heterogeneity ($R^2$) of each metadata object as described in [10]. The metadata A, B, D, and E had a relatively high $R^2$, while the metadata C had a relatively low $R^2$. In other words, metadata C contains homogenous data, while the remaining metadata objects may contain heterogeneous data (Figure 2). The distribution of unbiased standardized mean differences of gene expression in the GSE5281 dataset, which is different from the other datasets, is presented in Figure 3. The percent of present calls and the 3′:5′ GAPDH ratio of the heterogeneous dataset is presented in Figure 4.

Figure 1.

Study profile for meta-analysis in Alzheimer’s gene expression datasets.

Figure 2.

Barplots on strength of study heterogeneity measuring by random effects coefficient determination (R²) in meta-analysis in Alzheimer’s gene expression datasets. The R² of 10,000 genes were categorized into five groups. Tentatively, R² close to 0.25, 0.50, and 0.75 indicate low, moderate, and high heterogeneity, respectively. The y-axis presents the R² and the x-axis presents the number of genes in the meta-analysis.

Figure 3.

Distribution of unbiased standardized mean difference of gene expression (x-axis) between Alzheimer’s and control groups in GSE1297, GSE5281, GSE29378, and GSE48350 datasets.

Figure 4.

Percentage of present calls and 3′/5′ GAPDH ratio of GSE5281samples.

In this meta-analysis on the metadata of the four AD gene expression datasets, 1766 DE genes were identified by the classical RE model, while 466 DE genes were identified by the weighted BRE model. Almost all the DE genes identified by the weighted BRE model were genes among the significant DE genes identified by the classical RE model. Figure 5 presents the heatmap of 1766 DE up-regulated and down-regulated genes detected in the AD samples. There was no trend apparently toward more up-regulated genes or down-regulated genes on the AD samples as compared to the control samples. Meanwhile, there was a trend toward more down-regulated genes on the AD samples as compared to the control samples in the heatmap of the 466 identified DE genes (Figure 6). The 446 genes could potentially be down-regulated genes that may contribute to the good classification of Alzheimer’s samples (Table 1).

Figure 5.

Heatmaps of expression patterns of 1766 differentially expressed genes in hippocampus in Alzheimer’s and control samples. The differentially expressed genes were detected by the classical random-effect meta-analysis model on the metadata of four Alzheimer’s gene expression datasets: GSE1297, GSE5281, GSE29378, and GSE48350.

Figure 6.

Heatmaps of expression patterns of 446 differentially expressed genes in hippocampus in Alzheimer’s and control samples. The differentially expressed genes were detected by the classical random-effect meta-analysis model on the metadata of four Alzheimer’s gene expression datasets: GSE1297, GSE5281, GSE29378, and GSE48350.

AACS, AASDHPPT, ABCA1, ACLY, ACOT7, ADAM22, ADAM23, ADARB1,AFF2, AGK, AMPH, ANGPT1, ANP32C, AP2S1, AP3B2, AP3D1, AP3M2, APBA2, APMAP, ARFGEF1, ARHGDIG, ARHGEF9, ARPC5L, ASIC2, ASNS, ASPHD1, ATAT1, ATP1A1, ATP1A3, ATP2A2, ATP2B2, ATP5B, ATP5C1, ATP5D, ATP5G1, ATP5H, ATP5L, ATP6AP1, ATP6V0B, ATP6V0E1, ATP6V1B2, ATP6V1E1, ATP6V1G2, ATP8A2, ATPIF1, ATR, ATRN, ATRNL1, ATXN7L3B, BCL2, BEX1, BEX4, BPGM, BSN, C10orf88, C12orf10, C14orf2, C16orf45, C1orf216, C2CD5, C2orf47, C5orf22, CA10, CABYR, CACNA2D3, CADPS, CALY, CAMK1, CAMK2N1, CAMKV, CAPRIN2, CCK, CDC40, CDC42EP4, CDK5, CDKN2D, CGREF1, CHGB, CHN1, CISD1, CLIP3, CLTA, CNR1, COPS3, COPS7A, COPZ2, COQ6, COX4I1, COX6C, CP, CREBBP, CRYM, CS, CUL2, CYCS, CYP4F12, DAP3, DCTN1, DDX41, DEAF1, DGUOK, DHRS11, DHRS3, DIRAS3, DLEC1, DLG2, DLGAP2, DMXL2, DNASE2, DNM1, DNM1L, DNM3, DOCK3, DOPEY1, DROSHA, DYNC1H1, DYNC1I1, ECM2, EEF1A2, EGFR, EHD3, ELF1, ELOVL4, ELOVL6, ENC1, ENO2, ENTPD2, ENTPD3, EPB41L1, EPS15, ERC2, FAM111A, FAM127A, FAM162A, FAM174B, FAM188A, FAM216A, FAM60A, FAM98A, FAR2, FGF12, FH, FHL2, FIBP, FKBP3, FMO2, FOCAD, FOXJ1, FOXO1, FRMPD4, FSD1, FXN, FYCO1, FYN, GABBR2, GABRG2, GAD, GAD2, GCC2, GLS2, GNAI2, GNG3, GNG4,GOT1, GPHN,GPI, GPRASP1, GRIA2, GRIN1, GRM1, GSTA4, GUCY1B3, GUK1v GYPC, HAGH, HARS, HERC1, HMGCR, HMP19, HN1, HNRNPUL1, HPRT1, HSPA12A, IGF1R, IMMT, IMP3, IMP4, INA, INPP5F, ITPKB, ITSN1, KAT6A, KCNN3, KCNQ2, KIAA0513, KIAA1324, KIF21B, KIFAP3, LARGE, LCMT1, LDB2, LEMD3, LGALS8, LPAR4, LPCAT4, LPIN1, LPP, LRPPRC, LRRC8B, LY6H, MAK16, MAP1A, MAP2K1, MAP2K4, MAP3K9, MAPK9, MAST3, MCF2, MCTS1, MDH1, MDH2, MICU1,MKKS, MLLT11, MOAP1, MPP1, MPPED2, MRPL15, MRPL17, MRPL35, MRPS11, MRPS17, MRPS22, MTMR11, MTSS1L, MTX2, MXI1, MYL12B, MYT1L, NAP1L2, NAP1L3, NCALD, NDN, NDRG3, NDRG4, NDUFA10, NDUFA3, NDUFA4, NDUFA8, NDUFA9, NDUFS3, NDUFS5, NDUFV2, NECAP1, NEDD8, NEFL, NEFM, NELL1, NETO2, NFIB, NIPSNAP3B, NLK, NME1, NMNAT2, NOVA1, NREP, NRGN, NRIP3, NRN1, NSF, NSG1, NUPL2, OGDHL, OPA1, ORC5, P4HTM, PAGE1, PAX6, PDCD1LG2, PEX11B, PIN1, PLCD1, PLCE1, PLCL2, PLD3, PLEC, PLEKHA4, PLK2, PLSCR4, PLXNB2, PMFBP1, PNMAL1, PNO1, PODXL2, POLB, POLRMT, POP7, PPFIA4, PPIA, PPIP5K1, PPM1H, PPME1, PPP1R13L, PPP2CA, PPP3CB, PREP, PREPL, PRKCZ, PRMT1, PRPF40A, PSD4, PSMD8, PTDSS1, PTGES2, PTPRE, PTPRR, PTRH2, PTS, PVRL3, RAB11A, RAB26, RAB27A, RAB2A, RAB6A, RAD51C, RAP1GDS1, RARS, RBFOX2, RGS17, RGS7, RHOQ, RIMBP2, RIT2, RND2, RNF123, RNF41, RNFT2, RNMT, RNPS1, RPH3A, RPP40, RPS6KC1, RUNDC3B, RWDD2A, RXRA, SCAMP2, SCG5, SCN2A, SCN3B, SDHA, SEC22A, SEC61A2, SEH1L, SEPT6, SERPINI2, SEZ6L2, SLC12A5, SLC25A11, SLC25A12, SLC25A14, SLC25A4, SLC4A1AP, SLIRP, SLITRK3, SMARCA4, SMO, SMOX, SMYD2, SNAP25, SNAP91, SNCB, SOX2, SPAG7, SPIN2A, SPINT2, SRM, SRPR, SS18L1, SSPN, STAU2, STMN2, STX6, STXBP1, SULT4A1, SUSD4, SV2B, SYDE1, SYN1, SYN2, SYNGR1, SYNJ1, SYT1, TAGLN3, TAZ, TBC1D31, TBC1D9, TBCC, TBCE, TBL1X, TBPL1, TCEA2, TCF7L2, TERF2IP, TGFBR3, THOC5, TMEM151B, TMEM160, TMEM246, TMEM59L, TMEM70, TMEM97, TNPO1, TOMM34, TOMM70A, TOR1A, TPD52, TPI1, TRAP1, TRAPPC2, TRIM37, TRIM9, TRIOBP, TSPAN13, TSPAN7, TSSC1, TUBA1B, TUBA4A, TUBB3, TUBG1, TUBG2, TXNDC9, UBE2M, UBE2S, UCHL1, UCHL3, UQCC1, UTP11L, VSNL1, WDR47, WDR7, WFDC1, XK, YWHAH, ZFP36L1, ZNF365, ZNHIT3

Table 1.

List of 446 significantly differentially expressed genes in Alzheimer’s gene expression datasets.

Note: The differentially expressed genes were detected by the weighted Bayesian random-effect meta-analysis models on the metadata of four Alzheimer’s gene expression datasets: GSE1297, GSE5281, GSE29378, and GSE48350.

We then performed gene network analysis using a publicly available web interface, GeneMANIA [34]. The 446 DE genes identified by the weighted BRE model participate in 146 significant pathways at a false discovery rate of 1%. The first-thirty highly significant pathways with more than twenty identified DE genes in each network included cellular respiration, oxidative phosphorylation, mitochondrial protein complex, inner mitochondrial membrane, protein complex, ATP metabolic process, respiratory electron transport chain, ATP synthesis coupled electron transport, electron transport chain, mitochondrial ATP synthesis coupled, electron transport, mitochondrial inner membrane, energy derivation by oxidation of organic compounds, respiratory chain complex, respirasome NADH dehydrogenase activity, NADH dehydrogenase (quinone) activity, NAD(P)H dehydrogenase (quinone) activity, mitochondrial respirasome, oxidoreductase activity, acting on NAD(P)H, quinone or similar compound as acceptor, respiratory chain complex I, NADH dehydrogenase complex, proton transmembrane transporter activity, aerobic respiration, presynapse, postsynapse, NADH dehydrogenase, complex assembly, oxidoreductase complex, proton-transporting two-sector ATPase complex, mitochondrial proton-transporting ATP synthase complex, ATPase-coupled cation transmembrane transporter activity, synaptic vesicle recycling, inner mitochondrial membrane organization, and cellular response to peptide. GeneMANIA overall retrieved the genes with known coexpression (51.98%), consolidated pathways (25.08%), physical interactions (27.73%), colocalization (10.79%), genetic interactions (5.79%), predicted interactions (2.65%), pathway (0.86%), and share protein domain (0.20%). More details can be found on www.genemania.org.

7. Discussion

In this study, we developed a meta-analytic approach for incorporating sample-quality information into the BRE meta-analysis model using an efficient weight identified by a series of simulation studies [10, 13] to adjust the study heterogeneity in the model. We illustrated the weighted Bayesian approach as compared to the classical RE model through an application in Alzheimer’s gene expression studies. We have seen the results of Alzheimer’s gene expression varied by the sample qualities [13]. The variation of sample quality restricted meta-analysis techniques to properly detect DE genes [35, 36]. Meanwhile, the BRE meta-analysis model allows flexibility in calculating $y_{\mathit{ig}}$ and its variance ${\sigma }_{\mathit{ig}}^2$ as well as study-specific adjustments [37]. We therefore can up-weight good-quality samples and down-weight borderline-quality samples in the model. This developed approach utilizes sample-quality information in the meta-analysis of high-dimensional microarray studies in detecting DE genes.

Additionally, the classical RE model tended to estimate ${\tau }_g^2$ as being zero and the variance of ${\widehat{\beta }}_g$was underestimated, while the BRE meta-analysis model can allow for the uncertainty of the parameter estimates in the model. The BRE model used the marginal posterior distribution of ${\tau }_g^2$ for ${\widehat{\beta }}_g$ estimation, which does not reply on the point estimate of ${\tau }_g^2$. The BRE model can therefore, in turn, increase the fitness of the model [38].

8. Conclusions

This meta-analytic approach with the sample-quality weight can increase the precision and accuracy of the Bayesian random-effects models in gene expression meta-analysis. The performance of the weighted Bayesian random-effects model may be varied depending on data feature, levels of sample quality, and adjustment of parameter estimates.