This script filters out unreliable and outlying probes and data points.
Load packages
library(tidyverse)
library(DNAmArray)Load data
load("../GOTO_Data/Processing/GOTO_targets-unfiltered.Rdata")
load("../GOTO_Data/Processing/GOTO_RGset-unfiltered.Rdata")
load("../GOTO_Data/Processing/GOTO_GRset-unfiltered.Rdata")
load("../GOTO_Data/Processing/GOTO_methData-unfiltered.Rdata")
methData_unfiltered## class: SummarizedExperiment
## dim: 865859 534
## metadata(0):
## assays(1): beta
## rownames(865859): cg14817997 cg26928153 ... cg07587934 cg16855331
## rowData names(4): cpg chr start end
## colnames(534): 203527980082_R01C01 203527980082_R02C01 ...
## 203550300093_R07C01 203550300093_R08C01
## colData names(19): DNA_labnr IOP2_ID ... Basename smokeRGset## class: RGChannelSetExtended
## dim: 1051815 534
## metadata(0):
## assays(5): Green Red GreenSD RedSD NBeads
## rownames(1051815): 1600101 1600111 ... 99810990 99810992
## rowData names(0):
## colnames(534): 203527980082_R01C01 203527980082_R02C01 ...
## 203550300093_R07C01 203550300093_R08C01
## colData names(19): DNA_labnr IOP2_ID ... Basename filenames
## Annotation
## array: IlluminaHumanMethylationEPIC
## annotation: ilm10b4.hg19Load annotation
anno <- read_tsv(
"/exports/molepi/users/ljsinke/LLS/Shared_Data/Manifests/EPIC.hg19.manifest.tsv.gz")Remove probes with zero intensity or based on less than 3 beads
RGset <- probeFiltering(RGset)## Filtering on number of beads...
## On average 0.2 % of the probes ( 1051815 ) have number of beads below 3
## Filtering on zero intensities...
## On average 0.003 % of the Type II probes ( 723910 ) have Red and/or Green intensity below 1
## On average 0.037 % of the Type I probes ( 99908 ), measured in Green channel, have intensity below 1
## On average 0.042 % of the Type I probes ( 184454 ), measured in Red channel, have intensity below 1
## Set filtered probes in Red and/or Green channels to NA...
## ... done 100 out of 534 ...
## ... done 200 out of 534 ...
## ... done 300 out of 534 ...
## ... done 400 out of 534 ...
## ... done 500 out of 534 ...Remove probes where the detection p-value is less than 0.01, so they are not significantly different from background noise.
This function additionally removes samples and probes with a success rate below 5%
betas <- reduce(GRset, RGset, what="beta")## Calculate and filter on detection P-value...
## On average 0.37 % of the CpGs ( 866091 ) have detection P-value above the threshold 0.01
## Transform to beta -values...
## On average 0.23 % of the probes ( 866091 ) were set to NA in the probe filtering step!
## Calculate success rates and reduce...
## Percentage of samples having success rate above 0.95 is 100 %
## Percentage of CpGs having success rate above 0.95 is 98.94 %anno <- anno %>%
dplyr::select(
cpg = probeID,
chr = CpG_chrm,
start = CpG_beg,
end = CpG_end,
everything()
) %>%
dplyr::filter(
cpg %in% rownames(betas)
)
dim(anno)## [1] 856881 57Save CpGs on sex chromosomes
xy_cpgs <- (anno %>%
dplyr::filter(chr %in% c("X", "Y")))$cpgProbes from the X chromosome and the Y chromosome are also excluded
dim(betas)## [1] 856881 534betas <- betas[!rownames(betas) %in% xy_cpgs,]
dim(betas)## [1] 856881 534Load in the CpGs within ENCODE blacklist regions
load("../../../LLS/Shared_Data/ENCODE_Blacklist/ENCODEBlacklist_CpGomit-EPIC.RData")These 8,552 probes are excluded from the methylation data
betas <- betas[!rownames(betas) %in% cpg_blacklist,]
dim(betas)## [1] 848029 534For general use, there are 97,194 probes to be excluded
maskProbes <- anno[anno$MASK_general == TRUE,]$cpg
length(maskProbes)## [1] 97194Remove 92,252 of these from the GOTO data
betas <- betas[!rownames(betas) %in% maskProbes,]
dim(betas)## [1] 755777 534Any data points outside the range of LQ - 3 x IQR to UQ + 3 x IQR are excluded. This included 1,124,127 values that were converted to NA.
xtabs(~is.na(betas))## is.na(betas)
## FALSE TRUE
## 402728308 856610iqr_dnam <- apply(betas, 1, function(x){
iqr <- IQR(x, na.rm = TRUE)
q1 <- quantile(x, na.rm=TRUE)[[2]]
q3 <- quantile(x, na.rm=TRUE)[[4]]
x <- ifelse((x <= q1 - (3*iqr) | x >= q3 + (3*iqr)), NA, x)
})
xtabs(~is.na(iqr_dnam))## is.na(iqr_dnam)
## FALSE TRUE
## 401604181 1980737dim(iqr_dnam)## [1] 534 755777We replace the betas with the new beta values.
betas <- t(iqr_dnam)Keep only the CpGs where less than 5% of sample values are missing
This is all CpGs, representing high quality data
perc_na <- rowSums(is.na(iqr_dnam))*100/ncol(iqr_dnam)
summary(perc_na)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1692 0.2850 0.3610 0.4908 0.5072 8.2690betas <- betas[,perc_na <= 95]
dim(betas)## [1] 755777 534Keep only the samples where less than 5% of CpGs are missing
This is all samples, representing high quality data
perc_na <- colSums(is.na(iqr_dnam))*100/nrow(iqr_dnam)
summary(perc_na)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.1873 0.4908 0.5618 41.9476betas <- betas[perc_na <= 95,]
dim(betas)## [1] 755777 534Annotation data
anno <- anno %>%
dplyr::filter(cpg %in% rownames(betas))Sample information
targets <- targets %>%
dplyr::filter(
Basename %in% colnames(betas))Make summarized experiment object
methData <- SummarizedExperiment(
assays=SimpleList(beta=betas),
rowData=anno,
colData=targets)
methData## class: SummarizedExperiment
## dim: 755777 534
## metadata(0):
## assays(1): beta
## rownames(755777): cg18478105 cg09835024 ... cg10633746 cg12623625
## rowData names(57): cpg chr ... MASK_extBase MASK_general
## colnames(534): 203527980082_R01C01 203527980082_R02C01 ...
## 203550300093_R07C01 203550300093_R08C01
## colData names(19): DNA_labnr IOP2_ID ... Basename smokeSave RGset
RGset <- RGset[, colnames(RGset) %in% colnames(methData)]
RGset## class: RGChannelSet
## dim: 1051815 534
## metadata(0):
## assays(2): Green Red
## rownames(1051815): 1600101 1600111 ... 99810990 99810992
## rowData names(0):
## colnames(534): 203527980082_R01C01 203527980082_R02C01 ...
## 203550300093_R07C01 203550300093_R08C01
## colData names(19): DNA_labnr IOP2_ID ... Basename filenames
## Annotation
## array: IlluminaHumanMethylationEPIC
## annotation: ilm10b4.hg19Save clean data
save(targets, file = "../GOTO_Data/GOTO_targets-filtered.Rdata")
save(methData, file="../GOTO_Data/GOTO_methData-filtered.Rdata")
save(RGset, file="../GOTO_Data/GOTO_RGset-filtered.Rdata")sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.10 (Green Obsidian)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [2] snpStats_1.44.0
## [3] survival_3.5-5
## [4] ggrepel_0.9.1
## [5] ggfortify_0.4.14
## [6] irlba_2.3.5.1
## [7] Matrix_1.5-4.1
## [8] omicsPrint_1.14.0
## [9] MASS_7.3-60
## [10] DNAmArray_2.0.0
## [11] pls_2.8-2
## [12] FDb.InfiniumMethylation.hg19_2.2.0
## [13] org.Hs.eg.db_3.14.0
## [14] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [15] GenomicFeatures_1.46.5
## [16] AnnotationDbi_1.56.2
## [17] IlluminaHumanMethylationEPICmanifest_0.3.0
## [18] minfi_1.40.0
## [19] bumphunter_1.36.0
## [20] locfit_1.5-9.8
## [21] iterators_1.0.14
## [22] foreach_1.5.2
## [23] Biostrings_2.62.0
## [24] XVector_0.34.0
## [25] SummarizedExperiment_1.24.0
## [26] Biobase_2.58.0
## [27] MatrixGenerics_1.10.0
## [28] matrixStats_1.0.0
## [29] GenomicRanges_1.46.1
## [30] GenomeInfoDb_1.34.9
## [31] IRanges_2.32.0
## [32] S4Vectors_0.36.2
## [33] BiocGenerics_0.44.0
## [34] BiocParallel_1.32.6
## [35] MethylAid_1.28.0
## [36] forcats_0.5.2
## [37] stringr_1.5.0
## [38] dplyr_1.1.3
## [39] purrr_0.3.4
## [40] readr_2.1.2
## [41] tidyr_1.2.1
## [42] tibble_3.2.1
## [43] ggplot2_3.4.3
## [44] tidyverse_1.3.2
## [45] rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.3 tidyselect_1.2.0
## [3] RSQLite_2.2.17 grid_4.2.2
## [5] munsell_0.5.0 codetools_0.2-19
## [7] preprocessCore_1.60.2 withr_2.5.0
## [9] colorspace_2.1-0 filelock_1.0.2
## [11] highr_0.10 knitr_1.43
## [13] rstudioapi_0.14 labeling_0.4.2
## [15] GenomeInfoDbData_1.2.9 farver_2.1.1
## [17] bit64_4.0.5 rhdf5_2.42.1
## [19] vctrs_0.6.3 generics_0.1.3
## [21] xfun_0.39 timechange_0.2.0
## [23] BiocFileCache_2.2.1 R6_2.5.1
## [25] illuminaio_0.40.0 bitops_1.0-7
## [27] rhdf5filters_1.10.1 cachem_1.0.8
## [29] reshape_0.8.9 DelayedArray_0.24.0
## [31] assertthat_0.2.1 vroom_1.5.7
## [33] promises_1.2.0.1 BiocIO_1.8.0
## [35] scales_1.2.1 googlesheets4_1.0.1
## [37] gtable_0.3.3 rlang_1.1.1
## [39] genefilter_1.76.0 splines_4.2.2
## [41] rtracklayer_1.54.0 gargle_1.5.0
## [43] GEOquery_2.62.2 htm2txt_2.2.2
## [45] hexbin_1.28.3 broom_1.0.1
## [47] yaml_2.3.7 reshape2_1.4.4
## [49] RaggedExperiment_1.18.0 modelr_0.1.9
## [51] backports_1.4.1 httpuv_1.6.11
## [53] tools_4.2.2 gridBase_0.4-7
## [55] nor1mix_1.3-0 ellipsis_0.3.2
## [57] jquerylib_0.1.4 RColorBrewer_1.1-3
## [59] siggenes_1.68.0 MultiAssayExperiment_1.20.0
## [61] Rcpp_1.0.10 plyr_1.8.8
## [63] sparseMatrixStats_1.10.0 progress_1.2.2
## [65] zlibbioc_1.44.0 RCurl_1.98-1.12
## [67] prettyunits_1.1.1 openssl_2.0.6
## [69] haven_2.5.1 fs_1.6.2
## [71] magrittr_2.0.3 data.table_1.14.8
## [73] reprex_2.0.2 googledrive_2.0.0
## [75] hms_1.1.2 mime_0.12
## [77] evaluate_0.21 xtable_1.8-4
## [79] XML_3.99-0.14 mclust_6.0.0
## [81] readxl_1.4.1 gridExtra_2.3
## [83] compiler_4.2.2 biomaRt_2.50.3
## [85] crayon_1.5.2 htmltools_0.5.5
## [87] later_1.3.1 tzdb_0.4.0
## [89] lubridate_1.9.2 DBI_1.1.3
## [91] dbplyr_2.2.1 rappdirs_0.3.3
## [93] cli_3.6.1 quadprog_1.5-8
## [95] pkgconfig_2.0.3 GenomicAlignments_1.30.0
## [97] xml2_1.3.4 annotate_1.72.0
## [99] bslib_0.5.0 rngtools_1.5.2
## [101] multtest_2.50.0 beanplot_1.3.1
## [103] rvest_1.0.3 doRNG_1.8.6
## [105] scrime_1.3.5 digest_0.6.31
## [107] base64_2.0.1 cellranger_1.1.0
## [109] DelayedMatrixStats_1.16.0 restfulr_0.0.15
## [111] curl_5.0.1 shiny_1.7.2
## [113] Rsamtools_2.10.0 rjson_0.2.21
## [115] lifecycle_1.0.3 nlme_3.1-162
## [117] jsonlite_1.8.5 Rhdf5lib_1.20.0
## [119] askpass_1.1 limma_3.54.2
## [121] fansi_1.0.4 pillar_1.9.0
## [123] lattice_0.21-8 KEGGREST_1.34.0
## [125] fastmap_1.1.1 httr_1.4.6
## [127] glue_1.6.2 png_0.1-8
## [129] bit_4.0.5 stringi_1.7.12
## [131] sass_0.4.6 HDF5Array_1.22.1
## [133] blob_1.2.4 memoise_2.0.1Clear
rm(list=ls())