Gene set testing for Illumina HumanMethylation Arrays
Comparing the false discovery rate of different methods
Jovana Maksimovic, Alicia Oshlack and Belinda Phipson
July 21, 2020
Last updated: 2020-07-21
Checks: 7 0
Knit directory: methyl-geneset-testing/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it's best to always run the code in an empty environment.
The command set.seed(20200302) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 7fd8e1c. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/figures.nb.html
Ignored: code/.job/
Ignored: code/old/
Ignored: data/
Ignored: output/.DS_Store
Ignored: output/450K.rds
Ignored: output/CD4vCD8.GO.csv
Ignored: output/CD4vCD8.KEGG.csv
Ignored: output/EPIC.rds
Ignored: output/FDR-analysis/
Ignored: output/compare-methods/
Ignored: output/figures/
Ignored: output/random-cpg-sims/
Unstaged changes:
Modified: .gitignore
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/fdrAnalysis.Rmd) and HTML (docs/fdrAnalysis.html) files. If you've configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 976b2b5 | JovMaksimovic | 2020-07-17 | Build site. |
| Rmd | cfcac1b | Jovana Maksimovic | 2020-07-09 | Updated paper figures. |
| html | 8454e85 | Jovana Maksimovic | 2020-06-23 | Build site. |
| Rmd | 6aa0ab0 | Jovana Maksimovic | 2020-06-23 | wflow_publish(c("analysis/fdrAnalysis.Rmd")) |
| Rmd | 823b793 | Jovana Maksimovic | 2020-06-16 | UPdated analyses with label translation and new colour palette. |
| html | 22f00e9 | Jovana Maksimovic | 2020-05-19 | Build site. |
| Rmd | 63b0011 | Jovana Maksimovic | 2020-05-19 | wflow_publish(c("analysis/exploreArrayBias450.Rmd", "analysis/exploreArrayBiasEPIC.Rmd", |
| html | 7af677c | Jovana Maksimovic | 2020-04-17 | Build site. |
| Rmd | 0fec41a | Jovana Maksimovic | 2020-04-17 | Completed repeat of FDR analysis using TCGA data |
| Rmd | f42dc82 | JovMaksimovic | 2020-04-14 | Updated all FDR analysis code. |
| Rmd | 1ea3283 | JovMaksimovic | 2020-04-14 | Updated null analysis with KIRC TCGA data |
library(here)
library(minfi)
library(paletteer)
library(limma)
library(reshape2)
library(missMethyl)
library(ggplot2)
library(glue)
library(tibble)
library(dplyr)
library(curatedTCGAData)
library(MultiAssayExperiment)
library(TCGAutils)
library(DMRcate)
library(patchwork)
source(here("code/utility.R"))Load data
We are using publicly available kidney clear cell carcinoma (KIRC) 450k data from The Cancer Genome Atlas (TCGA). We are using only the normal samples to look at false discovery rate control by various methylation gene set testing methods.
Download the data using the curatedTCGAData Bioconductor package and then extract the normal samples. The data seems to have already been normalised(?) so we will only remove poor probes (>1 NA), SNP, multi-mapping and sex-chromosome probes.
kirc <- curatedTCGAData(diseaseCode = "KIRC", assays = "Methylation_methyl450",
dry.run = FALSE)
kirc <- splitAssays(kirc, c("11")) # extract only the normal samples
exp <- experiments(kirc)[[1]]
meta <- colData(kirc)
betas <- as.matrix(assay(exp))
colnames(betas) <- substr(colnames(betas), 1, 12)
m <- match(colnames(betas), meta$patientID)
meta <- meta[m, ]
head(meta[, 1:5])betasNoNA <- betas[rowSums(is.na(betas)) == 0, ]
mds <- plotMDS(betasNoNA, gene.selection = "common", plot = FALSE)
dat <- tibble(x = mds$x, y = mds$y, gender = meta$gender)
ggplot(dat, aes(x = x, y = y, colour = gender)) +
geom_point() +
labs(x = "Principal Component 1", y = "Principal Component 2",
colour = "Sex")
betasFlt <- rmSNPandCH(betasNoNA, rmXY = TRUE, rmcrosshyb = TRUE)snapshotDate(): 2019-10-22see ?DMRcatedata and browseVignettes('DMRcatedata') for documentationloading from cachesee ?DMRcatedata and browseVignettes('DMRcatedata') for documentationloading from cachesee ?DMRcatedata and browseVignettes('DMRcatedata') for documentationloading from cachedim(betasFlt)[1] 364602 160mds <- plotMDS(betasFlt, gene.selection = "common", plot = FALSE)
dat <- tibble(x = mds$x, y = mds$y, gender = meta$gender)
pal <- paletteer::paletteer_d("wesanderson::Moonrise3", 2)
cols <- c("female" = pal[2], "male" = pal[1])
p <- ggplot(dat, aes(x = x, y = y, colour = gender)) +
geom_point(size = 2) +
labs(x = "Principal Component 1", y = "Principal Component 2",
colour = "Sex") +
scale_color_manual(values = cols)
p
fig <- here("output/figures/Fig-3C.rds")
saveRDS(p, fig, compress = FALSE)dat$age <- meta$years_to_birth
dat$race <- sub(" or", "\nor", meta$race)
dat$ethnicity <- sub(" or", "\nor", meta$ethnicity)
a <- ggplot(dat, aes(x = x, y = y, colour = age)) +
geom_point() +
labs(x = "Principal Component 1", y = "Principal Component 2",
colour = "Age") +
viridis::scale_color_viridis(direction = -1)
b <- ggplot(dat, aes(x = x, y = y, colour = race)) +
geom_point() +
labs(x = "Principal Component 1", y = "Principal Component 2",
colour = "Race") +
theme(legend.text = element_text(size = 7))
c <- ggplot(dat, aes(x = x, y = y, colour = ethnicity)) +
geom_point() +
labs(x = "Principal Component 1", y = "Principal Component 2",
colour = "Ethnicity") +
theme(legend.text = element_text(size = 7))
(b + c) / a
dat <- as_tibble(melt(betasFlt))
colnames(dat) <- c("cpg", "ID", "beta")
p <- ggplot(dat, aes(x = beta)) +
geom_line(aes(color = ID), stat="density", size=0.5, alpha=0.5,
show.legend = FALSE)
p + labs(x = "Beta value", y = "Density")
FDR analysis
Run 100 null simulations comparing 5, 10, 20, 40, 80 randomly selected samples per "group". There should be no significant differential methylation between the groups as the data contains no signal. Consequently, we expect about 5% false gene set discoveries across the 100 simulations from methods that are "holding their size". We will compare gometh (with top 1000 and 5000 ranked CpGs selected), MethylGSA:GLM, MethylGSM:ORA, MethylGSA:GSEA and ChAMP:ebGSEA for the complete list of BROAD gene sets provided in the ChAMP package.
outFile <- here("data/TCGA.KIRC.rds")
if(!file.exists(outFile)){
saveRDS(betasFlt, file = outFile)
}Load the results of the FDR simulations.
inFile <- here("output/FDR-analysis/FDR-analysis.rds")
if(file.exists(inFile)) dat <- as_tibble(readRDS(inFile))The plots below shows that MethylGSA-ORA does not control the false discovery rate very well as the median proportion of p-value 0.05 for the 100 simulations is greater than 0.5. MethylGSA-GSEA does better, although its median false discovery rate is still relatively high at around 0.25. MethylGSA-GLM has good false discovery rate control with median FDR at around 0.05. ChAMP-ebGSEA is only slightly anti-conservative using both the wicox test (WT) and known population median test (KPMT) with a median FDR at around 0.06-0.08. MissMethyl-gometh is very consistent although somewhat conservative with a median FDR at around 0.02-0.03.
dat %>% mutate(method = unname(dict[method])) %>%
group_by(simNo, sampleNo, method) %>%
summarise(psig = sum(pvalue < 0.05)/length(pvalue)) %>%
mutate(sampleOrd = as.integer(sampleNo)) -> sigDat`summarise()` regrouping output by 'simNo', 'sampleNo' (override with `.groups` argument)p <- ggplot(sigDat, aes(x = reorder(sampleNo, sampleOrd), y = psig,
fill = method)) +
geom_violin(scale = "width", width = 0.8, size = 0.3) +
stat_summary(geom = "point", size = 0.5, color = "white",
position = position_dodge(0.8),
show.legend = FALSE, fun = median) +
geom_hline(yintercept=0.05, linetype="dashed", color = "red") +
labs(y="Prop. gene sets with p-value < 0.05", x="No. samples per group",
fill="Method") +
scale_fill_manual(values = methodCols)
p
| Version | Author | Date |
|---|---|---|
| 976b2b5 | JovMaksimovic | 2020-07-17 |
fig <- here("output/figures/Fig-3D.rds")
saveRDS(p, fig, compress = FALSE)
sessionInfo()R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] DMRcatedata_2.2.1 ExperimentHub_1.12.0
[3] AnnotationHub_2.18.0 BiocFileCache_1.10.2
[5] dbplyr_1.4.4 rhdf5_2.30.1
[7] patchwork_1.0.1 DMRcate_2.0.7
[9] TCGAutils_1.6.2 curatedTCGAData_1.8.1
[11] MultiAssayExperiment_1.12.6 dplyr_1.0.0
[13] tibble_3.0.3 glue_1.4.1
[15] ggplot2_3.3.2 missMethyl_1.20.4
[17] reshape2_1.4.4 limma_3.42.2
[19] paletteer_1.2.0 minfi_1.32.0
[21] bumphunter_1.28.0 locfit_1.5-9.4
[23] iterators_1.0.12 foreach_1.5.0
[25] Biostrings_2.54.0 XVector_0.26.0
[27] SummarizedExperiment_1.16.1 DelayedArray_0.12.3
[29] BiocParallel_1.20.1 matrixStats_0.56.0
[31] Biobase_2.46.0 GenomicRanges_1.38.0
[33] GenomeInfoDb_1.22.1 IRanges_2.20.2
[35] S4Vectors_0.24.4 BiocGenerics_0.32.0
[37] here_0.1 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] R.utils_2.9.2
[2] tidyselect_1.1.0
[3] htmlwidgets_1.5.1
[4] RSQLite_2.2.0
[5] AnnotationDbi_1.48.0
[6] grid_3.6.3
[7] munsell_0.5.0
[8] codetools_0.2-16
[9] preprocessCore_1.48.0
[10] statmod_1.4.34
[11] withr_2.2.0
[12] colorspace_1.4-1
[13] knitr_1.29
[14] rstudioapi_0.11
[15] labeling_0.3
[16] git2r_0.27.1
[17] GenomeInfoDbData_1.2.2
[18] farver_2.0.3
[19] bit64_0.9-7.1
[20] rprojroot_1.3-2
[21] vctrs_0.3.2
[22] generics_0.0.2
[23] xfun_0.15
[24] biovizBase_1.34.1
[25] R6_2.4.1
[26] illuminaio_0.28.0
[27] AnnotationFilter_1.10.0
[28] bitops_1.0-6
[29] reshape_0.8.8
[30] assertthat_0.2.1
[31] bsseq_1.22.0
[32] promises_1.1.1
[33] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
[34] scales_1.1.1
[35] nnet_7.3-14
[36] gtable_0.3.0
[37] methylumi_2.32.0
[38] ensembldb_2.10.2
[39] rlang_0.4.7
[40] genefilter_1.68.0
[41] splines_3.6.3
[42] lazyeval_0.2.2
[43] rtracklayer_1.46.0
[44] DSS_2.34.0
[45] acepack_1.4.1
[46] GEOquery_2.54.1
[47] dichromat_2.0-0
[48] prismatic_0.2.0
[49] checkmate_2.0.0
[50] BiocManager_1.30.10
[51] yaml_2.2.1
[52] GenomicFeatures_1.38.2
[53] backports_1.1.8
[54] httpuv_1.5.4
[55] Hmisc_4.4-0
[56] tools_3.6.3
[57] nor1mix_1.3-0
[58] ellipsis_0.3.1
[59] RColorBrewer_1.1-2
[60] siggenes_1.60.0
[61] Rcpp_1.0.5
[62] plyr_1.8.6
[63] base64enc_0.1-3
[64] progress_1.2.2
[65] zlibbioc_1.32.0
[66] purrr_0.3.4
[67] RCurl_1.98-1.2
[68] BiasedUrn_1.07
[69] prettyunits_1.1.1
[70] rpart_4.1-15
[71] openssl_1.4.2
[72] IlluminaHumanMethylationEPICmanifest_0.3.0
[73] viridis_0.5.1
[74] cluster_2.1.0
[75] fs_1.4.2
[76] magrittr_1.5
[77] data.table_1.12.8
[78] whisker_0.4
[79] ProtGenerics_1.18.0
[80] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
[81] hms_0.5.3
[82] mime_0.9
[83] evaluate_0.14
[84] xtable_1.8-4
[85] XML_3.99-0.3
[86] jpeg_0.1-8.1
[87] readxl_1.3.1
[88] mclust_5.4.6
[89] gridExtra_2.3
[90] compiler_3.6.3
[91] biomaRt_2.42.1
[92] crayon_1.3.4
[93] R.oo_1.23.0
[94] htmltools_0.5.0
[95] later_1.1.0.1
[96] Formula_1.2-3
[97] tidyr_1.1.0
[98] DBI_1.1.0
[99] MASS_7.3-51.6
[100] rappdirs_0.3.1
[101] Matrix_1.2-18
[102] readr_1.3.1
[103] permute_0.9-5
[104] R.methodsS3_1.8.0
[105] quadprog_1.5-8
[106] Gviz_1.30.3
[107] pkgconfig_2.0.3
[108] GenomicAlignments_1.22.1
[109] foreign_0.8-76
[110] IlluminaHumanMethylation450kmanifest_0.4.0
[111] xml2_1.3.2
[112] annotate_1.64.0
[113] rngtools_1.5
[114] multtest_2.42.0
[115] beanplot_1.2
[116] ruv_0.9.7.1
[117] rvest_0.3.5
[118] doRNG_1.8.2
[119] scrime_1.3.5
[120] VariantAnnotation_1.32.0
[121] stringr_1.4.0
[122] digest_0.6.25
[123] cellranger_1.1.0
[124] rmarkdown_2.3
[125] base64_2.0
[126] htmlTable_2.0.1
[127] edgeR_3.28.1
[128] DelayedMatrixStats_1.8.0
[129] curl_4.3
[130] gtools_3.8.2
[131] shiny_1.5.0
[132] Rsamtools_2.2.3
[133] lifecycle_0.2.0
[134] nlme_3.1-148
[135] GenomicDataCommons_1.10.0
[136] jsonlite_1.7.0
[137] Rhdf5lib_1.8.0
[138] viridisLite_0.3.0
[139] BSgenome_1.54.0
[140] askpass_1.1
[141] pillar_1.4.6
[142] lattice_0.20-41
[143] fastmap_1.0.1
[144] httr_1.4.2
[145] survival_3.2-3
[146] GO.db_3.10.0
[147] interactiveDisplayBase_1.24.0
[148] png_0.1-7
[149] BiocVersion_3.10.1
[150] bit_1.1-15.2
[151] stringi_1.4.6
[152] HDF5Array_1.14.4
[153] rematch2_2.1.2
[154] blob_1.2.1
[155] org.Hs.eg.db_3.10.0
[156] latticeExtra_0.6-29
[157] memoise_1.1.0