Predicting gene expression using DNA methylation in three human populations

Datasets

Dataset cleaning and filtering

To assess the DNA methylation effect in prediction gene expression, we defined the “methylation probes” as the 344,303 probes in Table S1 of Grundberg’s (2012) paper. The probes on the methylation array but excluded from Table S1, which have potential SNP effects or cross hybridization effects, are termed “S&C probes”. The combinations of these two types of probes are termed “all probes”. The Adipose dataset has 32,478 missing values in the DNA methylation data. Samples with missing values were excluded from regression analysis of the respective gene. Among the 485,679 probes in the dataset, 344,201 probes remained in Adipose dataset after filtering. For the PBMC dataset, 344,180 out of the 485,461 probes in the dataset remained after filtering and 344,202 out of 485,578 probes remained for the LCL dataset. In order to make the method comparable and the analysis consistent, only genes that have the LASSO models were used, 8040 genes and 149,152 CpG sites in Adipose dataset, 4,252 genes and 73,553 CpG sites in PBMC dataset and 7514 genes and 143,599 CpG sites in LCL dataset (Table S1).

Modelling the relationship between gene expression and DNA methylation

CpG probes were mapped to genes using UCSC RefGene annotation. Gene expression and DNA methylation data for each gene were extracted using in-house perl script. Since there was no missing value for methylation of the PBMC and LCL dataset, all samples were used in the regression analysis.

We used three types of regressions, single linear regression, multiple regression, and least absolute shrinkage and selection operator (LASSO) regression (Tibshirani, 1996), to model the linear relationship between gene expression and DNA methylation. Squared correlation (R²) between predicted and observed data was used to compare the three types of regressions. In the single linear regression, each CpG site was modeled separately to predict gene expression level. The CpG site that provides maximum R² was used to represent each gene. In multiple regressions, all the CpG sites in each gene were used as predictors and the R² was calculated. In the LASSO regression with default parameters, all CpG sites were used to predict the gene expression. We used the GLMNET package in R to fit the LASSO model in which penalized parameters were obtained using 10-fold cross-validation to minimize the mean squared error, while the predictors and response variables were all standardized.

Cross validation

In addition to calculating the R² from fitting the models (fitting R²), we also conducted five-fold cross validation to compare the prediction power of the three regression models using the validation R² (R².cv). Specifically, the samples were randomly separated into training set (4/5 of data) and testing set (1/5 of data). The procedure was iterated 10 times and the mean R² of the 10 five-fold cross validations was used as our final cross validation R² for each model. For single regression, cross validation was conducted for each CpG site and the maximum R² was used for each gene. For LASSO regression analysis, we first obtained the optimal penalty parameter using ten-fold cross validation and then used another five-fold cross validation to evaluate the predictive performance of the model.

Note that we calculated fitting R² in the LASSO cross validation models. We used the entire datasets as testing in the LASSO cross validation models in order to obtain the fitting R² in a fashion consistent with the multiple and single cross-validation models. In this case, all the R² values in the paper are squared correlation of the predicted and the true values in the training set.

Model comparisons on significant genes

We first identified genes that showed overall model prediction p values less than 0.0001 in multiple regressions and then compared the three regression models on these genes.

Gene Ontology (GO) and pathway enrichment analysis

For top 2,000 genes with highest R², we use The Database for Annotation, Visualization and Integrated Discovery (DAVID ) at https://david.ncifcrf.gov/ (Huang, Sherman & Lempicki, 2008) to conduct GO term enrichment analysis based on modified Fisher Exact Test. The background genes were set to be the genes on the expression array, HumanHT-12_V3_0_ R2_11283641_A. The significantly overrepresented GO terms were selected based on the EASE Score, which is the geometric mean of p-values on logarithm scale for the member terms. We applied medium classification stringency in the DAVID website to our data. “GOTERM_BP_FAT” was used to obtain more information in biological processes of the Gene Ontology enrichment analysis. “KEGG_PATHWAY” was selected for pathway enrichment analysis in the same fashion. The most enriched GO terms and pathways with low p-value or FDR were shown in the results.

Gene expression prediction using different type of probes on the methylation microarray

