Expression and Methylation Analysis with MEAL
23 May 2018
Package
MEAL 1.10.1
1 Introduction
Experiments comprising different omic data analyses are becoming common. In this cases, the researchers are interested in analyzing the different omic types independently as well as finding association between the different layers.
In this vignette, we present how address to an analysis of methylation and expression data using MEAL. We will show how to analyze each data type independently and how to integrate the results of both analyses. We will also exemplify how to run an analysis when we are interested in our target region.
We will use the data from brgedata package. This package contains data from a spanish children cohort. We will work with the methylation and the expression data encapsulated in a GenomicRatioSet and in an ExpressionSet respectively. In our analysis, we will evaluate the effect of sex on methylation and expression. We will also deeply explore changes in a region of chromosome 11 around MMP3, a gene differently expressed in men and women (chr11:102600000-103300000).
Let’s start by loading the required packages and data:
library(MEAL)
library(brgedata)
library(MultiDataSet)
library(missMethyl)
library(minfi)
library(GenomicRanges)
library(ggplot2)From the six objects of brgedata, we will focus on those containing methylation and expression data: brge_methy and brge_gexp.
brge_methy is a GenomicRatioSet with methylation data corresponding to the Illumina Human Methylation 450K. It contains 392277 probes and 20 samples. The object contains 9 phenotypic variables (age, sex and cell types proportions):
data(brge_methy)
brge_methy## class: GenomicRatioSet
## dim: 392277 20
## metadata(0):
## assays(1): Beta
## rownames(392277): cg13869341 cg24669183 ... cg26251715 cg25640065
## rowData names(14): Forward_Sequence SourceSeq ...
## Regulatory_Feature_Group DHS
## colnames(20): x0017 x0043 ... x0077 x0079
## colData names(9): age sex ... Mono Neu
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19
## Preprocessing
## Method: NA
## minfi version: NA
## Manifest version: NAcolData(brge_methy)## DataFrame with 20 rows and 9 columns
## age sex NK Bcell
## <numeric> <factor> <numeric> <numeric>
## x0017 4 Female -5.81386004885394e-19 0.173507848676514
## x0043 4 Female 0.00164826459258012 0.182017163135389
## x0036 4 Male 0.0113186314286605 0.169017297114852
## x0041 4 Female 0.00850822150291592 0.0697157716842748
## x0032 4 Male 0 0.113977975025216
## ... ... ... ... ...
## x0018 4 Female 0.0170284419499822 0.0781697956157343
## x0057 4 Female 0 0.0797774108319114
## x0061 4 Female 0 0.164026603367578
## x0077 4 Female 0 0.11227306633587
## x0079 4 Female 0.0120147521727839 0.091365002951843
## CD4T CD8T Eos
## <numeric> <numeric> <numeric>
## x0017 0.204068686112657 0.100974056859719 0
## x0043 0.161712715578477 0.128772235710367 0
## x0036 0.155463683913806 0.12774174901185 0
## x0041 0.00732789124320035 0.0321860552286678 0
## x0032 0.222303988949171 0.0216089881163733 2.60503891663895e-18
## ... ... ... ...
## x0018 0.112735084403202 0.0679681631510725 0.00837769947260005
## x0057 0.111071701794334 0.00910488689520549 4.61887343424789e-18
## x0061 0.224202761557123 0.131252122776747 0.0085107393064941
## x0077 0.168056145550346 0.078405927598852 0.0613396950863682
## x0079 0.205830330328636 0.114753887947108 0
## Mono Neu
## <numeric> <numeric>
## x0017 0.0385653506040031 0.490936252606603
## x0043 0.049954190999866 0.491822214741198
## x0036 0.102710482176672 0.459631855532962
## x0041 0.0718279616580959 0.807749226994262
## x0032 0.0567246148936481 0.575614453750378
## ... ... ...
## x0018 0.0579534708330382 0.658600075099212
## x0057 0.101526022615966 0.686010320431362
## x0061 0.0382646576955094 0.429574723586544
## x0077 0.0583411264640133 0.515284487718098
## x0079 0.0750535276475847 0.513596840714489brge_gexp is an ExpressionSet with the expression data corresponding to an Affimetrix GeneChip Human Gene 2.0 ST Array. It contains 67528 features and 100 samples as well as two phenotypic variables (age and sex):
data(brge_gexp)
brge_gexp## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 67528 features, 100 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: x0001 x0002 ... x0139 (100 total)
## varLabels: age sex
## varMetadata: labelDescription
## featureData
## featureNames: TC01000001.hg.1 TC01000002.hg.1 ...
## TCUn_gl000247000001.hg.1 (67528 total)
## fvarLabels: transcript_cluster_id probeset_id ... notes (11 total)
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:lapply(pData(brge_gexp), table)## $age
##
## 4
## 99
##
## $sex
##
## Female Male
## 48 51Finally, we will create a GenomicRanges with our target region:
targetRange <- GRanges("chr11:102600000-103300000")Let’s illustrate the use of this package by analyzing the effect of the sex in methylation and expression. First, a genome wide analysis will be performed and then the region funcionalities will be introduced. It should be noticed that methylation and expression analyses will include hypothesis testing and visualitzation.
2 Methylation Analysis
2.1 Running the analyses
This demonstration will show the main functions needed to perform a methylation analysis. A more exhaustive description of these and other auxiliary functions can be found at MEAL vignette.
The function runPipeline can run different analyses to the same dataset: limma, DiffVar, bumphunter, blockFinder, DMRcate and RDA. This function accepts a GenomicRatioSet, an ExpressionSet or a SummarizedExperiment, so it can analyze methylation and gene expression data. We can pass our variables of interest to function runPipeline with the parameter variable_names. runPipeline also allows the inclusion of adjustment variables (such as cell count) with the parameter covariable_names. When we are interested in a target region, we can apply RDA by passing a GenomicRanges to the argument range. Users can find a deeper explanation of this function in the Methylation Analysis with MEAL vignette.
As we are interested in evaluating the effect of sex, we will set variable_names to “sex”. We will include cell counts as covaribles to adjust their effect and we will pass our target region to range to compute RDA:
cellCounts <- colnames(colData(brge_methy))[3:9]
methRes <- runPipeline(set = brge_methy, variable_names = "sex", covariable_names = cellCounts, range = targetRange)
methRes## Object of class 'ResultSet'
## . created with: runPipeline
## . sva: no
## . #results: 6 ( error: 1 )
## . featureData: 392277 probes x 19 variablesnames(methRes)## [1] "DiffMean" "DiffVar" "bumphunter" "blockFinder" "dmrcate"
## [6] "RDA"The analysis generates a ResultSet object containing the results of six different analyses: DiffMean, DiffVar, bumphunter, blockFinder, DMRcate and RDA. The six methods are based in the same linear model that we specified with the parameters variable_names and covariable_names. Consequently, the first step after running the pipeline will be evaluate how this model fits the data.
We will rely on the results of limma (DiffMean), as this method run the same linear model to all CpGs. A common way to evaluate the goodness of fit of these models is by making a QQplot on the p-values of the linear regression. This plot compares the distribution of p-values obtained in the analysis with a theoretical distribution. Our plot also shows the lambda, a measure of inflation. If lambda is higher than 1, the p-values are smaller than expected. If it is lower than 1, p-values are bigger than expected. In both causes, the most common cause is that we should include more variables in the model.
We can get a QQplot with the function plot. When applied to a ResultSet, we can make different plots depending on the method. We can select the name of the result with the parameter rid and the kind of plot with the parameter type. We will set rid to “DiffMean” to use limma results and type to “qq” to get a QQplot:
plot(methRes, rid = "DiffMean", type = "qq")
When there are few CpGs modicated, we expect points in the QQplot lying on the theoretical line. However, as all the CpGs from the sexual chromosomes will be differentially methylated between males and females, we obtained a S-shape. In this case, lambda is lower than 1 (0.851).
There are different ways to address this issue. An option is to include other covariates in our model. When we do not have more covariates, as it is our case, we can apply Surrogate Variable Analysis (SVA) to the data. SVA is a statistical technique that tries to determine hidden covariates or SV (Surrogate Variables) based on the measurements. This method is very useful to correct inflation when we do not have the variables that are missing in the model.
Our function runPipeline includes the parameter sva that runs SVA and includes these variables as covariates in the linear model. We do not need to include cell counts as covariates as their effect is be included in the SVs :
methRes <- runPipeline(set = brge_methy, variable_names = "sex", range = targetRange, sva = TRUE)
methResWe have included two methods to analyze individual CpGs: DiffMean (limma) and DiffVar. DiffMean tests, for each CpG, if there is a difference between the means of the different groups. In our case, we tested, for each CpG, if the mean methylation in boys was different to the mean methylation in girls. On the other hand, DiffVar tests, for each CpG, if there is a difference in variance between the groups. In our case, we tested for each CpG if the variance was different in boys than in girls.
As both methods return results for each CpG, we can get an overview of the results using a Manhattan plot. A Manhattan plot represents the p-value of each CpG versus their location in the genome. This representation is useful to find genomic regions exhibiting differentiated behaviours.
We can get a Manhattan plot from a ResultSet using the function plot and setting type to “manhattan”. We will make two plots, one for DiffMean results and the other for DiffVar. We will highlight the CpGs of our target region by passing a GenomicRanges to highlight argument. GenomicRanges passed to plot should have chromosomes coded as number:
targetRangeNum <- GRanges("11:102600000-103300000")
plot(methRes, rid = "DiffMean", type = "manhattan", main = "Differences in Means", highlight = targetRangeNum)
plot(methRes, rid = "DiffVar", type = "manhattan", main = "Differences in Variances", highlight = targetRangeNum)
There are no genomic regions with a lot of CpGs with differences in means. Our target region does not seem to be very different from the rest. On the other hand, chromosome X has a lot of CpGs with big differences in variance. Our target region does not look different.
2.2 Exploring a target region
There are two options in MEAL to explore a target region. If we run RDA in our pipeline, we can plot the RDA components with the function plotRDA. We can color the samples by a given phenotype passing a data.frame to pheno argument:
plotRDA(object = methRes,
pheno = colData(brge_methy)[, "sex", drop = FALSE])
We do not see a different distribution of female and male samples in the first RDA component. Our model explained very few variability (R2 = 0.064) and the p-value is big, suggesting that the variance explained by our model is not bigger than what we can expect by chance.
We can get more estimates of the RDA with the function getRDAresults:
getRDAresults(methRes)## R2 pval global.R2 global.pval
## 0.06568320 0.08439156 0.12967579 0.03079692R2 and pval are the same values that are printed in the RDA plot. global.R2 is the genome wide variance explained by the model. global.pval is the probability of finding a region in the genome where our model explains more variability than in our target region.
If we did not have a priori a target region and we did not run RDA in the pipeline, we have another option. We can plot the results of the other methods (DiffMean, DiffVar, bumphunter, blockFinder and DMRcate) using the function plotRegion. We should pass a GenomicRanges with our target region to the argument range:
plotRegion(rset = methRes, range = targetRange)
This plot has three parts. At the top, we can find the transcripts and CpGs annotations. In the middle, the results of DMR detection methods. In the bottom, the results of DiffMean (blue) and DiffVar (green). For both results, there are the coefficients (lines) and the p-values (dots). We do not see changes in methylation in any of the results.
3 Expression Analysis
3.1 Running the analysis
The expression analysis can be performed in the same way than methylation. The main difference is that an ExpressionSet will be used instead of a GenomicRatioSet. We will run the same model than for methylation. We will check the effect of sex on gene expression and we will further evaluate our target region:
targetRange <- GRanges("chr11:102600000-103300000")
gexpRes <- runPipeline(set = brge_gexp, variable_names = "sex", range = targetRange)
names(gexpRes)## [1] "DiffMean" "DiffVar" "RDA"runPipeline only runs DMR methods when the input object is a GenomicRatioSet. Therefore, our object has three results: differences of means (DiffMean), differences of variances (DiffVar) and redundancy analysis (RDA).
As we did with methylation data, the first step will be check the goodness of our model. To do so, we will run a QQ-plot using the p-values of the differences of means:
plot(gexpRes, rid = "DiffMean", type = "qq")
Most of the p-values lie on the theoretical line and the lambda is close to 1 (0.94). Therefore, we will not need to include other covariates in our model.
We will now get a general overview of the results. In gene expression data, it is very common to use a Volcano plot which shows the change in expression of a probe against their p-value. We will make a Volcano plot for our two methods (DiffMean and DiffVar). We can easily make these plots using the function plot and setting the argument type to volcano:
plot(gexpRes, rid = "DiffMean", type = "volcano") + ggtitle("Differences in Means")
plot(gexpRes, rid = "DiffVar", type = "volcano") + ggtitle("Differences in Variances")
Several transcripts have a big difference in mean gene expression and a very small p-value. Transcripts with a positive fold change have a higher expression in boys than in girls and are placed in chromosome Y (name starts by TC0Y). Transcripts with a negative fold change, have a lower expression in boys than in girls and are placed in chromosome X (name starts by TC0X). On the other hand, only few transcripts have big differences in variance and the p-values are small.
We can get an overview of the genome distribution of these transcripts by plotting a Manhattan plot. We will also highlight the transcripts of our target region:
targetRange <- GRanges("chr11:102600000-103300000")
plot(gexpRes, rid = "DiffMean", type = "manhattan", main = "Differences in Means", highlight = targetRangeNum)
plot(gexpRes, rid = "DiffVar", type = "manhattan", main = "Differences in Variances", highlight = targetRangeNum)
There are several transcripts with highly significant differences in means in chromosomes X and Y. However, transcripts in our target region are not significantly changed. There are only few transcripts with differences in variance between boys and girls and all are placed in chromosome Y.
3.2 Exploring a target region
We will explore our target region using the RDA plot. As with methylation, we will color the samples by a given phenotype passing a data.frame to pheno argument:
brge_gexp$sex <- as.factor(brge_gexp$sex)
plotRDA(object = gexpRes, pheno = pData(brge_gexp)[, "sex", drop = FALSE])
We do not see a different distribution of female and male samples in the first RDA component. Our model explained very few variability (R2 = 0.009) and the p-value is big, suggesting that the variance explained by our model is not bigger than what we can expect by chance.
We can get more estimates of the RDA with the function getRDAresults:
getRDAresults(gexpRes)## R2 pval global.R2 global.pval
## 0.008512207 0.541845815 0.022901595 0.434756524The R2 of our target region is smaller than the mean R2, so our region is not particularly changed.
Finally, we will also make the region plot in our target region. In this case, we will visualize DiffMean and DiffVar results:
plotRegion(rset = gexpRes, range = targetRange)
There are no transcripts differentially expressed in the region.
4 Methylation and gene expression integration
We have included in MEAL two ways to assess association between gene expression and DNA methylation. The first option is to simultaneously examine the gene expression and DNA methylation results of a target region. We can do this with the function plotRDA and passing two ResultSets.
We will examine the methylation and gene expression results in our target region. For the sake of clarity, we will only print DiffMean and bumphunter results.
plotRegion(rset = methRes, rset2 = gexpRes, range = targetRange, results = c("bumphunter", "DiffMean"))
Differences in methylation are represented in black and differences in gene expression are represented in blue. We can see that CpGs with high differences are near transcripts with highest differences expression.
We can get a quantitative estimate of the association between methylation and gene expression using the function correlationMethExprs. This function applies a method previously described in (Steenaard et al. 2015). CpGs and expression probes are paired by proximity. Expression probes are paired to a CpG if they are completely inside a range of \(\pm\) 250Kb from the CpG. This distance of 250Kb is set by default but it can be changed with the parameter flank (whose units are bases). The correlation between methylation and expression is done by a linear regression.
To account for technical (e.g. batch) or biological (e.g. sex, age) artifacts, a model including those variables (z) is fitted: \[ x_{ij} = \sum_{k = 1}^K{\beta_{ik} z_{kj}} + r_{ij}, i = 1, ..., P \] where \(x_{ij}\) is the methylation or expression level of probe i (P is the total number of probes) for individual j, \(\beta_{ik}\) is the effect of variable k at probe i, \(z_{kj}\) is the value of variable k (K is the total number of covariates) for indivudal j and \(r_{ij}\) is the residual value of probe i and individual j. The residuals of methylation and expression are used to assess the correlation.
Therefore, let’s create a new MultiDataSet with expression and methylation data. The function createMultiDataSet creates an empty MultiDataSet. Then, we can add gene expression and methylation datasets with the function add_genexp and add_methy:
multi <- createMultiDataSet()
multi <- add_genexp(multi, brge_gexp)
multi <- add_methy(multi, brge_methy)In this analysis, we will only take into account the probes in our target region. We can subset our MultiDataSet passing a GenomicRanges to the third position of [:
multi.filt <- multi[, , targetRange]Now, we will compute the correlation between the CpGs and the expression probes.
methExprs <- correlationMethExprs(multi.filt)
head(methExprs)## cpg exprs Beta se P.Value adj.P.Val
## 1 cg08804111 TC11002235.hg.1 0.7496396 0.2143482 0.002572113 0.7462925
## 2 cg09615874 TC11002239.hg.1 -2.2805946 0.6419809 0.002275597 0.7462925
## 3 cg02344527 TC11003311.hg.1 -3.8018045 1.1008235 0.002834022 0.7462925
## 4 cg01766252 TC11002232.hg.1 7.9602956 2.4425111 0.004357524 0.8538027
## 5 cg16466334 TC11002239.hg.1 5.8514095 1.8510037 0.005403815 0.8538027
## 6 cg13759446 TC11000962.hg.1 -7.6725222 2.8959170 0.016308600 0.8747493The results shows the CpG name, the expression probe, the change of the relationship and se, and the p-value and adjusted p-value (using B&H). Results are ordered by the adjusted p-value, so we can see than in our data there are no correlated CpGs-expression probes.
5 Conclusions
This case example shows the main functionalities of MEAL package. It has been shown how to evaluate the effect of a phenotype on the methylation and on the expression at genome wide and at region level. Finally, the integration between methylation and expression is tested by checking if there are expression probes correlated to the methylated probes.
sessionInfo()## R version 3.5.0 Patched (2018-05-03 r74699)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=C
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] missMethyl_1.14.0
## [2] brgedata_1.2.0
## [3] doRNG_1.6.6
## [4] rngtools_1.3.1
## [5] pkgmaker_0.22
## [6] registry_0.5
## [7] ggplot2_2.2.1
## [8] minfiData_0.26.0
## [9] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
## [10] IlluminaHumanMethylation450kmanifest_0.4.0
## [11] minfi_1.26.0
## [12] bumphunter_1.22.0
## [13] locfit_1.5-9.1
## [14] iterators_1.0.9
## [15] foreach_1.4.4
## [16] Biostrings_2.48.0
## [17] XVector_0.20.0
## [18] SummarizedExperiment_1.10.1
## [19] DelayedArray_0.6.0
## [20] BiocParallel_1.14.1
## [21] matrixStats_0.53.1
## [22] GenomicRanges_1.32.3
## [23] GenomeInfoDb_1.16.0
## [24] IRanges_2.14.10
## [25] S4Vectors_0.18.2
## [26] MEAL_1.10.1
## [27] MultiDataSet_1.8.0
## [28] Biobase_2.40.0
## [29] BiocGenerics_0.26.0
## [30] BiocStyle_2.8.1
##
## loaded via a namespace (and not attached):
## [1] rms_5.1-2
## [2] R.utils_2.6.0
## [3] tidyselect_0.2.4
## [4] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
## [5] RSQLite_2.1.1
## [6] AnnotationDbi_1.42.1
## [7] htmlwidgets_1.2
## [8] grid_3.5.0
## [9] munsell_0.4.3
## [10] codetools_0.2-15
## [11] preprocessCore_1.42.0
## [12] statmod_1.4.30
## [13] colorspace_1.3-2
## [14] knitr_1.20
## [15] rstudioapi_0.7
## [16] labeling_0.3
## [17] GenomeInfoDbData_1.1.0
## [18] bit64_0.9-7
## [19] rhdf5_2.24.0
## [20] rprojroot_1.3-2
## [21] TH.data_1.0-8
## [22] xfun_0.1
## [23] qqman_0.1.4
## [24] biovizBase_1.28.0
## [25] R6_2.2.2
## [26] illuminaio_0.22.0
## [27] AnnotationFilter_1.4.0
## [28] bitops_1.0-6
## [29] reshape_0.8.7
## [30] assertthat_0.2.0
## [31] scales_0.5.0
## [32] bsseq_1.16.0
## [33] multcomp_1.4-8
## [34] nnet_7.3-12
## [35] gtable_0.2.0
## [36] methylumi_2.26.0
## [37] sandwich_2.4-0
## [38] ensembldb_2.4.1
## [39] rlang_0.2.0
## [40] MatrixModels_0.4-1
## [41] genefilter_1.62.0
## [42] calibrate_1.7.2
## [43] splines_3.5.0
## [44] rtracklayer_1.40.2
## [45] lazyeval_0.2.1
## [46] acepack_1.4.1
## [47] DSS_2.28.0
## [48] GEOquery_2.48.0
## [49] dichromat_2.0-0
## [50] checkmate_1.8.5
## [51] yaml_2.1.19
## [52] snpStats_1.30.0
## [53] GenomicFeatures_1.32.0
## [54] backports_1.1.2
## [55] Hmisc_4.1-1
## [56] tools_3.5.0
## [57] bookdown_0.7
## [58] nor1mix_1.2-3
## [59] RColorBrewer_1.1-2
## [60] siggenes_1.54.0
## [61] SNPassoc_1.9-2
## [62] Rcpp_0.12.17
## [63] plyr_1.8.4
## [64] base64enc_0.1-3
## [65] progress_1.1.2
## [66] zlibbioc_1.26.0
## [67] purrr_0.2.4
## [68] RCurl_1.95-4.10
## [69] BiasedUrn_1.07
## [70] prettyunits_1.0.2
## [71] rpart_4.1-13
## [72] openssl_1.0.1
## [73] zoo_1.8-1
## [74] ggrepel_0.8.0
## [75] cluster_2.0.7-1
## [76] magrittr_1.5
## [77] data.table_1.11.2
## [78] SparseM_1.77
## [79] mvtnorm_1.0-7
## [80] ProtGenerics_1.12.0
## [81] hms_0.4.2
## [82] evaluate_0.10.1
## [83] xtable_1.8-2
## [84] XML_3.98-1.11
## [85] mclust_5.4
## [86] gridExtra_2.3
## [87] compiler_3.5.0
## [88] biomaRt_2.36.1
## [89] tibble_1.4.2
## [90] R.oo_1.22.0
## [91] htmltools_0.3.6
## [92] mgcv_1.8-23
## [93] Formula_1.2-3
## [94] tidyr_0.8.1
## [95] DBI_1.0.0
## [96] MASS_7.3-50
## [97] DMRcate_1.16.0
## [98] Matrix_1.2-14
## [99] readr_1.1.1
## [100] permute_0.9-4
## [101] quadprog_1.5-5
## [102] R.methodsS3_1.7.1
## [103] Gviz_1.24.0
## [104] bindr_0.1.1
## [105] pkgconfig_2.0.1
## [106] GenomicAlignments_1.16.0
## [107] haplo.stats_1.7.9
## [108] foreign_0.8-70
## [109] xml2_1.2.0
## [110] annotate_1.58.0
## [111] multtest_2.36.0
## [112] beanplot_1.2
## [113] ruv_0.9.7
## [114] stringr_1.3.1
## [115] VariantAnnotation_1.26.0
## [116] digest_0.6.15
## [117] vegan_2.5-2
## [118] rmarkdown_1.9
## [119] base64_2.0
## [120] htmlTable_1.11.2
## [121] DelayedMatrixStats_1.2.0
## [122] curl_3.2
## [123] Rsamtools_1.32.0
## [124] gtools_3.5.0
## [125] quantreg_5.35
## [126] nlme_3.1-137
## [127] bindrcpp_0.2.2
## [128] Rhdf5lib_1.2.1
## [129] limma_3.36.1
## [130] BSgenome_1.48.0
## [131] pillar_1.2.2
## [132] lattice_0.20-35
## [133] httr_1.3.1
## [134] survival_2.42-3
## [135] GO.db_3.6.0
## [136] glue_1.2.0
## [137] DMRcatedata_1.16.0
## [138] bit_1.1-13
## [139] stringi_1.2.2
## [140] HDF5Array_1.8.0
## [141] blob_1.1.1
## [142] polspline_1.1.12
## [143] org.Hs.eg.db_3.6.0
## [144] latticeExtra_0.6-28
## [145] memoise_1.1.0
## [146] dplyr_0.7.5References
Steenaard, Rebecca V, Symen Ligthart, Lisette Stolk, Marjolein J Peters, Joyce B van Meurs, Andre G Uitterlinden, Albert Hofman, Oscar H Franco, and Abbas Dehghan. 2015. “Tobacco smoking is associated with methylation of genes related to coronary artery disease.” Clinical Epigenetics 7 (1):54. https://doi.org/10.1186/s13148-015-0088-y.