Computational protocol: Cell-composition effects in the analysis of DNA methylation array data: a mathematical perspective

Similar protocols

Protocol publication

[…] To evaluate and demonstrate the concepts proposed above, we applied the proposed method to nine DNA methylation data sets obtained from GEO, each described in Table . Table  also indicates he regression model used in each analysis. For all but two female-only data sets, sex chromosome loci were omitted. For each of the 450K data sets, 166,314 CpGs were removed due to known problems with cross-reactivity or polymorphisms [] and loci with greater than 5 missing values were also omitted. For 27K data sets we considered a range of k from 1 to 25; for three of the 450K data sets we considered a range of k from 1 to 50, but due to small residual degrees of freedom, we considered only a range up to 30 for the Reinius 450K blood reference data set.For each data set and for each value of k, we computed B k, Δ k, and their element-wise bootstrap standard errors (100 bootstrap samples for each analysis) using our previously published reference-free algorithm []. From estimates and standard errors we computed p-values, and across each column of B k and of Δ k we computed the corresponding q-values using the Bioconductor qvalue package (version 1.34.0). Additionally, we applied the two new proposed methods for estimating k, as well as two variants of the previously proposed random matrix theory approach, one that scales the rows of R, and one that does not.To assess the biological significance of results, we selected several gene sets and investigated the enrichment of significant (q < 0.05) loci within each set. For each data set, we used the value of k selected by our proposed method, except for the artery data set where we used k = 10 due difficulty in estimating k, caused by small sample size. The first gene set consisted of known DMRs for leukocytes: for 27K, the 500 CpG sites published for 27K in our earlier paper [], and for 450K, 417 CpGs not excluded from the top 600 CpGs reported by Jaffe & Irrizary [] based on the Reinius data set []. The second was a set of CpGs mapped to polycomb target genes, compiled from four separate published articles [-] and used extensively in our previous work. The remaining gene sets were selected KEGG pathways, obtained via Bioconductor annotation packages for the 27K and 450K platforms. Please see Table  for a summary of pathways investigated. To test over-representation while circumventing known problems with the application of such gene-set analysis to DNA methylation data, we used the exact Mantel-Haenzel test to stratify by genomic context. For 27K we stratified by CpG Island status, and for 450K we stratified by strata defined by Infinium chemistry type, relation to CpG Island, and gene region. Since a single CpG may be mapped to different region designations due to splice variants or gene adjacency, we used the following rule to establish precedence: promoter (“TSS1500” or “TSS200”) were combined as “TSS” and had highest precedence, and the remaining chain of precedence was as follows; TSS > 1stExon > Body > 5’UTR > 3’UTR > null. In cases where sparsity prevented stratification (a few analyses with PcG gene set) we used the unstratified Fisher test. We also evaluated the behavior of the reference-free method in null scenarios. Although our original paper conducted simulations using artificial data sets with only 1000 CpGs, we devised simple but more realistic approach based on real 27K data sets. From the HNSCC and ovarian cancer case/control data sets, we obtained three null data sets as follows. First, to simulate a completely null result, we permuted the phenotype information (case status and age) with respect to the array data. Second, to simulate a null linear result with some non-null, nonlinear effects, we selected CpGs with q < 0.0001 for the case coefficient in a limma analysis, unadjusted for cell-type, conducted on an M-value (logit-beta) scale (458 CpGs for the HNSCC data set and 817 for the ovarian cancer data set), we then multiplied the corresponding coefficients on M-values by random scalars drawn from a normal distribution with standard deviation 2, and finally applied the corresponding effects to M-values whose case effect had been removed (on a beta scale so that any linear age effect would be maintained), converting the result back to beta scale. Note that this approach would produce a data set with weak linear age effects, no linear case effects, and nonlinear case effects at a small fraction of CpGs. The third approach was identical to the second except that a threshold of q < 0.001 was used (resulting in 643 CpGs for the HNSCC data set and 1163 for the ovarian cancer data set). Each of these three approaches was applied five separate times, with subsequent analysis by the reference-free method across values of k ranging from 1 to 25, and tabulation of the resulting number of significant (q < 0.05) coefficients. Our hypothesis was that the first set of analyses would produce very few significant q-values for either Δ k or B k, and that the latter two sets would produce some significant q-values for the B k coefficients but few to none for the Δ k coefficients.In order to illustrate the principle of recursive partitioning above, we conducted an additional analysis on the two reference data sets. For each value of k, we applied hierarchical clustering to ΔkT (i.e. clustered the columns of Δ k), using Euclidean metric and Ward’s method of clustering. We have demonstrated that Ward’s method performs similarly to the recursively partitioned mixture model []. We hypothesized that even with relatively small values of k, phylogenetic relationships between cell types will be evident in the resulting clusterings.Finally, as an evaluation of the hypothesis that the left-singular vectors u k of the SVD in (4) concentrate a majority of information about major cell types within the initial values of k, we computed the cross products of the left-singular vectors u k across all five 27K data sets and across all four 450K data sets. We also clustered the intercept columns of each Δ k across 27K data sets and across 450K data sets, and computed correlation coefficients stratified over CpG Island status (27K) or more general genomic context (450K) as in the aforementioned Mantel-Haenzel test. We performed a similar analysis for the B k coefficients.The core engine of the methods discussed in this article and published in our previous work are available in the CRAN/R package RefFreeEWAS. The novel elements proposed in this article are available in the Additional file . […]

Pipeline specifications

Software tools limma, RefFreeEWAS
Application Phylogenetics