The probes excluded by Table S1 of Grundberg’s (2012) paper were treated as probes with SNP and/or hybridization effects (S&C probes). We compared these probes, the methylation probes, and the combination of these two types of probes in predicting gene expression.

Analysis codes

We wrapped up our major analysis codes into a package at https://github.com/dorothyzh/MethylXcan. It includes all three regressions and calculates the squared correlation for each model. The package is written in R and Perl, and has been tested under linux or MACSOX system. Users can use this package on the datasets described here or on their own data after formatting their methylation data, expression profiling data, and annotation files as specified by the package.

Multiple regressions using all methylation CpGs from a gene predict gene expression the best in model fitting

As a reference, we first conducted association analysis on each methylation probe in predicting gene expression using single regression. Most of the genes with valid LASSO models have at least one significant CpG at nominal significance level of 0.05. For example, in the Adipose dataset, 7,326 out of 8,040 genes have at least one CpG site significant at 0.05 level and 3460 out of 8040 genes have at least one CpG site significant at 0.0001. However, the prediction power represented by the largest R² in each gene is generally low. Only 19 genes have R² greater than 0.3 and 486 genes have R² larger than 0.1 when the most predictive CpG site is considered (Table 1). Similar results were obtained from the PBMC and LCL datasets, except that the PBMC dataset has a substantially smaller number of genes with a CpG significant at 0.0001 level (582 out of 4,252 genes) although the distribution of the estimated R² is similar to that from the Adipose dataset. This could be due to the smaller sample size or the nature of the PBMC tissue source. On the other hand, the LCL dataset has a larger number of genes with higher R² from single CpG regression analysis, which could be related to the homogeneous nature of cell lines.

Since multiple CpG sites from each gene were assayed on the methylation microarray, we applied multiple linear regression to utilize all methylation CpG sites in the gene region as predictors simultaneously. The R² explained by the regression model did improve substantially for the majority of genes compared with that from the single linear regression (Fig. 1). As expected, the significant genes from multiple regression analyses tend to have larger R² compared with the non-significant genes. The improvement of R² from the multiple regression over single regress in the PBMC and LCL datasets is similar to that in the Adipose dataset (Fig. 1).

Compared with multiple regressions, LASSO regression did not generate R² quite as high in all datasets (Figs. 1D–1F), which is also indicated by smaller number of genes with R² exceeding each threshold (Table 1). For example, the number of genes with R² greater than 0.2 decreased from 476 to 360 for the Adipose dataset, from 1,163 to 561 for the PBMC datasets, and from 2,170 to 1,704 for the LCL dataset. Similar trend was observed at the other two thresholds for all three datasets.

It is widely recognized that gene expression is negatively correlated with DNA methylation level in the promoter region but often positively correlated with DNA methylation level in gene body (Ball et al., 2009; Wu et al., 2010; Jones, 2012). The different directions of correlation among CpG sites in the same gene may lead to perceptions that combining all CpG sites is not advantageous in the prediction of gene expression. However, the multiple and lasso regressions can accommodate coefficients in different directions without affecting prediction power. Nevertheless, we tested CpG sites from promoter region and those from gene body for prediction separately in the LCL dataset. As expected, neither performs as well as combined (Fig. S1).

LASSO regression shows better prediction in cross-validation

To better assess the accuracy of the predictive models, we performed 5-fold cross validation on single, multiple regressions, and LASSO regressions to estimate the prediction R². The results showed that the LASSO regression produced much larger R² values than the single regression and less dramatic but discernible increases over multiple regressions (Fig. 2). These differences are also reflected in the number of genes with R² exceeding the three thresholds. For example, 827 genes (10.29%) from the Adipose dataset have R² greater than 0.1 from LASSO regression, while 722 genes (8.98%) and 106 genes (1.32%) have R² greater than 0.1 from multiple regression and single regression, respectively (Table 1). For genes with R² greater than 0.3, LASSO regression has 42 genes (0.52%) while multiple regression and single regression have 38 (0.47%) and 2 genes (0.02%), respectively. These results indicate that penalized regression has better prediction than multiple or single regressions in cross-validation. Cross validation tends to overcome bias and over-fitting issues. As expected, cross-validation R² values are generally lower than those from the model fittings, which is reflected by the smaller number of genes with R² values greater than the R² thresholds (Table 1). Similar results were obtained from the PBMC dataset and the LCL dataset.

