Epistasis Test in Meta-Analysis: A Multi-Parameter Markov Chain Monte Carlo Model for Consistency of Evidence

<div><p>Conventional genome-wide association studies (GWAS) have been proven to be a successful strategy for identifying genetic variants associated with complex human traits. However, there is still a large heritability gap between GWAS and transitional family studies. The “missing heritability” has been suggested to be due to lack of studies focused on epistasis, also called gene–gene interactions, because individual trials have often had insufficient sample size. Meta-analysis is a common method for increasing statistical power. However, sufficient detailed information is difficult to obtain. A previous study employed a meta-regression-based method to detect epistasis, but it faced the challenge of inconsistent estimates. Here, we describe a Markov chain Monte Carlo-based method, called “Epistasis Test in Meta-Analysis” (ETMA), which uses genotype summary data to obtain consistent estimates of epistasis effects in meta-analysis. We defined a series of conditions to generate simulation data and tested the power and type I error rates in ETMA, individual data analysis and conventional meta-regression-based method. ETMA not only successfully facilitated consistency of evidence but also yielded acceptable type I error and higher power than conventional meta-regression. We applied ETMA to three real meta-analysis data sets. We found significant gene–gene interactions in the renin–angiotensin system and the polycyclic aromatic hydrocarbon metabolism pathway, with strong supporting evidence. In addition, glutathione <i>S</i>-transferase (GST) mu 1 and theta 1 were confirmed to exert independent effects on cancer. We concluded that the application of ETMA to real meta-analysis data was successful. Finally, we developed an R package, etma, for the detection of epistasis in meta-analysis [etma is available via the Comprehensive R Archive Network (CRAN) at <a href="https://cran.r-project.org/web/packages/etma/index.html" target="_blank">https://cran.r-project.org/web/packages/etma/index.html</a>].</p></div

Inconsistent estimates of interaction effects in the same data.

<p>This figure describes a meta-regression analysis based on the data from Fang et al. [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0152891#pone.0152891.ref027" target="_blank">27</a>] (detailed data are shown in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0152891#pone.0152891.s004" target="_blank">S1 Table</a>). The upper plot describes an investigation of the association between proportions of null/null GSTT1 in cases and the odds ratios of GSTM1 in cancer, and the lower plot describes an investigation of the association between proportions of null/null GSTM1 in cases and the odds ratios of GSTT1 in cancer. The solid lines denote unbiased estimators of odds ratios, and the dashed lines show 95% confidence intervals of odds ratios. According to a previous article, the slopes in meta regression approximate interaction effects [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0152891#pone.0152891.ref016" target="_blank">16</a>]. However, the estimates of interaction effect were inconsistent when we exchanged the independent and moderator variables (0.1377 and 0.2338, respectively). This phenomenon does not occur in individual data analysis and leads to problems in interpretation.</p

The result of real data analysis using ETMA.

<p>The result of real data analysis using ETMA.</p

Summary of simulation conditions.

<p>Summary of simulation conditions.</p

Type I error of individual data analysis, ETMA and conventional meta-regression.

<p>Type I error of individual data analysis, ETMA and conventional meta-regression.</p

A typical analysis pipeline of ETMA function in 'etma' package.

<p>This figure summarized the pipeline of ETMA function. The main input is a meta-analysis dataset, which including the number of wild/mutation type of SNP1/SNP2 in case/control group. The main options include the length of chains in step 1/2, the maximum number of iterations, and the start seed. Main outputs include three matrixes. Matrix b includes the beta values (logarithmic ORs) of each SNP and interaction term, and VCOV is the variance covariance matrix of beta value. <i>P</i> is an n by 3 matrix describing three study-specific parameters (p1 = Disease risk in subjects with wild-type alleles of SNP1 and SNP2; p5 = Mutation frequency of SNP1; p6 = Mutation frequency of SNP2)</p

The proportion of individual with different status of disease/SNP1/SNP2 could be calculated by <i>p</i>_baseline, <i>MAF</i>₁, <i>MAF</i>₂, OR_y,SNP1, OR_y,SNP2 and OR_interaction.

<p>The proportion of individual with different status of disease/SNP1/SNP2 could be calculated by <i>p</i>_baseline, <i>MAF</i>₁, <i>MAF</i>₂, OR_y,SNP1, OR_y,SNP2 and OR_interaction.</p

The statistical power of individual data analysis, ETMA and conventional meta-regression.

<p>The <i>x</i>-axis describes three levels of interaction effect (OR_interaction = 1.2, 1.5 or 2.0), and the <i>y</i>-axis indicates the statistical power provided by individual data analysis (black), ETMA (red) and conventional meta-regression (blue), respectively. The details of these methods are described in the Method. The different subplots present comparisons using different simulation parameters, and the titles of these subplots show their detailed settings. Each data point was based on 1,000 simulations.</p

Prediction of novel target genes and pathways involved in irinotecan-resistant colorectal cancer

<div><p>Background</p><p>Acquired drug resistance to the chemotherapeutic drug irinotecan (the active metabolite of which is SN-38) is one of the significant obstacles in the treatment of advanced colorectal cancer (CRC). The molecular mechanism or targets mediating irinotecan resistance are still unclear. It is urgent to find the irinotecan response biomarkers to improve CRC patients’ therapy.</p><p>Methods</p><p>Genetic Omnibus Database GSE42387 which contained the gene expression profiles of parental and irinotecan-resistant HCT-116 cell lines was used. Differentially expressed genes (DEGs) between parental and irinotecan-resistant cells, protein-protein interactions (PPIs), gene ontologies (GOs) and pathway analysis were performed to identify the overall biological changes. The most common DEGs in the PPIs, GOs and pathways were identified and were validated clinically by their ability to predict overall survival and disease free survival. The gene-gene expression correlation and gene-resistance correlation was also evaluated in CRC patients using The Cancer Genomic Atlas data (TCGA).</p><p>Results</p><p>The 135 DEGs were identified of which 36 were upregulated and 99 were down regulated. After mapping the PPI networks, the GOs and the pathways, nine genes (GNAS, PRKACB, MECOM, PLA2G4C, BMP6, BDNF, DLG4, FGF2 and FGF9) were found to be commonly enriched. Signal transduction was the most significant GO and MAPK pathway was the most significant pathway. The five genes (FGF2, FGF9, PRKACB, MECOM and PLA2G4C) in the MAPK pathway were all contained in the signal transduction and the levels of those genes were upregulated. The FGF2, FGF9 and MECOM expression were highly associated with CRC patients’ survival rate but not PRKACB and PLA2G4C. In addition, FGF9 was also associated with irinotecan resistance and poor disease free survival. FGF2, FGF9 and PRKACB were positively correlated with each other while MECOM correlated positively with FGF9 and PLA2G4C, and correlated negatively with FGF2 and PRKACB after doing gene-gene expression correlation.</p><p>Conclusion</p><p>Targeting the MAPK signal transduction pathway through the targeting of the FGF2, FGF9, MECOM, PLA2G4C and PRKACB might increase tumor responsiveness to irinotecan treatment.</p></div

Protein–protein interaction (PPI) network of differentially expressed genes (A) upregulated genes and (B) downregulated genes.

<p>The PPI pairs were imported into Cytoscape software as described in Methods and Materials. Pink nodes represent up regulated genes while green nodes represent down regulated genes. The lines represent interaction relationship between nodes.</p