Estill 2016 reanalysis
The methylation in ART meta-analysis group
2024-01-03
Introduction
In this report, I will take you through a re-analysis methylation data first described by Estill et al (2016).
In their study, they analysed the DNA methylation patterns of 137 neonatal blood spots conceived naturally (NAT), through through insemination (IUI), or through ICSI using fresh or cryopreserved (frozen) embryo transfer (FH and FZ respectively).
The platform used in the study is the Illumina Infinium HumanMethylation450k BeadChip assay. The authors used a pipeline based on ChAMP (Tian et al, 2017), together with Adjacent Site Clustering (Sofer et al, 2013).
The methylation data have been deposited to NCBI GEO repository accession number GSE79257.
The main conclusions from the original study were:
The methylation profiles of assisted reproductive technology and IUI newborns were dramatically different from those of naturally (in vivo) conceived newborns.
Profiles of ICSI-frozen (FET) and IUI infants were strikingly similar, suggesting that cryopreservation may temper some of the epigenetic aberrations induced by IVF or ICSI.
The DNA methylation changes associated with IVF/ICSI culture conditions and/or parental infertility were detected at metastable epialleles, suggesting a lasting impact on a child’s epigenome.
Both infertility and ICSI alter DNA methylation at specific genomic loci, an effect that is mitigated to some extent by FET (freezing).
The aim of this work is to;
develop the analytical pipelines required for efficient re-analysis of 450K array data,
to confirm that we are able to obtain differential methylation results that are similar to those obtained in the original study, and
to critically evaluate the conclusions made in the original study.
In this report I will be using the missMethyl vignette as a guide to analyse this dataset (Phipson et al, 2015).
Previously I curated a list of good and bad probes using the script “filter_probes.Rmd”.
Loading packages and functions
These packackes will help us to perform vital steps such as normalisation, filtering, differential analysis, etc, and provide information about the array probe annotaions.
These functions provide shortcuts to help with charts and other analysis. They will eventually be shoved into another Rscript or package but can stay here for now.
source("meth_functions.R")
# Annotation
annofull <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
anno <- data.frame(annofull[,c("UCSC_RefGene_Name","Regulatory_Feature_Group",
"Islands_Name","Relation_to_Island")])Data import
dir.create("GSE79257")## Warning in dir.create("GSE79257"): 'GSE79257' already existsARRAY_SAMPLESHEET="GSE79257/GSE79257_Illumina_samplesheet.csv"
# only download it if it is not present on the system
if ( !file.exists(ARRAY_SAMPLESHEET ) ) {
DLFILE=paste(ARRAY_SAMPLESHEET,".gz",sep="")
options(timeout=1000000)
download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE79nnn/GSE79257/suppl/GSE79257_Illumina_samplesheet.csv.gz",
destfile = DLFILE)
gunzip(DLFILE, overwrite = TRUE )
}
ARRAY_DATA="GSE79257/GSE79257_RAW.tar"
# only download it if it is not present on the system
if ( !dir.exists("GSE79257/IDAT") ) {
dir.create("GSE79257/IDAT")
if ( !file.exists(ARRAY_DATA) ) {
system('wget -O idats.tar.gz "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE79257&format=file"')
file.rename("idats.tar.gz",ARRAY_DATA)
untar(exdir = "GSE79257/IDAT", tarfile = ARRAY_DATA)
}
}
baseDir <- "GSE79257"
targets <- read.metharray.sheet(baseDir)## [read.metharray.sheet] Found the following CSV files:## [1] "GSE79257/GSE79257_Illumina_samplesheet.csv"head(targets)## Sample_Name Sample_Plate Sample_Group Pool_ID Project_ID Sample_Well Array
## 1 20001_D 2 NAT_Male <NA> NA R04C01 R04C01
## 2 20028_D 2 NAT_Male <NA> NA R01C01 R01C01
## 3 20029_D 2 NAT_Male <NA> NA R02C02 R02C02
## 4 20030_D 2 NAT_Male <NA> NA R04C02 R04C02
## 5 20031_D 2 NAT_Male <NA> NA R01C02 R01C02
## 6 20032_D 1 NAT_Male <NA> NA R03C02 R03C02
## Slide Basename
## 1 9266441185 GSE79257/IDAT/GSM2089923_9266441185_R04C01
## 2 9266441185 GSE79257/IDAT/GSM2089924_9266441185_R01C01
## 3 9266441185 GSE79257/IDAT/GSM2089925_9266441185_R02C02
## 4 9266441186 GSE79257/IDAT/GSM2089926_9266441186_R04C02
## 5 9285451032 GSE79257/IDAT/GSM2089927_9285451032_R01C02
## 6 9297949059 GSE79257/IDAT/GSM2089928_9297949059_R03C02rgSet <- read.metharray.exp(targets = targets)Normalisation
mSet <- preprocessRaw(rgSet)
mSetSw <- SWAN(mSet,verbose=FALSE)## [SWAN] Preparing normalization subset
## 450k
## [SWAN] Normalizing methylated channel
## [SWAN] Normalizing unmethylated channelpar(mfrow=c(1,2), cex=0.8)
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")
Figure 1. Normalisation of bead-array data with SWAN.
Filter probes
Here we are running parallel analyses, both including and excluding sex chromosomes.
# include sex chromosomes
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetSw <- mSetSw[keep,]
# exclude SNP probes
mSetSw <- mapToGenome(mSetSw)
mSetSw_nosnp <- dropLociWithSnps(mSetSw)
dim(mSetSw)## [1] 443589 137dim(mSetSw_nosnp)## [1] 428936 137mSetSw <- mSetSw_nosnp
# exclude sex chromosomes
keep <- !(featureNames(mSetSw) %in% annofull$Name[annofull$chr %in% c("chrX","chrY")])
mSetFlt <- mSetSw[keep,]
head(mSetFlt)## class: GenomicMethylSet
## dim: 6 137
## metadata(0):
## assays(2): Meth Unmeth
## rownames(6): cg13869341 cg14008030 ... cg00381604 cg20253340
## rowData names(0):
## colnames(137): GSM2089923_9266441185_R04C01
## GSM2089924_9266441185_R01C01 ... GSM2090058_9611519080_R05C01
## GSM2090059_9611519080_R05C02
## colData names(10): Sample_Name Sample_Plate ... Basename filenames
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19
## Preprocessing
## Method: SWAN (based on a MethylSet
## preprocessed as 'Raw (no normalization or bg correction)')
## minfi version: 1.48.0
## Manifest version: 0.4.0dim(mSetFlt)## [1] 418833 137Extracting Beta and M-values
# include sex chromosomes
meth <- getMeth(mSetSw)
unmeth <- getUnmeth(mSetSw)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mSetSw)
# exclude sex chromosomes
meth <- getMeth(mSetFlt)
unmeth <- getUnmeth(mSetFlt)
Mval_flt <- log2((meth + 100)/(unmeth + 100))
beta_flt <- getBeta(mSetFlt)MDS analysis
[Multidimensional scaling(https://en.wikipedia.org/wiki/Multidimensional_scaling) plot is a method used to identify the major sources of variation in a dataset. In the MDS plots below, I will be plotting the first two dimensions (principal components [PCs]), with each sample label coloured either by ART classification, sex, ART and sex, and then array chip and then sample plate.
We wil begin with MDS analysis including the sex chromosomes and then exclude them.
First, let’s quantify the contribution of the major principal components. with a scree plot, we can see whether most of the variation is captured in the first two PCs or whether it is spread over more PCs. As we can see in Figure 2, the main source of variation is what is shown in PC1, and a much lesser extent on the other dimensions. Interestingly, excluding sex chromosomes does not seem to change the relative contributions of PCs very much.
par(mfrow=c(2,1))
myscree(Mval,main="incl sex chr")
myscree(Mval_flt,main="excl sex chr")
Figure 2. Screeplot shows contribution of top principal components to the overall variation in the experiment.
Here is the MDS plot by ART classification (Figure 3). You can see that there are four clusters of samples when MDS is projected this way. Natural birth infants (NAT) are mostly seen in the two clusters at the bottom of the chart and the IVF frozen, IVF fresh and in intrauterine insemination (IUI) appear mostly as two clusters at the top of the chart. IUI seem to have a few samples that appear as intermediate between the upper and lower clusters.
When excluding sex chromosomes, the samples tend to form two clusters with the NAT samples on the left of the chart and others on the right of the chart. IUI and FH have some samples which have intermediate methylation, which points to a slight effect of freezing.
targets$sex <- factor(sapply(strsplit(targets$Sample_Group,"_"),"[[",2))
targets$art <- factor(sapply(strsplit(targets$Sample_Group,"_"),"[[",1))
sample_groups <- factor(targets$art)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$art))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by ART type")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
Figure 3. MDS plot coloured by ART classification.
mydist <- plotMDS(Mval, labels=targets$Sample_Name,col=colours,main="sex chromosomes included")
Figure 3. MDS plot coloured by ART classification.
mydist_flt <- plotMDS(Mval_flt, labels=targets$Sample_Name,col=colours,main="sex chromosomes excluded")
Figure 3. MDS plot coloured by ART classification.
Next, we created an MDS plot by sex (Figure 4). The female samples appear in the two cluster on the right hand side of the chart and the male sample on the left. This could be simply due to the chromosomal makeup of the male and female samples or could also be due to autosomes. To check this, the MDS plot needs to be repeated with autosomes excluded.
When excluding the sex chromosomes, there is no clear clustering of the samples by sex. It appears that ART classification is the dominant source of variation.
sample_groups <- factor(targets$sex)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")## Warning in brewer.pal(n = length(levels(sample_groups)), name = "Paired"): minimal value for n is 3, returning requested palette with 3 different levelscolours <- colour_palette[as.integer(factor(targets$sex))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by sex")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
Figure 4. MDS plot coloured by sex.
plotMDS(mydist, labels=targets$Sample_Name,col=colours,main="sex chromosomes included")
Figure 4. MDS plot coloured by sex.
plotMDS(mydist_flt, labels=targets$Sample_Name,col=colours,main="sex chromosomes excluded")
Figure 4. MDS plot coloured by sex.
For completeness, we show the MDS by ART and sex (Figure 5). It confirms the trends seen above where the samples are split on the x axis by sex and on the y axis by ART classification.
When sex chromosomes are removed, it is clear that ART classification is the dominant source of variance.
sample_groups <- factor(targets$Sample_Group)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$Sample_Group))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by ART and sex")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
Figure 5. MDS plot coloured by ART classification and sex.
plotMDS(mydist, labels=targets$Sample_Name,col=colours,main="sex chromosomes included")
Figure 5. MDS plot coloured by ART classification and sex.
plotMDS(mydist_flt, labels=targets$Sample_Name,col=colours,main="sex chromosomes excluded")
Figure 5. MDS plot coloured by ART classification and sex.
To acertain whether technical factors like batch effects account for variance on the two main PCs, we create an MDS by array chip (Figure 6). There appears to be no relationship between array chip number and the top two PCs in the two MDS plots shown below.
sample_groups <- factor(targets$Array)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")
colours <- colour_palette[as.integer(factor(targets$Array))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by array chip")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
Figure 6. MDS plot coloured by array chip.
plotMDS(mydist, labels=targets$Sample_Name,col=colours,main="sex chromosomes included")
Figure 6. MDS plot coloured by array chip.
plotMDS(mydist_flt, labels=targets$Sample_Name,col=colours,main="sex chromosomes excluded")
Figure 6. MDS plot coloured by array chip.
This was also performed for sample plate (Figure 7). It appears that the laboratory did not randomise samples over plates 1 and 2. For example the samples on plate 2 appear to be mostly NAT. This does not appear to be a major concern because NAT samples from plate 2 seem to cluster with NAT samples on plate 1 in both MDS plots.
sample_groups <- factor(targets$Sample_Plate)
colour_palette=brewer.pal(n = length(levels(sample_groups)), name = "Paired")## Warning in brewer.pal(n = length(levels(sample_groups)), name = "Paired"): minimal value for n is 3, returning requested palette with 3 different levelscolours <- colour_palette[as.integer(factor(targets$Sample_Plate))]
plot(1,axes = FALSE,xlab="",ylab="",main="MDS by sample plate")
legend("center",legend=levels(sample_groups),pch=16,cex=1.2,col=colour_palette)
Figure 7. MDS plot coloured by sample plate.
plotMDS(mydist, labels=targets$Sample_Name,col=colours,main="sex chromosomes included")
Figure 7. MDS plot coloured by sample plate.
plotMDS(mydist_flt, labels=targets$Sample_Name,col=colours,main="sex chromosomes excluded")
Figure 7. MDS plot coloured by sample plate.
Differential analysis
There are several differential contrasts that would be of interest to us in this study:
NAT (43) vs FH (38): to see what probes are affected by IVF where the embryos are not frozen.
NAT (43) vs FZ (38): to see what probes are affected by IVF where the embryos are frozen.
NAT (43) vs IUI (18): It will be interesting to see if the probes affected by IVF are also sensitive to IUI, a process that involves in utero injection of sperm.
FH (38) vs FZ (38): to see whether freezing results in differential methylation.
NAT (43) vs FH+FZ (76): if freezing results in only a small effect, then these groups can be lumped together and provide more power to determine subtle effects of IVF.
IUI (18) vs FH+FZ (76): It will be interesting to see whether the process of extracting eggs, culturing and fertilising in vitro has an effect that is separate to IUI.
The differential analysis is centred around limma to identify differentially methylated probes. TopConfects was also run to obtain the probes with the largest confident effect (topconfect). There are five outputs below:
Volcano plot (limma result).
Beeswarm plot (top probes by limma p-value).
Heatmap (top probes by limma p-value).
Beeswarm plot (top probes by topconfect ranking).
Heatmap (top probes by topconfect ranking). (TODO)
NAT vs FH
In this first analysis I will look at the effect of natural conception versus IVF (fresh). I have conducted the differential analysis in parallel, both includng and excluding sex chromosomes. Even when excluding autosomes, there are ~2000 probes on the autosomes that show differential methylation between sexes. Based on this data and MDS analysis, we will now exclude sex chromosomes from downstream analysis.
In the comparison of NAT and FH samples there were ~150k probes with differential methylation (FDR<0.05), with ~120k loss of methylation and ~30k with increased methylation.
# include sex chromosomes
samplesheet <- subset(targets,art=="NAT" | art=="FH")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
groups <- factor(samplesheet$art,levels=c("NAT","FH"))
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
top_nat_vs_fh_inc <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval,name="top_nat_vs_fh_inc",
myann=anno ,beta= beta)## Your contrast returned 149589 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Fitting chrX...## Fitting chrY...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.## Tiles plot may use more than one track. Please select correct area for next track if necessary.
head(top_nat_vs_fh_inc$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 328980 cg20757019 Unclassified
## 141894 cg08260549 XPO7;XPO7 Promoter_Associated
## 274682 cg16905586 GMPPB;GMPPB;GMPPB;GMPPB Promoter_Associated
## 35094 cg01931797 KCTD21;USP35 Promoter_Associated
## 135935 cg07897831 CASKIN2;TSEN54 Unclassified
## 293640 cg18233405 TSPYL5;TSPYL5 Promoter_Associated
## Islands_Name Relation_to_Island logFC AveExpr
## 328980 chr2:95872891-95873548 Island -1.2241820 -4.118334
## 141894 chr8:21776842-21777942 Island -0.9930728 -4.160961
## 274682 chr3:49760863-49761895 Island -0.7419176 -3.588790
## 35094 chr11:77899564-77899990 Island -1.0175289 -4.484827
## 135935 chr17:73511016-73513176 Island -0.8928940 -4.323892
## 293640 chr8:98289604-98290404 Island -0.9340364 -4.615953
## t P.Value adj.P.Val B
## 328980 -22.31002 8.326397e-37 3.571491e-31 71.79336
## 141894 -19.02718 5.324395e-32 1.011174e-26 61.44805
## 274682 -18.94747 7.072199e-32 1.011174e-26 61.18019
## 35094 -18.47073 3.924551e-31 3.933862e-26 59.56082
## 135935 -18.40184 5.038632e-31 3.933862e-26 59.32435
## 293640 -18.37758 5.502726e-31 3.933862e-26 59.24095head(top_nat_vs_fh_inc$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr14 101290195-101294147 * | 35 4.82832e-107
## [2] chr11 2016747-2022386 * | 65 2.24831e-105
## [3] chr6 32154116-32167020 * | 142 1.16206e-50
## [4] chr6 31864574-31870840 * | 112 2.77165e-70
## [5] chr6 30709428-30713442 * | 82 5.29446e-115
## [6] chr6 33280571-33291947 * | 208 8.12625e-49
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 3.13237e-132 2.17113e-13 2.82212e-126 -0.0669987 -0.04120846
## [2] 4.21686e-85 5.34869e-08 1.70736e-110 -0.0756624 -0.02474392
## [3] 1.89802e-80 1.10779e-06 1.82761e-107 -0.0960859 -0.01090696
## [4] 7.42090e-97 1.80824e-11 6.79956e-106 -0.1410136 -0.01503305
## [5] 7.49631e-70 9.77437e-18 1.77128e-104 -0.0624100 -0.00831736
## [6] 1.12934e-68 4.58073e-06 2.19136e-102 -0.0552569 -0.00846611
## overlapping.genes
## <character>
## [1] MEG3
## [2] H19
## [3] PBX2, GPSM3, NOTCH4
## [4] C2, EHMT2, ZBTB12
## [5] XXbac-BPG252P9.10, F..
## [6] TAPBP, ZBTB22, DAXX
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengthstop_nat_vs_fh_inc$comp## $up_comp
## all up OR fisherPval
## Intergenic 266413 20788 1.6392397 5.672567e-308
## Gene_Associated 1206 88 1.0952089 4.192755e-01
## Gene_Associated_Cell_type_specific 1767 153 1.3203734 1.626757e-03
## NonGene_Associated 1293 79 0.9049633 4.353968e-01
## NonGene_Associated_Cell_type_specific 190 10 0.7727280 5.600792e-01
## Promoter_Associated 89336 4333 0.6575245 3.412314e-147
## Promoter_Associated_Cell_type_specific 5690 256 0.6522253 1.679106e-12
## Unclassified 29032 1371 0.6739359 9.982092e-49
## Unclassified_Cell_type_specific 34009 1689 0.7099240 5.058297e-44
## lowerCI upperCI
## Intergenic 1.5961167 1.6836860
## Gene_Associated 0.8710422 1.3620318
## Gene_Associated_Cell_type_specific 1.1107794 1.5599922
## NonGene_Associated 0.7111944 1.1372031
## NonGene_Associated_Cell_type_specific 0.3642795 1.4545250
## Promoter_Associated 0.6359030 0.6797499
## Promoter_Associated_Cell_type_specific 0.5728407 0.7399317
## Unclassified 0.6370239 0.7125219
## Unclassified_Cell_type_specific 0.6746072 0.7467107
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 266413 58801 0.4589520 0.000000e+00
## Gene_Associated 1206 299 0.8402879 9.385218e-03
## Gene_Associated_Cell_type_specific 1767 305 0.5308216 1.105283e-26
## NonGene_Associated 1293 474 1.4777826 3.701480e-11
## NonGene_Associated_Cell_type_specific 190 54 1.0125646 9.357899e-01
## Promoter_Associated 89336 36695 2.1168691 0.000000e+00
## Promoter_Associated_Cell_type_specific 5690 1959 1.3445715 3.841120e-25
## Unclassified 29032 10692 1.5339558 7.048211e-241
## Unclassified_Cell_type_specific 34009 11543 1.3430483 7.269885e-130
## lowerCI upperCI
## Intergenic 0.4527888 0.4652341
## Gene_Associated 0.7347243 0.9588405
## Gene_Associated_Cell_type_specific 0.4675606 0.6010473
## NonGene_Associated 1.3166304 1.6571361
## NonGene_Associated_Cell_type_specific 0.7245890 1.3976473
## Promoter_Associated 2.0842897 2.1499081
## Promoter_Associated_Cell_type_specific 1.2718573 1.4211088
## Unclassified 1.4961431 1.5726239
## Unclassified_Cell_type_specific 1.3116893 1.3750390top_nat_vs_fh_inc$cgi## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 146442 6341 0.5248715 0.000000e+00 0.5099721 0.5401447
## N_Shelf 20298 1741 1.3247364 6.783369e-26 1.2586043 1.3937074
## N_Shore 57958 3972 1.0272375 1.290749e-01 0.9919454 1.0635743
## OpenSea 140668 12200 1.5574725 1.309248e-272 1.5198977 1.5958473
## S_Shelf 18006 1556 1.3338608 2.026844e-24 1.2636736 1.4071843
## S_Shore 45564 2957 0.9614608 5.100796e-02 0.9240551 1.0001266
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 146442 59758 2.4998119 0.000000e+00 2.4654871 2.5343182
## N_Shelf 20298 4051 0.6232173 6.772444e-167 0.6015852 0.6454968
## N_Shore 57958 15707 0.9402596 7.430645e-10 0.9218899 0.9589293
## OpenSea 140668 25166 0.4387279 0.000000e+00 0.4318615 0.4456545
## S_Shelf 18006 3555 0.6160495 1.594968e-155 0.5933469 0.6395283
## S_Shore 45564 12585 0.9700307 5.984964e-03 0.9491010 0.9913656names(top_nat_vs_fh_inc)## [1] "dma" "dm_up" "dm_dn" "confects" "dmr" "comp" "cgi"
## [8] "fit"head(top_nat_vs_fh_inc$confects$table$name)## [1] "cg15415945" "cg06834912" "cg06545389" "cg04922029" "cg20757019"
## [6] "cg12697337"# exclude sex chromosomes
samplesheet <- subset(targets,art=="NAT" | art=="FH")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
groups <- factor(samplesheet$art,levels=c("NAT","FH"))
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
top_nat_vs_fh <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_nat_vs_fh",
myann=anno, beta=beta_flt)## Your contrast returned 146555 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.## Tiles plot may use more than one track. Please select correct area for next track if necessary.
head(top_nat_vs_fh$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 321265 cg20757019 Unclassified
## 138436 cg08260549 XPO7;XPO7 Promoter_Associated
## 268273 cg16905586 GMPPB;GMPPB;GMPPB;GMPPB Promoter_Associated
## 34194 cg01931797 KCTD21;USP35 Promoter_Associated
## 132599 cg07897831 CASKIN2;TSEN54 Unclassified
## 286730 cg18233405 TSPYL5;TSPYL5 Promoter_Associated
## Islands_Name Relation_to_Island logFC AveExpr
## 321265 chr2:95872891-95873548 Island -1.2241820 -4.118334
## 138436 chr8:21776842-21777942 Island -0.9930728 -4.160961
## 268273 chr3:49760863-49761895 Island -0.7419176 -3.588790
## 34194 chr11:77899564-77899990 Island -1.0175289 -4.484827
## 132599 chr17:73511016-73513176 Island -0.8928940 -4.323892
## 286730 chr8:98289604-98290404 Island -0.9340364 -4.615953
## t P.Value adj.P.Val B
## 321265 -22.31356 8.337821e-37 3.492155e-31 71.80781
## 138436 -19.03079 5.309772e-32 9.677671e-27 61.46163
## 268273 -18.95590 6.931883e-32 9.677671e-27 61.20999
## 34194 -18.47361 3.921392e-31 3.824710e-26 59.57181
## 132599 -18.40628 5.005783e-31 3.824710e-26 59.34068
## 286730 -18.38140 5.479096e-31 3.824710e-26 59.25513 head(top_nat_vs_fh$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr14 101290195-101294147 * | 35 3.82425e-107
## [2] chr11 2016747-2022386 * | 65 1.73506e-105
## [3] chr6 32154116-32167020 * | 142 1.08956e-50
## [4] chr6 31864574-31870840 * | 112 2.61816e-70
## [5] chr6 30709428-30713442 * | 82 4.92840e-115
## [6] chr6 33280571-33291947 * | 208 7.48769e-49
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 1.93516e-132 2.07706e-13 1.71275e-126 -0.0669987 -0.04120846
## [2] 2.87349e-85 5.23314e-08 1.00142e-110 -0.0756624 -0.02474392
## [3] 1.38988e-80 1.09131e-06 1.19444e-107 -0.0960859 -0.01090696
## [4] 5.40851e-97 1.77345e-11 4.71505e-106 -0.1410136 -0.01503305
## [5] 6.04678e-70 9.55132e-18 1.31151e-104 -0.0624100 -0.00831736
## [6] 8.29144e-69 4.50494e-06 1.39691e-102 -0.0552569 -0.00846611
## overlapping.genes
## <character>
## [1] MEG3
## [2] H19
## [3] PBX2, GPSM3, NOTCH4
## [4] C2, EHMT2, ZBTB12
## [5] XXbac-BPG252P9.10, F..
## [6] TAPBP, ZBTB22, DAXX
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_nat_vs_fh$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 20639 1.6358334 1.283211e-302
## Gene_Associated 1196 87 1.0725853 5.271778e-01
## Gene_Associated_Cell_type_specific 1756 153 1.3063033 2.353077e-03
## NonGene_Associated 1278 78 0.8881848 3.444281e-01
## NonGene_Associated_Cell_type_specific 189 10 0.7635452 4.722967e-01
## Promoter_Associated 87082 4292 0.6571073 3.537608e-146
## Promoter_Associated_Cell_type_specific 5171 243 0.6712047 1.886931e-10
## Unclassified 28583 1370 0.6724366 3.593532e-49
## Unclassified_Cell_type_specific 33088 1680 0.7143489 3.942191e-42
## lowerCI upperCI
## Intergenic 1.5926441 1.6803496
## Gene_Associated 0.8518288 1.3354277
## Gene_Associated_Cell_type_specific 1.0988658 1.5434842
## NonGene_Associated 0.6968896 1.1177612
## NonGene_Associated_Cell_type_specific 0.3599170 1.4374288
## Promoter_Associated 0.6353693 0.6793850
## Promoter_Associated_Cell_type_specific 0.5873456 0.7641008
## Unclassified 0.6356062 0.7110042
## Unclassified_Cell_type_specific 0.6786691 0.7514985
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 57526 0.4586547 0.000000e+00
## Gene_Associated 1196 297 0.8418263 1.000090e-02
## Gene_Associated_Cell_type_specific 1756 305 0.5346665 4.462711e-26
## NonGene_Associated 1278 470 1.4848386 2.837249e-11
## NonGene_Associated_Cell_type_specific 189 54 1.0197463 9.355504e-01
## Promoter_Associated 87082 35804 2.1198081 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 1783 1.3468836 2.769267e-23
## Unclassified 28583 10522 1.5326960 3.614593e-236
## Unclassified_Cell_type_specific 33088 11242 1.3447524 1.429856e-127
## lowerCI upperCI
## Intergenic 0.4524029 0.4650207
## Gene_Associated 0.7357146 0.9610726
## Gene_Associated_Cell_type_specific 0.4709550 0.6054586
## NonGene_Associated 1.3221738 1.6660278
## NonGene_Associated_Cell_type_specific 0.7295098 1.4081089
## Promoter_Associated 2.0868398 2.1532335
## Promoter_Associated_Cell_type_specific 1.2706326 1.4273773
## Unclassified 1.4946375 1.5716444
## Unclassified_Cell_type_specific 1.3129724 1.3772510top_nat_vs_fh$cgi## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 142197 6295 0.5294023 0.000000e+00 0.5142934 0.5448632
## N_Shelf 19854 1728 1.3226360 2.099786e-25 1.2563091 1.3918233
## N_Shore 56505 3929 1.0249248 1.670289e-01 0.9895069 1.0613828
## OpenSea 138314 12128 1.5454977 2.118283e-261 1.5080731 1.5837641
## S_Shelf 17602 1545 1.3332679 3.497647e-24 1.2628101 1.4068861
## S_Shore 44361 2927 0.9616895 5.329602e-02 0.9240699 1.0005860
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 142197 57988 2.4855321 0.000000e+00 2.4511768 2.5205746
## N_Shelf 19854 3971 0.6247484 9.466645e-162 0.6028117 0.6473590
## N_Shore 56505 15410 0.9493474 2.776891e-07 0.9305876 0.9684489
## OpenSea 138314 24804 0.4391975 0.000000e+00 0.4322887 0.4462233
## S_Shelf 17602 3498 0.6210415 1.750459e-147 0.5979318 0.6449091
## S_Shore 44361 12332 0.9794111 6.389525e-02 0.9580349 1.0011932saveRDS(top_nat_vs_fh,file="estill_nat_vs_fh.rds")NAT vs FZ
In the comparison of NAT and FZ samples, there were ~130k probes with lower methylation and ~37k with higher methylation.
samplesheet <- subset(targets,art=="NAT" | art=="FZ")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
groups <- factor(samplesheet$art,levels=c("NAT","FZ"))
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
top_nat_vs_fz <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_nat_vs_fz",
myann=anno, beta=beta_flt)## Your contrast returned 161455 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.## Tiles plot may use more than one track. Please select correct area for next track if necessary.
head(top_nat_vs_fz$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 321265 cg20757019 Unclassified
## 115170 cg06834912 FOXF1
## 51218 cg02916102 CYP11A1;CYP11A1;CYP11A1 Unclassified_Cell_type_specific
## 234473 cg14572967 CD72 Promoter_Associated
## 107155 cg06347499 SLC25A33 Promoter_Associated
## 267523 cg16849041 DOCK5;DOCK5
## Islands_Name Relation_to_Island logFC AveExpr t
## 321265 chr2:95872891-95873548 Island -1.099916 -4.063048 -18.77180
## 115170 chr16:86549069-86550512 N_Shore -3.140598 -4.173200 -18.18771
## 51218 chr15:74658038-74658574 Island -1.036568 -4.419304 -18.13909
## 234473 chr9:35616980-35617324 Island -1.049560 -4.658904 -17.97386
## 107155 chr1:9599350-9600444 Island -1.095213 -4.335596 -17.42775
## 267523 chr8:25041795-25042820 Island -1.126325 -4.502949 -17.29963
## P.Value adj.P.Val B
## 321265 1.406048e-31 5.888995e-26 60.60529
## 115170 1.162966e-30 1.939453e-25 58.60074
## 51218 1.389183e-30 1.939453e-25 58.43184
## 234473 2.546919e-30 2.666834e-25 57.85554
## 107155 1.933929e-29 1.619987e-24 55.92494
## 267523 3.128329e-29 2.183746e-24 55.46621head(top_nat_vs_fz$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr11 2016513-2022386 * | 67 3.80002e-152
## [2] chr14 101290195-101294430 * | 36 5.61629e-120
## [3] chr6 33128825-33139083 * | 99 4.21939e-81
## [4] chr6 32062528-32065890 * | 64 2.16478e-115
## [5] chr20 57425515-57431303 * | 77 1.17336e-133
## [6] chr6 30709428-30713318 * | 81 1.83148e-117
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 1.93630e-122 3.17440e-14 9.14143e-151 -0.0718469 -0.02693608
## [2] 9.16775e-155 6.71799e-10 2.06213e-150 -0.0646443 -0.04581542
## [3] 6.85833e-116 1.38744e-07 1.42506e-137 -0.0699756 -0.01855617
## [4] 3.47984e-136 2.54515e-06 1.42892e-128 -0.0974059 -0.04686334
## [5] 9.19425e-100 8.95856e-09 4.07801e-122 -0.0767107 -0.02790275
## [6] 2.20198e-93 1.01418e-13 1.63684e-116 0.0618130 -0.00591081
## overlapping.genes
## <character>
## [1] H19
## [2] MEG3
## [3] COL11A2
## [4] TNXB
## [5] GNAS, GNAS-AS1
## [6] XXbac-BPG252P9.10, F..
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_nat_vs_fz$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 21503 0.9875438 2.773572e-01
## Gene_Associated 1196 95 0.9543112 7.131338e-01
## Gene_Associated_Cell_type_specific 1756 165 1.1478638 9.916312e-02
## NonGene_Associated 1278 120 1.1467678 1.547791e-01
## NonGene_Associated_Cell_type_specific 189 9 0.5529776 8.539800e-02
## Promoter_Associated 87082 7738 1.1013601 1.260402e-12
## Promoter_Associated_Cell_type_specific 5171 419 0.9750217 6.478420e-01
## Unclassified 28583 2279 0.9554428 4.430339e-02
## Unclassified_Cell_type_specific 33088 2397 0.8538335 2.568840e-13
## lowerCI upperCI
## Intergenic 0.9653861 1.0102278
## Gene_Associated 0.7653295 1.1780135
## Gene_Associated_Cell_type_specific 0.9715042 1.3487630
## NonGene_Associated 0.9419135 1.3854222
## NonGene_Associated_Cell_type_specific 0.2488036 1.0734798
## Promoter_Associated 1.0724741 1.1309321
## Promoter_Associated_Cell_type_specific 0.8796946 1.0783395
## Unclassified 0.9136687 0.9987793
## Unclassified_Cell_type_specific 0.8174698 0.8914776
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 68538 0.6144925 0.000000e+00
## Gene_Associated 1196 322 0.8488056 1.167364e-02
## Gene_Associated_Cell_type_specific 1756 414 0.7100859 4.582624e-10
## NonGene_Associated 1278 463 1.3105926 4.656136e-06
## NonGene_Associated_Cell_type_specific 189 63 1.1525229 3.835920e-01
## Promoter_Associated 87082 33713 1.6212200 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 1738 1.1691989 1.727699e-07
## Unclassified 28583 10065 1.2746400 1.328336e-77
## Unclassified_Cell_type_specific 33088 11414 1.2349504 5.568689e-67
## lowerCI upperCI
## Intergenic 0.6062923 0.6228460
## Gene_Associated 0.7445272 0.9657573
## Gene_Associated_Cell_type_specific 0.6343571 0.7935899
## NonGene_Associated 1.1664940 1.4710053
## NonGene_Associated_Cell_type_specific 0.8377822 1.5717789
## Promoter_Associated 1.5961621 1.6467431
## Promoter_Associated_Cell_type_specific 1.1026356 1.2395112
## Unclassified 1.2427105 1.3072247
## Unclassified_Cell_type_specific 1.2059581 1.2646634saveRDS(top_nat_vs_fz,file="estill_nat_vs_fz.rds")NAT vs IUI
In the comparison of NAT and IUI samples there were ~114k probes with reduced methylation and ~36k with higher methylation.
samplesheet <- subset(targets,art=="NAT" | art=="IUI")
groups <- factor(samplesheet$art,levels=c("NAT","IUI"))
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
table(groups)## groups
## NAT IUI
## 43 18top_nat_vs_iui <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_nat_vs_iui",
myann=anno, beta=beta_flt)## Your contrast returned 147096 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.## Tiles plot may use more than one track. Please select correct area for next track if necessary.
head(top_nat_vs_iui$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 321265 cg20757019 Unclassified
## 307250 cg19732469 PAPOLA Promoter_Associated
## 34194 cg01931797 KCTD21;USP35 Promoter_Associated
## 179618 cg10992925 SEMA3F Unclassified
## 364043 cg23963136 CD320;CD320 Promoter_Associated
## 196471 cg12159575 UHRF1;UHRF1 Promoter_Associated
## Islands_Name Relation_to_Island logFC AveExpr
## 321265 chr2:95872891-95873548 Island -1.2799680 -3.920992
## 307250 chr14:96968154-96969371 Island -0.9791723 -4.548118
## 34194 chr11:77899564-77899990 Island -1.0788716 -4.332283
## 179618 chr3:50191971-50193130 Island -0.8597185 -3.173986
## 364043 chr19:8372597-8373302 Island -1.2766668 -4.599663
## 196471 chr19:4909262-4910256 S_Shore -1.2368248 -4.442626
## t P.Value adj.P.Val B
## 321265 -18.21888 9.859975e-27 4.129683e-21 49.43563
## 307250 -17.75977 3.781829e-26 7.919774e-21 48.19201
## 34194 -16.90450 4.930228e-25 6.883141e-20 45.80467
## 179618 -16.45021 1.996148e-24 2.090131e-19 44.49829
## 364043 -16.02696 7.509535e-24 6.290482e-19 43.25669
## 196471 -15.93245 1.012478e-23 7.067655e-19 42.97618head(top_nat_vs_iui$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr6 33128825-33143276 * | 127 7.94582e-55
## [2] chr11 2016513-2022386 * | 67 2.80308e-84
## [3] chr6 33280436-33292029 * | 214 1.87687e-55
## [4] chr6 30709428-30713442 * | 82 7.92512e-115
## [5] chr6 31864574-31870600 * | 110 1.18371e-69
## [6] chr6 33165404-33173927 * | 110 1.35441e-65
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 2.59977e-91 9.51236e-07 5.03523e-111 -0.0846547 -0.02037885
## [2] 1.38467e-71 3.93570e-06 1.15498e-85 -0.0836684 -0.03000679
## [3] 1.24527e-60 3.73066e-06 9.26277e-85 -0.0662654 -0.01095948
## [4] 1.59370e-58 4.55345e-15 2.50918e-82 -0.0598442 -0.01161075
## [5] 1.91909e-67 1.33225e-08 3.96561e-81 -0.0793614 -0.01432553
## [6] 3.61962e-56 1.66693e-06 3.79943e-69 -0.0550753 -0.00856993
## overlapping.genes
## <character>
## [1] COL11A2
## [2] H19
## [3] TAPBP, ZBTB22, DAXX
## [4] XXbac-BPG252P9.10, F..
## [5] C2, EHMT2, ZBTB12
## [6] RNY4P10, SLC39A7, HS..
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_nat_vs_iui$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 29067 3.2394581 0.000000e+00
## Gene_Associated 1196 116 1.1793388 9.386833e-02
## Gene_Associated_Cell_type_specific 1756 201 1.4209801 8.039071e-06
## NonGene_Associated 1278 35 0.3083256 1.538827e-16
## NonGene_Associated_Cell_type_specific 189 13 0.8105338 5.976200e-01
## Promoter_Associated 87082 2584 0.2825927 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 209 0.4590162 3.709789e-35
## Unclassified 28583 1044 0.3980764 4.320977e-238
## Unclassified_Cell_type_specific 33088 1708 0.5766686 1.919564e-119
## lowerCI upperCI
## Intergenic 3.1476310 3.3345463
## Gene_Associated 0.9649447 1.4298511
## Gene_Associated_Cell_type_specific 1.2201815 1.6476176
## NonGene_Associated 0.2137296 0.4312532
## NonGene_Associated_Cell_type_specific 0.4231138 1.4226732
## Promoter_Associated 0.2712280 0.2943967
## Promoter_Associated_Cell_type_specific 0.3975936 0.5274709
## Unclassified 0.3735100 0.4238807
## Unclassified_Cell_type_specific 0.5482631 0.6062671
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 54454 0.4614319 0.000000e+00
## Gene_Associated 1196 277 0.8241244 4.431569e-03
## Gene_Associated_Cell_type_specific 1756 308 0.5807100 7.801784e-20
## NonGene_Associated 1278 487 1.6872333 9.110923e-19
## NonGene_Associated_Cell_type_specific 189 56 1.1518902 3.670509e-01
## Promoter_Associated 87082 34745 2.1825617 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 1753 1.4094211 5.691658e-30
## Unclassified 28583 9812 1.4711968 3.041747e-187
## Unclassified_Cell_type_specific 33088 10227 1.2462900 2.246206e-68
## lowerCI upperCI
## Intergenic 0.4550495 0.4679104
## Gene_Associated 0.7177698 0.9438103
## Gene_Associated_Cell_type_specific 0.5117415 0.6573093
## NonGene_Associated 1.5037074 1.8917708
## NonGene_Associated_Cell_type_specific 0.8273984 1.5855362
## Promoter_Associated 2.1482773 2.2173044
## Promoter_Associated_Cell_type_specific 1.3292035 1.4940830
## Unclassified 1.4340035 1.5091964
## Unclassified_Cell_type_specific 1.2161194 1.2770908saveRDS(top_nat_vs_iui,file="estill_nat_vs_iui.rds")FH vs FZ
In the comparison of FH and FZ, only ~13k probes were differentially methylated with roughly half lower and higher. This suggests that freezing has only a small effect on methylation.
samplesheet <- subset(targets,art=="FH" | art=="FZ")
groups <- factor(samplesheet$art,levels=c("FH","FZ"))
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
table(groups)## groups
## FH FZ
## 38 38top_fh_vs_fz <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_fh_vs_fz",
myann=anno, beta=beta_flt)## Your contrast returned 12595 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.
head(top_fh_vs_fz$dma)## Row.names UCSC_RefGene_Name
## 216964 cg13601636 FCAR;FCAR;FCAR;FCAR;FCAR;FCAR;FCAR;FCAR;FCAR
## 175512 cg10704137
## 338404 cg22078451 PDE4D;PDE4D;PDE4D
## 201535 cg12538597
## 274520 cg17353079 HSPA1L
## 102728 cg06039171 TUBB1
## Regulatory_Feature_Group Islands_Name Relation_to_Island
## 216964 OpenSea
## 175512 chr6:25732692-25732957 N_Shore
## 338404 chr5:58334836-58335881 Island
## 201535 OpenSea
## 274520 chr6:31782872-31785190 N_Shelf
## 102728 chr20:57599202-57599408 S_Shore
## logFC AveExpr t P.Value adj.P.Val B
## 216964 -0.7068996 1.6085890 -9.316046 2.633176e-14 1.102861e-08 21.52679
## 175512 -0.5647422 1.6851697 -9.060052 8.252249e-14 1.728157e-08 20.48088
## 338404 -0.7742044 -3.3052215 -8.808686 2.537177e-13 3.542178e-08 19.45065
## 201535 -0.5809790 1.3316275 -8.261288 2.931316e-12 3.069330e-07 17.20032
## 274520 -0.5042415 0.8668139 -8.125325 5.379643e-12 4.506344e-07 16.64083
## 102728 -0.3542682 1.6325122 -8.080360 6.575225e-12 4.589869e-07 16.45581head(top_fh_vs_fz$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr6 32941825-32942808 * | 11 2.88566e-30
## [2] chr11 65190180-65191707 * | 11 5.21444e-22
## [3] chr16 3130216-3131034 * | 4 7.60225e-13
## [4] chr6 31777883-31780068 * | 34 3.88689e-29
## [5] chr16 66412840-66413098 * | 3 2.22594e-18
## [6] chr10 13513915-13516215 * | 11 3.57901e-10
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 1.47270e-13 2.44958e-03 8.65895e-13 0.0569789 0.0379123
## [2] 1.25651e-10 5.74394e-04 5.35585e-12 0.1313533 0.0475076
## [3] 1.13321e-06 1.44248e-05 4.35174e-10 -0.0499280 -0.0304902
## [4] 3.75513e-02 1.28434e-05 1.01854e-09 -0.0792462 -0.0129209
## [5] 2.87636e-06 2.60831e-06 1.22513e-08 -0.0777552 -0.0389962
## [6] 1.38163e-09 1.30036e-02 6.19373e-08 -0.0257035 -0.0174928
## overlapping.genes
## <character>
## [1] BRD2
## [2] NEAT1
## [3] IL32, RP11-473M20.9
## [4] HSPA1L
## [5] CDH5
## [6] BEND7
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_fh_vs_fz$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 1791 0.2563434 0.000000e+00
## Gene_Associated 1196 3 0.1739412 7.013617e-05
## Gene_Associated_Cell_type_specific 1756 4 0.1577364 4.685411e-07
## NonGene_Associated 1278 32 1.7848127 2.765940e-03
## NonGene_Associated_Cell_type_specific 189 0 0.0000000 1.210824e-01
## Promoter_Associated 87082 2792 3.4413004 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 121 1.6749497 2.310047e-07
## Unclassified 28583 676 1.7664787 8.488728e-38
## Unclassified_Cell_type_specific 33088 536 1.1556426 1.796700e-03
## lowerCI upperCI
## Intergenic 0.24231341 0.2710916
## Gene_Associated 0.03583583 0.5097680
## Gene_Associated_Cell_type_specific 0.04291445 0.4047831
## NonGene_Associated 1.21428946 2.5357522
## NonGene_Associated_Cell_type_specific 0.00000000 1.3662714
## Promoter_Associated 3.26722049 3.6240907
## Promoter_Associated_Cell_type_specific 1.38434071 2.0100693
## Unclassified 1.62664554 1.9159213
## Unclassified_Cell_type_specific 1.05478641 1.2640720
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 4830 1.6338447 1.507510e-74
## Gene_Associated 1196 8 0.4173025 7.286632e-03
## Gene_Associated_Cell_type_specific 1756 22 0.7868949 2.926549e-01
## NonGene_Associated 1278 22 1.0876249 6.532362e-01
## NonGene_Associated_Cell_type_specific 189 3 1.0012467 1.000000e+00
## Promoter_Associated 87082 1088 0.7433560 5.312671e-20
## Promoter_Associated_Cell_type_specific 5171 62 0.7510103 2.482790e-02
## Unclassified 28583 275 0.5859333 1.072379e-20
## Unclassified_Cell_type_specific 33088 330 0.6057816 2.613178e-21
## lowerCI upperCI
## Intergenic 1.5469736 1.7260205
## Gene_Associated 0.1796774 0.8254741
## Gene_Associated_Cell_type_specific 0.4914770 1.1962803
## NonGene_Associated 0.6786525 1.6555073
## NonGene_Associated_Cell_type_specific 0.2045689 2.9728617
## Promoter_Associated 0.6956573 0.7937564
## Promoter_Associated_Cell_type_specific 0.5742723 0.9656930
## Unclassified 0.5170743 0.6616281
## Unclassified_Cell_type_specific 0.5403432 0.6771528saveRDS(top_fh_vs_fz,file="estill_fh_vs_fz.rds")NAT vs FH/FZ
Because the effect of embryo freezing was found to be very small, we compared the NAT samples to all IVF samples (union of FH and FZ). There were ~200k probes with altered methylation ~157k with lower methylation and ~50k higher. This result points to a major effect of IVF on infant genome methylation.
samplesheet <- subset(targets,art=="NAT" | art=="FH" | art=="FZ")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$art <- gsub("FH","FX",samplesheet$art)
samplesheet$art <- gsub("FZ","FX",samplesheet$art)
groups <- factor(samplesheet$art,levels=c("NAT","FX"))
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
table(groups)## groups
## NAT FX
## 43 76top_nat_vs_fx <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_nat_vs_fx",
myann=anno, beta=beta_flt)## Your contrast returned 200995 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.## Tiles plot may use more than one track. Please select correct area for next track if necessary.
head(top_nat_vs_fx$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 321265 cg20757019 Unclassified
## 26062 cg01447112 Unclassified_Cell_type_specific
## 132599 cg07897831 CASKIN2;TSEN54 Unclassified
## 234473 cg14572967 CD72 Promoter_Associated
## 286730 cg18233405 TSPYL5;TSPYL5 Promoter_Associated
## 51218 cg02916102 CYP11A1;CYP11A1;CYP11A1 Unclassified_Cell_type_specific
## Islands_Name Relation_to_Island logFC AveExpr
## 321265 chr2:95872891-95873548 Island -1.1604770 -4.288225
## 26062 chr7:6703503-6704075 Island -1.0017265 -3.902576
## 132599 chr17:73511016-73513176 Island -0.8790666 -4.464762
## 234473 chr9:35616980-35617324 Island -1.0466745 -4.834420
## 286730 chr8:98289604-98290404 Island -0.9139453 -4.759295
## 51218 chr15:74658038-74658574 Island -1.0096745 -4.578506
## t P.Value adj.P.Val B
## 321265 -23.29733 1.466896e-46 6.143846e-41 94.42966
## 26062 -21.45706 4.525544e-43 7.359759e-38 86.71713
## 132599 -21.42299 5.271618e-43 7.359759e-38 86.57022
## 234473 -21.31572 8.531240e-43 8.932912e-38 86.10668
## 286730 -20.90859 5.371699e-42 4.499689e-37 84.33359
## 51218 -20.13552 1.871178e-40 1.306185e-35 80.90627head(top_nat_vs_fx$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr6 33128825-33143356 * | 128 2.19117e-120
## [2] chr14 101290195-101294430 * | 36 3.09033e-185
## [3] chr11 2016513-2022386 * | 67 1.14392e-205
## [4] chr6 33280052-33291947 * | 223 1.40929e-67
## [5] chr6 32153540-32167900 * | 147 1.21222e-81
## [6] chr6 33165404-33176673 * | 146 4.01795e-109
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 2.06246e-217 8.88682e-14 2.42499e-268 -0.0813805 -0.01761184
## [2] 4.85525e-243 1.40057e-17 1.15816e-242 -0.0642029 -0.04345260
## [3] 6.39933e-176 2.27513e-13 4.25572e-227 -0.0744775 -0.02584884
## [4] 3.73521e-148 1.37018e-09 1.08779e-206 -0.0479303 -0.00666412
## [5] 4.55315e-145 2.55205e-09 1.25607e-204 -0.0866510 -0.00956218
## [6] 1.07999e-138 1.18717e-14 6.72758e-202 -0.0497764 -0.00471404
## overlapping.genes
## <character>
## [1] COL11A2
## [2] MEG3
## [3] H19
## [4] TAPBP, ZBTB22, DAXX
## [5] PBX2, XXbac-BPG300A1..
## [6] RNY4P10, SLC39A7, HS..
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_nat_vs_fx$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 34003 1.6038946 0.000000e+00
## Gene_Associated 1196 141 1.0435597 6.158255e-01
## Gene_Associated_Cell_type_specific 1756 252 1.3097512 1.186881e-04
## NonGene_Associated 1278 121 0.8160349 3.392844e-02
## NonGene_Associated_Cell_type_specific 189 17 0.7715675 3.588797e-01
## Promoter_Associated 87082 7310 0.6637023 1.229075e-223
## Promoter_Associated_Cell_type_specific 5171 444 0.7308016 6.300443e-11
## Unclassified 28583 2389 0.6968157 2.604407e-66
## Unclassified_Cell_type_specific 33088 2879 0.7275512 2.214346e-60
## lowerCI upperCI
## Intergenic 1.5705549 1.6379737
## Gene_Associated 0.8688510 1.2454671
## Gene_Associated_Cell_type_specific 1.1411770 1.4981810
## NonGene_Associated 0.6707851 0.9850479
## NonGene_Associated_Cell_type_specific 0.4392214 1.2720259
## Promoter_Associated 0.6465302 0.6813302
## Promoter_Associated_Cell_type_specific 0.6612224 0.8060439
## Unclassified 0.6672193 0.7274942
## Unclassified_Cell_type_specific 0.6991522 0.7568800
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 79145 0.4937364 0.000000e+00
## Gene_Associated 1196 396 0.8557773 1.156574e-02
## Gene_Associated_Cell_type_specific 1756 462 0.6163830 3.720453e-20
## NonGene_Associated 1278 587 1.4710594 1.082030e-11
## NonGene_Associated_Cell_type_specific 189 74 1.1130359 4.969227e-01
## Promoter_Associated 87082 43386 2.0002008 0.000000e+00
## Promoter_Associated_Cell_type_specific 5171 2253 1.3404234 7.152760e-25
## Unclassified 28583 12961 1.4750897 1.548927e-214
## Unclassified_Cell_type_specific 33088 14175 1.3264643 3.609859e-129
## lowerCI upperCI
## Intergenic 0.4873902 0.5001896
## Gene_Associated 0.7566229 0.9666976
## Gene_Associated_Cell_type_specific 0.5529628 0.6861030
## NonGene_Associated 1.3154074 1.6448694
## NonGene_Associated_Cell_type_specific 0.8195684 1.5036392
## Promoter_Associated 1.9701072 2.0306933
## Promoter_Associated_Cell_type_specific 1.2677133 1.4172051
## Unclassified 1.4397645 1.5114399
## Unclassified_Cell_type_specific 1.2965564 1.3571061saveRDS(top_nat_vs_fx,file="estill_nat_vs_fx.rds")IUI vs FH/FZ
We know that IVF conceived infants have a major shift in genome methylation, but we don’t know exactly whether it is caused by the IVF process itself. To eliminate whether it might be potentially caused by subfertility we looked at whether the differences between IUI and FH+FZ were large.
We observe only ~1500 differentiallymethylated probes which is drastically smaller than the NAT vs FH/FZ comparison which suggests that the process of IVF accounts for only a small fraction (1%) of the variance in methylation. The difference in methylomes could be due to other factors such as subfertility or even the effect of hormone therapy.
samplesheet <- subset(targets,art=="IUI" | art=="FH" | art=="FZ")
samplesheet$Basename <- sapply(strsplit(samplesheet$Basename, "/"), "[[", 3)
samplesheet$art <- gsub("FH","FX",samplesheet$art)
samplesheet$art <- gsub("FZ","FX",samplesheet$art)
groups <- factor(samplesheet$art,levels=c("IUI","FX"))
sex <- factor(samplesheet$sex,levels=c("Male","Female"))
table(groups)## groups
## IUI FX
## 18 76top_iui_vs_fx <- dm_analysis(samplesheet=samplesheet,
sex=sex,groups=groups,mx=Mval_flt,name="top_iui_vs_fx",
myann=anno, beta=beta_flt)## Your contrast returned 2549 individually significant probes. We recommend the default setting of pcutoff in dmrcate().## Fitting chr1...## Fitting chr2...## Fitting chr3...## Fitting chr4...## Fitting chr5...## Fitting chr6...## Fitting chr7...## Fitting chr8...## Fitting chr9...## Fitting chr10...## Fitting chr11...## Fitting chr12...## Fitting chr13...## Fitting chr14...## Fitting chr15...## Fitting chr16...## Fitting chr17...## Fitting chr18...## Fitting chr19...## Fitting chr20...## Fitting chr21...## Fitting chr22...## Demarcating regions...## Done!## snapshotDate(): 2023-10-24## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation## loading from cache
##
## RCircos.Core.Components initialized.
## Type ?RCircos.Reset.Plot.Parameters to see how to modify the core components.
head(top_iui_vs_fx$dma)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group
## 276287 cg17476421 ZNF14 Promoter_Associated
## 171967 cg10471644 TPI1;TPI1;TPI1 Promoter_Associated
## 180034 cg11022944 PGGT1B Promoter_Associated_Cell_type_specific
## 311754 cg20040772 MEPCE;ZCWPW1 Promoter_Associated
## 321952 cg20811730 PRPF3 Promoter_Associated
## 152708 cg09197288 IFI27L1;IFI27L1
## Islands_Name Relation_to_Island logFC AveExpr
## 276287 chr19:19843482-19843943 S_Shore 0.4720359 -4.343799
## 171967 chr12:6976392-6977393 Island 0.4727934 -4.661678
## 180034 chr5:114598301-114598657 Island 0.6629178 -2.998964
## 311754 chr7:100025804-100028002 Island 0.4485672 -4.283353
## 321952 chr1:150293835-150294140 Island 0.5150811 -3.722163
## 152708 chr14:94547001-94547665 S_Shelf -0.4165592 1.780222
## t P.Value adj.P.Val B
## 276287 7.013256 3.242711e-10 0.0001358154 12.40438
## 171967 6.761548 1.059350e-09 0.0001964143 11.35292
## 180034 6.685985 1.507433e-09 0.0001964143 11.03940
## 311754 6.601775 2.230007e-09 0.0001964143 10.69125
## 321952 6.590957 2.344780e-09 0.0001964143 10.64662
## 152708 -6.492829 3.691887e-09 0.0002572800 10.24290head(top_iui_vs_fx$dmr)## GRanges object with 6 ranges and 8 metadata columns:
## seqnames ranges strand | no.cpgs min_smoothed_fdr
## <Rle> <IRanges> <Rle> | <integer> <numeric>
## [1] chr7 94284258-94287211 * | 88 2.12205e-87
## [2] chr11 49872169-49872414 * | 3 5.86804e-13
## [3] chr7 4353644-4353846 * | 2 4.34195e-09
## [4] chr4 175750109-175750802 * | 12 1.69949e-16
## [5] chr9 68298597-68298755 * | 2 7.97554e-09
## [6] chr15 52044314-52044428 * | 2 3.72393e-10
## Stouffer HMFDR Fisher maxdiff meandiff
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 4.77894e-04 0.012357588 6.03171e-08 -0.0549819 -0.0172863
## [2] 4.02800e-05 0.000797675 6.69981e-05 -0.0819196 -0.0503027
## [3] 1.82222e-04 0.003644612 3.50597e-04 -0.0602631 -0.0436041
## [4] 4.95056e-05 0.063028170 3.68171e-04 -0.0477983 -0.0229640
## [5] 1.97167e-03 0.000537735 7.62669e-04 -0.0708783 -0.0220389
## [6] 1.56697e-02 0.000514361 1.65704e-03 -0.0398125 -0.0170898
## overlapping.genes
## <character>
## [1] PEG10, SGCE
## [2] RP11-163O19.1
## [3] <NA>
## [4] GLRA3
## [5] <NA>
## [6] TMOD2
## -------
## seqinfo: 22 sequences from an unspecified genome; no seqlengthstop_iui_vs_fx$comp## $up_comp
## all up OR fisherPval
## Intergenic 260490 432 0.2616399 5.744634e-133
## Gene_Associated 1196 1 0.2435619 2.039970e-01
## Gene_Associated_Cell_type_specific 1756 1 0.1656212 3.673167e-02
## NonGene_Associated 1278 9 2.0754532 4.641704e-02
## NonGene_Associated_Cell_type_specific 189 0 0.0000000 1.000000e+00
## Promoter_Associated 87082 718 3.8599189 5.882006e-134
## Promoter_Associated_Cell_type_specific 5171 32 1.8350377 1.628974e-03
## Unclassified 28583 147 1.5659928 1.213337e-06
## Unclassified_Cell_type_specific 33088 91 0.7911211 3.062632e-02
## lowerCI upperCI
## Intergenic 0.233135635 0.2932149
## Gene_Associated 0.006160167 1.3618369
## Gene_Associated_Cell_type_specific 0.004191269 0.9254721
## NonGene_Associated 0.944851173 3.9636579
## NonGene_Associated_Cell_type_specific 0.000000000 5.7534509
## Promoter_Associated 3.474181905 4.2889913
## Promoter_Associated_Cell_type_specific 1.248846222 2.6060358
## Unclassified 1.310654410 1.8596259
## Unclassified_Cell_type_specific 0.632424271 0.9790893
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 260490 808 1.5861887 1.278279e-12
## Gene_Associated 1196 1 0.3120458 3.890428e-01
## Gene_Associated_Cell_type_specific 1756 5 1.0671745 8.142486e-01
## NonGene_Associated 1278 3 0.8787974 1.000000e+00
## NonGene_Associated_Cell_type_specific 189 0 0.0000000 1.000000e+00
## Promoter_Associated 87082 178 0.7208341 4.088191e-05
## Promoter_Associated_Cell_type_specific 5171 17 1.2359595 3.431049e-01
## Unclassified 28583 51 0.6519513 1.965145e-03
## Unclassified_Cell_type_specific 33088 55 0.6025648 9.828026e-05
## lowerCI upperCI
## Intergenic 1.389596883 1.8142393
## Gene_Associated 0.007889684 1.7454001
## Gene_Associated_Cell_type_specific 0.345400206 2.5025291
## NonGene_Associated 0.180707605 2.5806707
## NonGene_Associated_Cell_type_specific 0.000000000 7.3742216
## Promoter_Associated 0.610521667 0.8470960
## Promoter_Associated_Cell_type_specific 0.716803938 1.9901639
## Unclassified 0.482166194 0.8639093
## Unclassified_Cell_type_specific 0.450774625 0.7905907saveRDS(top_iui_vs_fx,file="estill_iui_vs_fx.rds")References
Estill MS, Bolnick JM, Waterland RA, Bolnick AD, Diamond MP, Krawetz SA. Assisted reproductive technology alters deoxyribonucleic acid methylation profiles in bloodspots of newborn infants. Fertil Steril. 2016 Sep 1;106(3):629-639.e10. doi: 10.1016/j.fertnstert.2016.05.006. Epub 2016 Jun 8. PubMed PMID: 27288894.
Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. 2016 Jan 15;32(2):286-8. doi: 10.1093/bioinformatics/btv560. Epub 2015 Sep 30. PubMed PMID: 26424855.
Sofer T, Schifano ED, Hoppin JA, Hou L, Baccarelli AA. A-clustering: a novel method for the detection of co-regulated methylation regions, and regions associated with exposure. Bioinformatics. 2013 Nov 15;29(22):2884-91. doi: 10.1093/bioinformatics/btt498. Epub 2013 Aug 29. PubMed PMID: 23990415; PubMed Central PMCID: PMC3810849.
Tian Y, Morris TJ, Webster AP, Yang Z, Beck S, Feber A, Teschendorff AE. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics. 2017 Dec 15;33(24):3982-3984. doi: 10.1093/bioinformatics/btx513. PubMed PMID: 28961746; PubMed Central PMCID: PMC5860089.
Session info
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## 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
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] IlluminaHumanMethylationEPICmanifest_0.3.0
## [2] ENmix_1.38.01
## [3] doParallel_1.0.17
## [4] qqman_0.1.9
## [5] RCircos_1.2.2
## [6] beeswarm_0.4.0
## [7] forestplot_3.1.3
## [8] abind_1.4-5
## [9] checkmate_2.3.1
## [10] reshape2_1.4.4
## [11] gplots_3.1.3
## [12] eulerr_7.0.0
## [13] GEOquery_2.70.0
## [14] RColorBrewer_1.1-3
## [15] IlluminaHumanMethylation450kmanifest_0.4.0
## [16] topconfects_1.18.0
## [17] DMRcatedata_2.20.0
## [18] ExperimentHub_2.10.0
## [19] AnnotationHub_3.10.0
## [20] BiocFileCache_2.10.1
## [21] dbplyr_2.4.0
## [22] DMRcate_2.16.1
## [23] limma_3.58.1
## [24] missMethyl_1.36.0
## [25] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [26] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [27] minfi_1.48.0
## [28] bumphunter_1.44.0
## [29] locfit_1.5-9.8
## [30] iterators_1.0.14
## [31] foreach_1.5.2
## [32] Biostrings_2.70.1
## [33] XVector_0.42.0
## [34] SummarizedExperiment_1.32.0
## [35] Biobase_2.62.0
## [36] MatrixGenerics_1.14.0
## [37] matrixStats_1.2.0
## [38] GenomicRanges_1.54.1
## [39] GenomeInfoDb_1.38.2
## [40] IRanges_2.36.0
## [41] S4Vectors_0.40.2
## [42] BiocGenerics_0.48.1
## [43] R.utils_2.12.3
## [44] R.oo_1.25.0
## [45] R.methodsS3_1.8.2
## [46] plyr_1.8.9
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.34.0 bitops_1.0-7
## [3] httr_1.4.7 dynamicTreeCut_1.63-1
## [5] tools_4.3.2 doRNG_1.8.6
## [7] backports_1.4.1 utf8_1.2.4
## [9] R6_2.5.1 HDF5Array_1.30.0
## [11] lazyeval_0.2.2 Gviz_1.46.1
## [13] rhdf5filters_1.14.1 permute_0.9-7
## [15] withr_2.5.2 prettyunits_1.2.0
## [17] gridExtra_2.3 base64_2.0.1
## [19] preprocessCore_1.64.0 cli_3.6.2
## [21] sass_0.4.8 readr_2.1.4
## [23] genefilter_1.84.0 askpass_1.2.0
## [25] Rsamtools_2.18.0 foreign_0.8-85
## [27] siggenes_1.76.0 illuminaio_0.44.0
## [29] dichromat_2.0-0.1 scrime_1.3.5
## [31] BSgenome_1.70.1 readxl_1.4.3
## [33] impute_1.76.0 rstudioapi_0.15.0
## [35] RSQLite_2.3.4 generics_0.1.3
## [37] BiocIO_1.12.0 gtools_3.9.5
## [39] dplyr_1.1.4 Matrix_1.6-1.1
## [41] interp_1.1-5 fansi_1.0.6
## [43] lifecycle_1.0.4 yaml_2.3.8
## [45] edgeR_4.0.3 rhdf5_2.46.1
## [47] SparseArray_1.2.2 blob_1.2.4
## [49] promises_1.2.1 crayon_1.5.2
## [51] lattice_0.22-5 GenomicFeatures_1.54.1
## [53] annotate_1.80.0 KEGGREST_1.42.0
## [55] pillar_1.9.0 knitr_1.45
## [57] beanplot_1.3.1 rjson_0.2.21
## [59] codetools_0.2-19 glue_1.6.2
## [61] data.table_1.14.10 vctrs_0.6.5
## [63] png_0.1-8 cellranger_1.1.0
## [65] gtable_0.3.4 assertthat_0.2.1
## [67] cachem_1.0.8 xfun_0.41
## [69] S4Arrays_1.2.0 mime_0.12
## [71] survival_3.5-7 statmod_1.5.0
## [73] interactiveDisplayBase_1.40.0 ellipsis_0.3.2
## [75] nlme_3.1-163 bit64_4.0.5
## [77] bsseq_1.38.0 progress_1.2.3
## [79] filelock_1.0.3 bslib_0.6.1
## [81] nor1mix_1.3-2 irlba_2.3.5.1
## [83] KernSmooth_2.23-22 rpart_4.1.21
## [85] colorspace_2.1-0 DBI_1.1.3
## [87] Hmisc_5.1-1 nnet_7.3-19
## [89] tidyselect_1.2.0 bit_4.0.5
## [91] compiler_4.3.2 curl_5.2.0
## [93] htmlTable_2.4.2 xml2_1.3.6
## [95] RPMM_1.25 DelayedArray_0.28.0
## [97] rtracklayer_1.62.0 caTools_1.18.2
## [99] scales_1.3.0 quadprog_1.5-8
## [101] rappdirs_0.3.3 stringr_1.5.1
## [103] digest_0.6.33 rmarkdown_2.25
## [105] htmltools_0.5.7 pkgconfig_2.0.3
## [107] jpeg_0.1-10 base64enc_0.1-3
## [109] sparseMatrixStats_1.14.0 highr_0.10
## [111] fastmap_1.1.1 ensembldb_2.26.0
## [113] rlang_1.1.2 htmlwidgets_1.6.4
## [115] shiny_1.8.0 DelayedMatrixStats_1.24.0
## [117] jquerylib_0.1.4 jsonlite_1.8.8
## [119] BiocParallel_1.36.0 mclust_6.0.1
## [121] VariantAnnotation_1.48.1 RCurl_1.98-1.13
## [123] magrittr_2.0.3 Formula_1.2-5
## [125] GenomeInfoDbData_1.2.11 Rhdf5lib_1.24.1
## [127] munsell_0.5.0 Rcpp_1.0.11
## [129] stringi_1.8.3 zlibbioc_1.48.0
## [131] MASS_7.3-60 org.Hs.eg.db_3.18.0
## [133] deldir_2.0-2 splines_4.3.2
## [135] multtest_2.58.0 hms_1.1.3
## [137] rngtools_1.5.2 geneplotter_1.80.0
## [139] biomaRt_2.58.0 BiocVersion_3.18.0
## [141] XML_3.99-0.16 evaluate_0.23
## [143] calibrate_1.7.7 latticeExtra_0.6-30
## [145] biovizBase_1.50.0 BiocManager_1.30.22
## [147] tzdb_0.4.0 httpuv_1.6.13
## [149] tidyr_1.3.0 openssl_2.1.1
## [151] purrr_1.0.2 reshape_0.8.9
## [153] ggplot2_3.4.4 xtable_1.8-4
## [155] restfulr_0.0.15 AnnotationFilter_1.26.0
## [157] later_1.3.2 tibble_3.2.1
## [159] memoise_2.0.1 AnnotationDbi_1.64.1
## [161] GenomicAlignments_1.38.0 cluster_2.1.4END OF REPORT