To make sure that the prediction R² is larger than thoses from random chance, we compared the cumulative R² from the three datasets with those from the the null distribution of correlations based on Fisher z-transfromation in quantile–quantile plots (Fig. 3). All datasets showed that the observed R² values are much larger than the expected R² values from random chance. In addition, the departure is the largest in the LCL dataset followed by the Adipose dataset and the PBMC dataset when methylation probes were considered, indicating that the LASSO models capture a larger proportion of the transcriptome variability in the LCL dataset than in the other two datasets. This is potentially due to the combination of sample size and nature of different tissues.

To rule out the possibility that prediction R² is mainly driven by the variability of gene expression and the variability of DNA methylation across individuals in the study population, we first examined the correlation between the variability of gene expression with R² from LASSO regression. No obvious correlation was observed (Fig. S2). For assessing the correlation between DNA methylation variability and prediction R², we took the CpG with maximum R² from single regression and examined the correlation between its variability with the prediction R² from LASSO regression. We only observed a potential positive correlation in the Adipose dataset when the R² is greater than 0.5, where there are a small number of genes (Fig. S3).

Using all probes improves prediction power for gene expression

In order to evaluate DNA methylation power in predicting gene expression, we first left out a large proportion of probes potentially affected by genetic or cross-hybridization effects (Table S1). However, using all probes on the array is preferred if our goal is to achieve better prediction accuracy of gene expression. To evaluate the prediction power from all probes, we included all available probes in LASSO regression and found that the overall prediction power did increase compared to the models using only the methylation probes (Fig. 3). We observed more genes with R² values exceeding the thresholds (Table 1). In addition, the largest R² value is much larger when all probes are used compared to that from only the methylation probes. For example, the largest R² is 0.92 from all probes compared to 0.74 from only methylation probes in the Adipose dataset. Similarly, the largest R² increases from 0.71 to 0.88 in the PBMC dataset and from 0.76 to 0.87 in the LCL dataset. Furthermore, valid LASSO models are available for more genes when all probes are used (Fig. S4).

The increase of prediction power on gene expression from all types of probes on the methylation microarray suggests that there is contribution from the probes with potential SNP effects or cross hybridization effects. To further assess the size and nature of their contribution, we separately estimated the prediction power of the methylation probes, S&C probes, and the combination of them (all probes). The results showed that the S&C probes have independent prediction power from the methylation probes and the combination of both has increased prediction power over the methylation probes alone (red points vs black line in Figs. 4A, 4C and 4E). The prediction power from the S&C probes was also estimated for genes with enough SNP probes to form a LASSO model and their prediction power are mostly above zero (blue points in Fig. 4). The fact that the blue points are randomly distributed instead of following the black line suggests that the two sources of R² are not correlated; therefore, the genetic effect and epigenetic effect do not seem to coexist in the same genes. Figure 5 shows some examples of genes with large prediction powers from either methylation probes or S&C probes. As expected, the methylation probes tend to show continuous methylation values while the S&C probes tend to show categorical values due to limited genotypes of the samples.

GO term analysis of better predicted genes

To examine the potential biological function of the genes showing relatively higher predictability, we conducted gene ontology (GO) enrichment analysis using DAVID on genes with R² larger than 0.2 from LASSO regression of methylation CpGs. At false discovery rate (FDR) of 0.01, cell adhesion, lipids metabolism, and regulation of immune system are among the most significantly enriched terms in the Adipose dataset (Table S2), which seem to be consistent with the previous findings for subcutaneous fat cells (Berg & Scherer, 2005). For the PBMC dataset, the most significantly enriched terms are mostly related to defense and immune functions, lymphocyte aggregation, T cell activation, inflammatory response, as well as cell adhesion. These results appear to be reasonable for atopy and persistent asthma blood cells. For the LCL data, some terms related to cell adhesion, migration, communication, and morphogenesis are highly enriched. Same GO term analyses were also conducted for R² from all CpG probes and similar results were obtained.