Methylation Normalization ← Return to Omics

Single variable description

Descriptive statistics is provided for each variable in the metadata. The summary can be downloaded as an Excel file using this link. Here is how the summary statistics were calculated:

• For categorical variables we provide the number of categories. We aslo calculate the minimum, mean, median and the maximum values for the distribution of the cases among the categories and the number of the missing cases.
• For numeric variables we provide minimum, mean, median, maximum values and the number of missing cases. In the column with the number of categories numeric variables have blank.

Relationship between the variables in the metadata

To understand the relationship between the variables (describing the subjects and the experiments) we performed pairwise tests:

• between numeric variables - Pearson correlation
• between categorical variables - Chi-square test using Monte Carlo simulation with with 2000 replications
• between numeric and categorical variables - Kruskal-Wallis rank sum test

For each test we recorded the P value, the test statistic and the type of the test. If a test produced an error or a warning, we recorded the message and the type of the test. No tests were performed if two variables were the same (tested against itself). Prior to calculations we removed all sample identification variables (such as projid, order_id and Sample_ID, Sample) as we discovered that they significantly slow down the calculations (especially for the Monte Carlo simulations). If the variables were of any other type combination rather than numeric-numeric, numeric-categorical, categorical-numeric or categorical-categorical no test was performed. We didn't perform any multiple testing adjustment. Results can be downloaded as an Excel file using this link.

For visualization purposes we divided the P values into significant (less or equal than 0.05) and non-significant (larger than 0.05) and plotted the results as a heat map below. Gray squares indicate tests against itself (main diagonal) and tests which produced errors or warnings. Refer to the Excel file for the explanation of errors and warnings.

Relationship between the batch and variables describing arrays

Illumina slides (Sentrix_ID) generally contain 6, 8 or 12 arrays (Sentrix_Position). Prior to hybridization DNA samples are PCR amplified. It is common for amplification or hybridization to occur in batches because of the number of samples. In this section we investigate the relationship between slides, arrays, batches and amplification plates/rows/wells.

Figure 1: Slides were hybridized in two batches

Each slide (Sentrix ID) contains 12 arrays (Sentrix Position), some arrays are missing. There are total of 65 slides used for this project. Arrays on each slide are grouped by batch.

Figures 2 and 3: Samples were amplified in two batches

Figure on the left shows the relationship between the amplification sample plate (Sample_Plate, x axis) and the Sentrix ID. This shows that the plates were amplified in two batches and the samples from the plates were distributed across multiple slides. The samples from the same amplification batch were not randomly distributed across the slides (all arrays on the same slide come from the same batch).

The figure on the right shows that the problem of array hybridization from the same row of an amplification plate was avoided. Each slide received samples from different rows on the plates.

Explore the effect of the variables on DNA methylation data

We used Singular Value Decomposition (SVD) to understand the effect of the variables on the intensity normalized data (using snm package). Figure on the left below shows the relative importance of every eigen vector. Figure on the right shows the relationship between two arrays which correspond to the most extreme values of the first eigen vector.

Next we take the first 5 largest eigen vectors and calculate their relationship with the clinical and the technical variables (with the exception of projid, order_id and Sample_ID). For each eigen vector we fitted a linear model with a clinical or a technical variable as a dependent variable. Here we report the overall P values for ANOVA F-tests for each model. Table below provides the P values while the figure highlights the most significant relationships. No multiple testing adjustment was performed.

We also plotted the relationship of some of the most influential clinical and technical variables with the eigen vectors. The scatter plot of the eigenvector 1 vs eigenvector 2 colored by the batch shows a distinct separation between the arrays. The effect of the Sentrix ID and the Sample Plate are similar because these variables are highly correlated. However their effect is most significant in the second eigenvector rather than the first eigenvector. The first eigenvector shows a strong correlation with NNLS (as well as CETS, although we are not showing it here because of the missing values present in CETS).

Remove subjects with missing values in clinical variables

We are interested in the relation of DNA methylation to the following variables: cogn_global_lb, amyloid and braaks. We also want to explicitly remove the effects of CETS, pmi, sex and age at death. Variables CETS, pmi, cogn_global_lv and amyloid have missing values for some of the subjects. First, we removed all subjects with missing values for these variables, the remaining number of subjects is 702.

Normalization using supervised normalization of microarrays (snm)

SNM approach identifies two "buckets" for data normalization: biological (variables of interest) and adjustment (variables which effect is to be removed). In this case the biological "bucket" is represented by a the following model: ~ amyloid + braaksc + cogn_global_lv . The adjustment "bucket" is represented by the following model: ~ Sentrix_ID + CETS + pmi + sex + age_death . At the same time we also perform intensity normalization of the data.

The figures below demonstrate percent variance explained by the eigen vectors after applying snm (left) and the relationship between the arrays corresponding to the most extreme points of the first eigen vector (right).

Relation of the first 5 eigen vectors with the variables after normalization can be downloaded in an Excel format. The plot below highlights the most significant variables. We observe some correlation (P value 0.02) of the Sample Plate and Sample Row with the second eigen vector (~ 2.5% variance explained).

The figures below show the relationship of the first eigen vectors with the technical variables after fitting the model.

Normalization using snm with extended adjustment model

Added on February 4, 2016: Dr Chris Gaiteri performed initial network analysis on the data normalized using the model above (after selecting CpGs near genes and then collapsing these selected CpGs into regions using minfi::cpGCollapse function). Using averaged methylation values extracted from the identified networks we showed that these vectors are correlated with the Sentrix Position (figure on the left below). We then regressed out the effect of the Sentrix Position from each of the vectors and calculated correlation between the variables (figure on the right).

Based on these results we included Sentrix Position in the adjustment model. Figures below show the percent variance explained and the outliers after normalization. These look very similar to the normalization results above.

Relation of the first 5 eigen vectors with the variables after the normalization can be downloaded in an Excel format. Including Sentrix Position removed all residual correlation with it and other technical variables.

The figures below show the relationship of the first eigen vectors with the technical variables after fitting the model.

These normalized data were then again subjected to filtering based on proximity to a known gene and collapsed to regions. Dr Chris Gaiteri built a network and extracted average methylation vectors for which we again calculated correlation with the metadata. While correlation with the Sentrix Position has completely removed we observed a strong correlation with the NNLS variable. The figure below demonstrates these correlations.

We therefore decided to include NNLS as another adjustment variable in the normalization model. The section below (added on February 16, 2016) explores the normalized data.

Include NNLS as an adjustment variable

While the model with biological variables remains the same, the adjustment model includes NNLS, CETS, Sentrix ID, Sentrix Position, pmi, sex and age at death.

Relation of the first 5 eigen vectors with the variables after normalization can be downloaded in an Excel format.

The figures below show the relationship of the first eigen vectors with the technical variables after fitting the model.

These normalized data were then again subjected to filtering based on proximity to a known gene and collapsed to regions. Dr Chris Gaiteri built a network and extracted average methylation vectors for which we calculated correlation with the metadata. The figure below shows demonstrates these correlations.