GOmeth test

Introduction

Here we will run test enrichment analysis approaches for understanding the effect of elevated homocysteine on blood methylomes. Homocysteine is thought to be an inhibitor of methyltransferases and hyperhomocysteine is associated with overall reduced genomic methylation. Previously we analysed the B-PROOF study data which enabled us to identify probes whose methylation correlatesor anti-correlated with homocysteine level.

The two enrichment analysis approaches we will use are:

• gsameth: an ORA approach that addressed the biases inherent with infinium array analysis.

• gmea: an experimental application of functional class scoring for infinium analysis.

suppressPackageStartupMessages({
  library("missMethyl")
  library("kableExtra")
  library("mitch")
  library("eulerr")
  library("IlluminaHumanMethylation450kmanifest")
  library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
  library("tictoc")
})

CORES=12

Load the data.

dm <- read.table("dma3a.tsv")
rownames(dm) <- dm$Row.names
dm$Row.names=NULL

head(dm, 10) %>% kbl() %>% kable_paper("hover", full_width = F)

	UCSC_RefGene_Name	Regulatory_Feature_Group	logFC	AveExpr	t	P.Value	adj.P.Val	B
cg04905210	POLR3D		-0.2926788	-2.8124509	-5.380966	2.0e-07	0.1002013	5.906039
cg09338148		Unclassified	-0.2938597	-1.1554749	-5.114443	8.0e-07	0.1744820	4.880538
cg04247967			-0.2070451	3.1805934	-4.986000	1.5e-06	0.2090211	4.399278
cg06458106			-0.1906572	1.6757406	-4.793141	3.5e-06	0.2413324	3.693145
cg26425904	OCA2	Unclassified	-0.2540917	-4.0195146	-4.761482	4.0e-06	0.2413324	3.579161
cg19590707	SLC16A3;SLC16A3;SLC16A3	Promoter_Associated	-0.2471850	-1.2377016	-4.716409	4.9e-06	0.2413324	3.417844
cg10746622	PRSSL1		-0.1995897	1.5613271	-4.709610	5.1e-06	0.2413324	3.393607
cg09762242	SIPA1;SIPA1	Promoter_Associated	-0.2092570	-2.5239893	-4.707037	5.1e-06	0.2413324	3.384441
cg11257888	RCAN3		0.2775245	0.4917805	4.677660	5.8e-06	0.2413324	3.280067
cg03622371	P2RY6;P2RY6;P2RY6;P2RY6	Unclassified	-0.2422774	-1.6613396	-4.670158	6.0e-06	0.2413324	3.253486

Reactome pathways were downloaded on the 14th Sept 2023 from MsigDB.

gs_entrez <- gmt_import("c2.cp.reactome.v2023.1.Hs.entrez.gmt")

gs_symbols <- gmt_import("c2.cp.reactome.v2023.1.Hs.symbols.gmt")

ORA with GSAMETH - parth of the MissMethyl package

Now perform ORA analysis using gsameth with the top 500 probes.

sigup <- rownames(head(subset(dm,logFC>0),1000))
sigdn <- rownames(head(subset(dm,logFC<0),1000))

tic()
gsaup <- gsameth(sig.cpg=head(sigup,500), all.cpg=rownames(dm), collection=gs_entrez)

## All input CpGs are used for testing.

gsadn <- gsameth(sig.cpg=head(sigdn,500), all.cpg=rownames(dm), collection=gs_entrez)

## All input CpGs are used for testing.

toc()

## 21.951 sec elapsed

gsaup <- gsaup[order(gsaup$P.DE),]
head(gsaup,10) %>% kbl(caption="Top hypermethylated pathways") %>% kable_paper("hover", full_width = F)

Top hypermethylated pathways

	N	DE	P.DE	FDR
REACTOME_INTERLEUKIN_21_SIGNALING	9	4	0.0000411	0.0446841
REACTOME_INTERLEUKIN_20_FAMILY_SIGNALING	26	5	0.0000729	0.0446841
REACTOME_INTERLEUKIN_6_SIGNALING	11	4	0.0000810	0.0446841
REACTOME_TNF_RECEPTOR_SUPERFAMILY_TNFSF_MEMBERS_MEDIATING_NON_CANONICAL_NF_KB_PATHWAY	16	4	0.0002746	0.0938804
REACTOME_TNFR2_NON_CANONICAL_NF_KB_PATHWAY	94	8	0.0003248	0.0938804
REACTOME_INTERFERON_ALPHA_BETA_SIGNALING	64	7	0.0003406	0.0938804
REACTOME_INTERLEUKIN_9_SIGNALING	7	3	0.0005343	0.1262406
REACTOME_SIGNALING_BY_SCF_KIT	43	6	0.0008023	0.1658711
REACTOME_INTERLEUKIN_6_FAMILY_SIGNALING	24	4	0.0010462	0.1742183
REACTOME_CYTOKINE_SIGNALING_IN_IMMUNE_SYSTEM	687	26	0.0010533	0.1742183

gsadn <- gsadn[order(gsadn$P.DE),]
head(gsadn,10) %>% kbl(caption="Top hypomethylated pathways") %>% kable_paper("hover", full_width = F)

Top hypomethylated pathways

	N	DE	P.DE	FDR
REACTOME_SYNTHESIS_OF_BILE_ACIDS_AND_BILE_SALTS_VIA_27_HYDROXYCHOLESTEROL	15	3	0.0017137	1
REACTOME_THROMBIN_SIGNALLING_THROUGH_PROTEINASE_ACTIVATED_RECEPTORS_PARS	32	4	0.0057943	1
REACTOME_BLOOD_GROUP_SYSTEMS_BIOSYNTHESIS	21	3	0.0060334	1
REACTOME_SYNTHESIS_OF_BILE_ACIDS_AND_BILE_SALTS_VIA_7ALPHA_HYDROXYCHOLESTEROL	24	3	0.0060442	1
REACTOME_TRANSPORT_OF_NUCLEOTIDE_SUGARS	8	2	0.0083146	1
REACTOME_PASSIVE_TRANSPORT_BY_AQUAPORINS	13	2	0.0108103	1
REACTOME_LINOLEIC_ACID_LA_METABOLISM	7	2	0.0108640	1
REACTOME_PRE_NOTCH_PROCESSING_IN_GOLGI	18	3	0.0158030	1
REACTOME_CELL_SURFACE_INTERACTIONS_AT_THE_VASCULAR_WALL	131	7	0.0158418	1
REACTOME_G_ALPHA_12_13_SIGNALLING_EVENTS	71	6	0.0180446	1

GMEA

Now with GMEA.

Start with gathering the annotation set.

ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
myann <- data.frame(ann450k[,c("UCSC_RefGene_Name","Regulatory_Feature_Group")])
promoters <- grep("Prom",myann$Regulatory_Feature_Group)

Histogram of limma t-statistics.

dm <- dm[,c("UCSC_RefGene_Name","t")]
hist(dm$t,breaks=seq(from=-6,to=6,by=1)) ; grid()

Aggregate probe t-scores to gene level scores.

The object generated has two parts:

1. df: a simple table of the gene name and the median tvalue

2. gme_res_df: results of 1-sample Wilcox test.

agg <- function(dm,cores=1) {
  gn <- unique(unlist(strsplit( dm$UCSC_RefGene_Name ,";")))
  gnl <- strsplit( dm$UCSC_RefGene_Name ,";")
  gnl <- mclapply(gnl,unique,mc.cores=cores)
  dm$UCSC_RefGene_Name <- gnl
  l <- mclapply(1:nrow(dm), function(i) {
    a <- dm[i,]
    len <- length(a[[1]][[1]])
    tvals <- as.numeric(rep(a[2],len))
    genes <- a[[1]][[1]]
    data.frame(genes,tvals)
  },mc.cores=cores)
  df <- do.call(rbind,l)
  gme_res <- mclapply( 1:length(gn), function(i) {
    g <- gn[i]
    tstats <- df[which(df$genes==g),"tvals"]
    myn <- length(tstats)
    mymean <- mean(tstats)
    mymedian <- median(tstats)
    wtselfcont <- wilcox.test(tstats)
    res <- c("gene"=g,"nprobes"=myn,"mean"=mymean,"median"=mymedian,
      "p-value(sc)"=wtselfcont$p.value)
  } , mc.cores=cores )
  gme_res_df <- do.call(rbind, gme_res)
  rownames(gme_res_df) <- gme_res_df[,1]
  gme_res_df <- gme_res_df[,-1]
  tmp <- apply(gme_res_df,2,as.numeric)
  rownames(tmp) <- rownames(gme_res_df)
  gme_res_df <- as.data.frame(tmp)
  gme_res_df$sig <- -log10(gme_res_df[,4])
  gme_res_df <- gme_res_df[order(-gme_res_df$sig),]
  gme_res_df$`fdr(sc)` <- p.adjust(gme_res_df$`p-value(sc)`)
  out <- list("df"=df,"gme_res_df"=gme_res_df)
  return(out)
}

tic()
dmagg <- agg(dm,cores=CORES)
toc()

## 48.805 sec elapsed

head(dmagg$gme_res_df) %>% kbl(caption="top gmea genes (Wilcox)") %>% kable_paper("hover", full_width = F)

top gmea genes (Wilcox)

	nprobes	mean	median	p-value(sc)	sig	fdr(sc)
TNXB	531	0.3419073	0.4168675	0	13.96422	0e+00
PCDHA1	162	-0.7098305	-0.6269170	0	13.62040	0e+00
NNAT	49	-0.9402128	-0.9282012	0	13.30331	0e+00
PCDHA2	149	-0.7203157	-0.6349524	0	12.58353	0e+00
PCDHA3	141	-0.7202342	-0.6349524	0	11.67327	0e+00
KCNQ1DN	39	-1.6122306	-1.7196789	0	11.43914	1e-07

Does 1-sample t-test work better?

agg2 <- function(dm,cores=1) {
  gn <- unique(unlist(strsplit( dm$UCSC_RefGene_Name ,";")))
  gnl <- strsplit( dm$UCSC_RefGene_Name ,";")
  gnl <- mclapply(gnl,unique,mc.cores=cores)
  dm$UCSC_RefGene_Name <- gnl
  l <- mclapply(1:nrow(dm), function(i) {
    a <- dm[i,]
    len <- length(a[[1]][[1]])
    tvals <- as.numeric(rep(a["t"],len))
    genes <- a[[1]][[1]]
    data.frame(genes,tvals)
  },mc.cores=cores)
  df <- do.call(rbind,l)
  keep <- names(which(table(df$genes)>1))
  df <- df[df$genes %in% keep,]
  gn <- unique(df$genes)
  gme_res <- lapply( 1:length(gn), function(i) {
    g <- gn[i]
    tstats <- df[which(df$genes==g),"tvals"]
    myn <- length(tstats)
    mymean <- mean(tstats)
    mymedian <- median(tstats)
    tselfcont <- t.test(tstats)
    res <- c("gene"=g,"nprobes"=myn,"mean"=mymean,"median"=mymedian,
      "p-value(sc)"=tselfcont$p.value)
  } )
  gme_res_df <- do.call(rbind, gme_res)
  rownames(gme_res_df) <- gme_res_df[,1]
  gme_res_df <- gme_res_df[,-1]
  tmp <- apply(gme_res_df,2,as.numeric)
  rownames(tmp) <- rownames(gme_res_df)
  gme_res_df <- as.data.frame(tmp)
  gme_res_df$sig <- -log10(gme_res_df[,4])
  gme_res_df <- gme_res_df[order(-gme_res_df$sig),]
  gme_res_df$`fdr(sc)` <- p.adjust(gme_res_df$`p-value(sc)`)
  out <- list("df"=df,"gme_res_df"=gme_res_df)
  return(out)
}

tic()
dmagg2 <- agg2(dm,cores=CORES)
toc()

## 91.983 sec elapsed

head(dmagg2$gme_res_df) %>% kbl(caption="top gmea genes (t-test)") %>% kable_paper("hover", full_width = F)

top gmea genes (t-test)

	nprobes	mean	median	p-value(sc)	sig	fdr(sc)
KCNQ1DN	39	-1.6122306	-1.7196789	0	16.49357	0
NNAT	49	-0.9402128	-0.9282012	0	16.23510	0
PCDHA1	162	-0.7098305	-0.6269170	0	15.72873	0
MESTIT1	59	-0.7582612	-0.8084278	0	14.66640	0
PCDHA2	149	-0.7203157	-0.6349524	0	14.41496	0
TNXB	531	0.3419073	0.4168675	0	14.23694	0

This results shows that wilcox test executes faster but 1-sample t-test gives more significant findings.

Check that the distribution is even.

hist(dmagg2$gme_res_df$median,breaks=30) ; grid()

There could be a directional bias, so I will also conduct a 1-sample t-test.

ttenrich <- function(m,genesets,cores=1) {
  res <- mclapply( 1:length(genesets), function(i) {
    gs <- genesets[i]
    name <- names(gs)
    n_members <- length(which(rownames(m) %in% gs[[1]]))
    if ( n_members > 4 ) {
      tstats <- m[which(rownames(m) %in% gs[[1]]),]
      myn <- length(tstats)
      mymean <- mean(tstats)
      mymedian <- median(tstats)
      wt <- t.test(tstats)
      res <- c(name,myn,mymean,mymedian,wt$p.value)
    }
  } , mc.cores = cores)
  res_df <- do.call(rbind, res)
  rownames(res_df) <- res_df[,1]
  res_df <- res_df[,-1]
  colnames(res_df) <- c("n_genes","t_mean","t_median","p.value")
  tmp <- apply(res_df,2,as.numeric)
  rownames(tmp) <- rownames(res_df)
  res_df <- tmp
  res_df <- as.data.frame(res_df)
  res_df <- res_df[order(res_df$p.value),]
  res_df$logp <- -log10(res_df$p.value )
  res_df$fdr <- p.adjust(res_df$p.value,method="fdr")
  res_df[order(abs(res_df$p.value)),]
  return(res_df)
}

tic()
tres <- ttenrich(m=m,genesets=gs_symbols)
toc()

## 1.171 sec elapsed

head(tres,30) %>%
  kbl(caption = "Top significant genesets") %>%
  kable_paper("hover", full_width = F)

Top significant genesets

	n_genes	t_mean	t_median	p.value	logp	fdr
REACTOME_METABOLISM_OF_RNA	457	0.2559801	0.3620978	0.0000000	8.885050	0.0000017
REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE	82	0.5311678	0.4434600	0.0000000	8.651174	0.0000017
REACTOME_EUKARYOTIC_TRANSLATION_INITIATION	87	0.5107556	0.4849303	0.0000000	8.163993	0.0000027
REACTOME_NONSENSE_MEDIATED_DECAY_NMD	83	0.5044973	0.4749432	0.0000000	8.141000	0.0000027
REACTOME_TRANSLATION	197	0.3507717	0.4039695	0.0000000	7.950902	0.0000033
REACTOME_GPCR_LIGAND_BINDING	307	-0.3751365	-0.4780206	0.0000000	7.361843	0.0000107
REACTOME_EUKARYOTIC_TRANSLATION_ELONGATION	68	0.5338105	0.4434600	0.0000001	7.010496	0.0000207
REACTOME_RESPONSE_OF_EIF2AK4_GCN2_TO_AMINO_ACID_DEFICIENCY	73	0.5073936	0.4302307	0.0000004	6.371865	0.0000787
REACTOME_REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS	121	0.3869456	0.4419840	0.0000009	6.062010	0.0001428
REACTOME_INFLUENZA_INFECTION	115	0.3835674	0.3955125	0.0000014	5.851508	0.0002086
REACTOME_CLASS_A_1_RHODOPSIN_LIKE_RECEPTORS	220	-0.3924788	-0.5963885	0.0000029	5.539357	0.0003891
REACTOME_SIGNALING_BY_GPCR	476	-0.2408081	-0.2447987	0.0000095	5.021140	0.0011763
REACTOME_RNA_POLYMERASE_II_TRANSCRIPTION_TERMINATION	41	0.4897926	0.4702110	0.0000170	4.769846	0.0019367
REACTOME_CELL_CELL_COMMUNICATION	92	-0.4780895	-0.3100255	0.0000285	4.544617	0.0030207
REACTOME_NEURONAL_SYSTEM	295	-0.2729275	-0.2529811	0.0000376	4.425230	0.0037113
REACTOME_INSULIN_RECEPTOR_RECYCLING	27	-0.7806304	-0.7712514	0.0000485	4.314184	0.0043863
REACTOME_PEPTIDE_LIGAND_BINDING_RECEPTORS	136	-0.4173369	-0.6719794	0.0000503	4.298299	0.0043863
REACTOME_PROCESSING_OF_CAPPED_INTRON_CONTAINING_PRE_MRNA	186	0.2386616	0.2931526	0.0000791	4.101830	0.0065125
REACTOME_O_LINKED_GLYCOSYLATION	67	-0.5819935	-0.4289824	0.0000919	4.036790	0.0071665
REACTOME_SIGNALING_BY_RECEPTOR_TYROSINE_KINASES	383	-0.2229328	-0.0509501	0.0001091	3.962284	0.0080823
REACTOME_HEMOSTASIS	418	-0.2260932	-0.0610934	0.0001336	3.874296	0.0094261
REACTOME_PROCESSING_OF_CAPPED_INTRONLESS_PRE_MRNA	22	0.6435041	0.5868480	0.0001723	3.763765	0.0116054
REACTOME_REGULATION_OF_MRNA_STABILITY_BY_PROTEINS_THAT_BIND_AU_RICH_ELEMENTS	63	0.3923384	0.5109388	0.0002420	3.616203	0.0153683
REACTOME_CELL_CELL_JUNCTION_ORGANIZATION	43	-0.5828519	-0.6639309	0.0002489	3.604011	0.0153683
REACTOME_SARS_COV_1_MODULATES_HOST_TRANSLATION_MACHINERY	27	0.5844997	0.4419840	0.0003632	3.439908	0.0215278
REACTOME_CELL_JUNCTION_ORGANIZATION	65	-0.4634360	-0.4063725	0.0004446	3.352058	0.0253406
REACTOME_METABOLISM_OF_LIPIDS	504	-0.1616165	-0.0787298	0.0004867	3.312714	0.0266529
REACTOME_COPI_DEPENDENT_GOLGI_TO_ER_RETROGRADE_TRAFFIC	74	0.3638834	0.4431430	0.0005323	3.273882	0.0266529
REACTOME_ACTIVATION_OF_THE_MRNA_UPON_BINDING_OF_THE_CAP_BINDING_COMPLEX_AND_EIFS_AND_SUBSEQUENT_BINDING_TO_43S	42	0.4627260	0.4361074	0.0005479	3.261337	0.0266529
REACTOME_METABOLISM_OF_COFACTORS	15	0.7217364	0.7610752	0.0005508	3.259009	0.0266529

head(tres[order(-abs(tres$t_median)),],30) %>%
  kbl(caption = "Top effect size genesets") %>%
  kable_paper("hover", full_width = F)

Top effect size genesets

	n_genes	t_mean	t_median	p.value	logp	fdr
REACTOME_NR1H2_NR1H3_REGULATE_GENE_EXPRESSION_LINKED_TO_GLUCONEOGENESIS	5	-0.9453970	-1.989966	0.2810426	0.5512279	0.6487618
REACTOME_TERMINAL_PATHWAY_OF_COMPLEMENT	7	-1.0794163	-1.631113	0.0494754	1.3056103	0.3174139
REACTOME_MOLYBDENUM_COFACTOR_BIOSYNTHESIS	5	-1.1199263	-1.547689	0.1871293	0.7278582	0.5591242
REACTOME_LYSINE_CATABOLISM	5	-0.7951470	-1.546757	0.2653436	0.5761914	0.6404547
REACTOME_REMOVAL_OF_AMINOTERMINAL_PROPEPTIDES_FROM_GAMMA_CARBOXYLATED_PROTEINS	8	-1.1685494	-1.526791	0.0307423	1.5122641	0.2484514
REACTOME_GAMMA_CARBOXYLATION_TRANSPORT_AND_AMINO_TERMINAL_CLEAVAGE_OF_PROTEINS	8	-1.1685494	-1.526791	0.0307423	1.5122641	0.2484514
REACTOME_ELEVATION_OF_CYTOSOLIC_CA2_LEVELS	13	-1.0093864	-1.451307	0.0296922	1.5273580	0.2476962
REACTOME_PROSTANOID_LIGAND_RECEPTORS	7	-1.0460252	-1.409368	0.0973073	1.0118544	0.4305287
REACTOME_INTERLEUKIN_21_SIGNALING	7	1.2774037	1.394065	0.0457604	1.3395106	0.3035760
REACTOME_REGULATION_OF_LOCALIZATION_OF_FOXO_TRANSCRIPTION_FACTORS	7	-0.9500163	-1.383797	0.1451895	0.8380649	0.5057727
REACTOME_PLATELET_ADHESION_TO_EXPOSED_COLLAGEN	12	-1.3382173	-1.375614	0.0017449	2.7582262	0.0527747
REACTOME_PROTON_COUPLED_MONOCARBOXYLATE_TRANSPORT	6	-1.0664352	-1.240376	0.0453047	1.3438570	0.3024393
REACTOME_EICOSANOIDS	9	1.0281067	1.238228	0.0289823	1.5378672	0.2476962
REACTOME_LGI_ADAM_INTERACTIONS	11	-0.7742820	-1.233643	0.0418522	1.3782818	0.2939571
REACTOME_ACETYLCHOLINE_NEUROTRANSMITTER_RELEASE_CYCLE	12	-0.7631296	-1.226517	0.0648040	1.1883985	0.3672483
REACTOME_EXTRINSIC_PATHWAY_OF_FIBRIN_CLOT_FORMATION	5	-1.0570508	-1.212044	0.0338318	1.4706749	0.2652842
REACTOME_APOPTOTIC_CLEAVAGE_OF_CELL_ADHESION_PROTEINS	6	-1.0239564	-1.207050	0.0704990	1.1518171	0.3744785
REACTOME_CALCITONIN_LIKE_LIGAND_RECEPTORS	5	-1.1976169	-1.203856	0.0306966	1.5129092	0.2484514
REACTOME_DEFECTIVE_FACTOR_VIII_CAUSES_HEMOPHILIA_A	7	-0.9045173	-1.198755	0.0412643	1.3844255	0.2926014
REACTOME_RHO_GTPASES_ACTIVATE_RHOTEKIN_AND_RHOPHILINS	8	-0.7460229	-1.175043	0.1260596	0.8994242	0.4794900
REACTOME_SIGNALING_BY_RNF43_MUTANTS	6	-0.7226475	-1.155542	0.1543898	0.8113814	0.5200129
REACTOME_NEGATIVE_REGULATION_OF_ACTIVITY_OF_TFAP2_AP_2_FAMILY_TRANSCRIPTION_FACTORS	10	-1.1117187	-1.143387	0.0014320	2.8440523	0.0471610
REACTOME_METAL_SEQUESTRATION_BY_ANTIMICROBIAL_PROTEINS	5	-1.3105588	-1.111227	0.0799229	1.0973286	0.3989707
REACTOME_RETINOID_CYCLE_DISEASE_EVENTS	7	-0.5975461	-1.111026	0.2170089	0.6635225	0.5922784
REACTOME_CARNITINE_METABOLISM	6	-0.8816198	-1.105754	0.2667640	0.5738727	0.6417927
REACTOME_LYSOSPHINGOLIPID_AND_LPA_RECEPTORS	5	-0.8962023	-1.104353	0.3210665	0.4934050	0.6790726
REACTOME_UPTAKE_AND_FUNCTION_OF_DIPHTHERIA_TOXIN	5	-1.2892659	-1.101326	0.1006063	0.9973749	0.4359606
REACTOME_ALPHA_DEFENSINS	5	-0.6589538	-1.086872	0.1861568	0.7301212	0.5584703
REACTOME_NR1H2_NR1H3_REGULATE_GENE_EXPRESSION_LINKED_TO_LIPOGENESIS	8	-0.8077926	-1.049452	0.1609464	0.7933188	0.5265398
REACTOME_MECP2_REGULATES_TRANSCRIPTION_OF_NEURONAL_LIGANDS	5	-0.7466953	-1.046536	0.0795665	1.0992695	0.3989707

sig <- subset(tres,fdr<0.05)
nrow(sig)

## [1] 47

head(sig[order(-abs(sig$t_median)),],30) %>%
  kbl(caption = "Top effect size genesets after FDR filtering") %>%
  kable_paper("hover", full_width = F)

Top effect size genesets after FDR filtering

	n_genes	t_mean	t_median	p.value	logp	fdr
REACTOME_NEGATIVE_REGULATION_OF_ACTIVITY_OF_TFAP2_AP_2_FAMILY_TRANSCRIPTION_FACTORS	10	-1.1117187	-1.1433866	0.0014320	2.844052	0.0471610
REACTOME_ADHERENS_JUNCTIONS_INTERACTIONS	21	-0.7679970	-0.9805554	0.0008673	3.061827	0.0343566
REACTOME_SULFIDE_OXIDATION_TO_SULFATE	5	-0.9406028	-0.9141215	0.0006725	3.172311	0.0284753
REACTOME_INSULIN_RECEPTOR_RECYCLING	27	-0.7806304	-0.7712514	0.0000485	4.314184	0.0043863
REACTOME_ROS_AND_RNS_PRODUCTION_IN_PHAGOCYTES	31	-0.7279597	-0.7623457	0.0012045	2.919180	0.0435319
REACTOME_METABOLISM_OF_COFACTORS	15	0.7217364	0.7610752	0.0005508	3.259009	0.0266529
REACTOME_TRANSFERRIN_ENDOCYTOSIS_AND_RECYCLING	27	-0.5728570	-0.7325968	0.0015060	2.822183	0.0483742
REACTOME_PROCESSING_OF_INTRONLESS_PRE_MRNAS	16	0.7431768	0.7309230	0.0007286	3.137521	0.0299933
REACTOME_PEPTIDE_LIGAND_BINDING_RECEPTORS	136	-0.4173369	-0.6719794	0.0000503	4.298299	0.0043863
REACTOME_PROTEIN_PROTEIN_INTERACTIONS_AT_SYNAPSES	62	-0.5095257	-0.6668842	0.0005575	3.253742	0.0266529
REACTOME_CELL_CELL_JUNCTION_ORGANIZATION	43	-0.5828519	-0.6639309	0.0002489	3.604011	0.0153683
REACTOME_CLASS_A_1_RHODOPSIN_LIKE_RECEPTORS	220	-0.3924788	-0.5963885	0.0000029	5.539357	0.0003891
REACTOME_PROCESSING_OF_CAPPED_INTRONLESS_PRE_MRNA	22	0.6435041	0.5868480	0.0001723	3.763765	0.0116054
REACTOME_AMINO_ACIDS_REGULATE_MTORC1	34	-0.5095605	-0.5654216	0.0015341	2.814137	0.0483742
REACTOME_REGULATION_OF_MRNA_STABILITY_BY_PROTEINS_THAT_BIND_AU_RICH_ELEMENTS	63	0.3923384	0.5109388	0.0002420	3.616203	0.0153683
REACTOME_EUKARYOTIC_TRANSLATION_INITIATION	87	0.5107556	0.4849303	0.0000000	8.163993	0.0000027
REACTOME_GPCR_LIGAND_BINDING	307	-0.3751365	-0.4780206	0.0000000	7.361843	0.0000107
REACTOME_NONSENSE_MEDIATED_DECAY_NMD	83	0.5044973	0.4749432	0.0000000	8.141000	0.0000027
REACTOME_RNA_POLYMERASE_II_TRANSCRIPTION_TERMINATION	41	0.4897926	0.4702110	0.0000170	4.769846	0.0019367
REACTOME_GOLGI_TO_ER_RETROGRADE_TRANSPORT	101	0.3191950	0.4481900	0.0009018	3.044908	0.0343566
REACTOME_DOWNSTREAM_SIGNALING_EVENTS_OF_B_CELL_RECEPTOR_BCR	56	0.3972437	0.4473579	0.0012624	2.898792	0.0435319
REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE	82	0.5311678	0.4434600	0.0000000	8.651174	0.0000017
REACTOME_EUKARYOTIC_TRANSLATION_ELONGATION	68	0.5338105	0.4434600	0.0000001	7.010496	0.0000207
REACTOME_COPI_DEPENDENT_GOLGI_TO_ER_RETROGRADE_TRAFFIC	74	0.3638834	0.4431430	0.0005323	3.273882	0.0266529
REACTOME_REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS	121	0.3869456	0.4419840	0.0000009	6.062010	0.0001428
REACTOME_SARS_COV_1_MODULATES_HOST_TRANSLATION_MACHINERY	27	0.5844997	0.4419840	0.0003632	3.439908	0.0215278
REACTOME_ACTIVATION_OF_THE_MRNA_UPON_BINDING_OF_THE_CAP_BINDING_COMPLEX_AND_EIFS_AND_SUBSEQUENT_BINDING_TO_43S	42	0.4627260	0.4361074	0.0005479	3.261337	0.0266529
REACTOME_RESPONSE_OF_EIF2AK4_GCN2_TO_AMINO_ACID_DEFICIENCY	73	0.5073936	0.4302307	0.0000004	6.371865	0.0000787
REACTOME_O_LINKED_GLYCOSYLATION	67	-0.5819935	-0.4289824	0.0000919	4.036790	0.0071665
REACTOME_SELENOAMINO_ACID_METABOLISM	80	0.3742453	0.4149139	0.0006150	3.211133	0.0279356

plot(tres$t_median,-log10(tres$p.value),pch=19,col="gray")
grid()
points(sig$t_median,-log10(sig$p.value),pch=19,col="red")

tsig <- sig

Agg limma t-test

Determine the median of probe values

gn <- unique(unlist(strsplit( myann$UCSC_RefGene_Name ,";")))
gnl <- strsplit( myann$UCSC_RefGene_Name ,";")
gnl <- mclapply(gnl,unique,mc.cores=CORES)
myann$gnl <- gnl

mymed <- function(g) {
  probes <- rownames(myann[grep(g,myann$gnl),])
  rows <- which(rownames(mx) %in% probes)
  if ( length(rows) > 1 ) {
    b <- mx[rows,]
    med <- apply(b,2,median)
    med <- matrix(med,nrow=1)
    colnames(med) <- colnames(b)
    rownames(med) <- g
    return(med)
  }
}

mx <- readRDS("dm3_mx.Rds")
design <- readRDS("dm3_design.Rds")

med <- mclapply(gn,mymed,mc.cores=CORES)
med <- med[lapply(med,length)>0]
medf <- do.call(rbind,med)

Limma differential methylation

Then do a limma test on median gene methylation. See if it is more robust than other approaches.

fit.reduced <- lmFit(medf,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))

##        (Intercept)   age sexmale groups
## Down          9708   195     725      0
## NotSig        2681 19514   18535  19714
## Up            7325     5     454      0

dm <- topTable(fit.reduced,coef=4, number = Inf)

#write.table(dma3a,file="dma3a.tsv")

head(dm, 30) %>% kbl() %>% kable_paper("hover", full_width = F)

	logFC	AveExpr	t	P.Value	adj.P.Val	B
NHLRC1	-0.2420055	-1.7115786	-4.785634	0.0000036	0.0715680	4.0385303
IGF2	-0.1012726	-0.2004777	-3.973273	0.0001038	0.3526937	1.1059260
FAM153C	-0.2897096	0.7535631	-3.959580	0.0001094	0.3526937	1.0602332
CD177	-0.1651738	0.7673867	-3.930612	0.0001222	0.3526937	0.9639899
SLC4A8	-0.1836541	-0.0749776	-3.893643	0.0001407	0.3526937	0.8420105
MIR431	0.1269209	2.5871245	3.739056	0.0002506	0.3526937	0.3423244
MIR433	0.1297470	2.9560124	3.738265	0.0002514	0.3526937	0.3398124
PQLC2	-0.2045636	-0.7358737	-3.676370	0.0003152	0.3526937	0.1445278
NKG7	-0.1914152	-0.6171371	-3.660996	0.0003333	0.3526937	0.0964470
LZTS2	-0.1874573	-0.0999778	-3.649719	0.0003471	0.3526937	0.0612867
IL21	0.1983379	-1.6311602	3.623768	0.0003812	0.3526937	-0.0192739
LINGO3	-0.1318951	-0.6057878	-3.620456	0.0003857	0.3526937	-0.0295209
LGALS1	-0.1400368	-0.7523987	-3.601113	0.0004134	0.3526937	-0.0892062
GFRA3	0.1837116	0.1192049	3.597819	0.0004183	0.3526937	-0.0993409
FER	-0.1223290	-0.5724057	-3.594106	0.0004239	0.3526937	-0.1107594
RSPH9	-0.1267647	-0.0298640	-3.583882	0.0004397	0.3526937	-0.1421451
BCL9L	0.1456354	0.6773651	3.562159	0.0004750	0.3526937	-0.2085759
ARMC3	-0.1596330	-2.5107850	-3.558377	0.0004814	0.3526937	-0.2201077
C16orf58	-0.1777868	-1.6120849	-3.554632	0.0004878	0.3526937	-0.2315139
ITM2C	-0.1875109	-0.6986916	-3.550803	0.0004945	0.3526937	-0.2431673
TNFRSF6B	-0.2576646	-0.9539851	-3.545419	0.0005040	0.3526937	-0.2595340
CTSB	-0.1642150	-1.8742349	-3.544140	0.0005063	0.3526937	-0.2634194
SLC35E2	0.1040556	-2.7570462	3.534310	0.0005242	0.3526937	-0.2932399
SERINC5	-0.1359664	0.0749186	-3.519834	0.0005516	0.3526937	-0.3370233
IL19	0.1208098	2.0249626	3.518365	0.0005544	0.3526937	-0.3414558
TNFAIP6	-0.2394318	0.5163111	-3.516354	0.0005584	0.3526937	-0.3475251
ZNF853	0.1198794	0.0746964	3.515093	0.0005608	0.3526937	-0.3513302
KIF5C	-0.1367720	0.2252379	-3.511645	0.0005677	0.3526937	-0.3617216
MIR620	0.1831690	1.6362189	3.507591	0.0005758	0.3526937	-0.3739331
ANXA3	-0.1495810	-0.8393991	-3.496924	0.0005977	0.3526937	-0.4060001

1-sample t-test

ttenrich <- function(m,genesets,cores=1) {
  res <- mclapply( 1:length(genesets), function(i) {
    gs <- genesets[i]
    name <- names(gs)
    n_members <- length(which(rownames(m) %in% gs[[1]]))
    if ( n_members > 4 ) {
      tstats <- m[which(rownames(m) %in% gs[[1]]),]
      myn <- length(tstats)
      mymean <- mean(tstats)
      mymedian <- median(tstats)
      wt <- t.test(tstats)
      res <- c(name,myn,mymean,mymedian,wt$p.value)
    }
  } , mc.cores = cores)
  res_df <- do.call(rbind, res)
  rownames(res_df) <- res_df[,1]
  res_df <- res_df[,-1]
  colnames(res_df) <- c("n_genes","t_mean","t_median","p.value")
  tmp <- apply(res_df,2,as.numeric)
  rownames(tmp) <- rownames(res_df)
  res_df <- tmp
  res_df <- as.data.frame(res_df)
  res_df <- res_df[order(res_df$p.value),]
  res_df$logp <- -log10(res_df$p.value )
  res_df$fdr <- p.adjust(res_df$p.value,method="fdr")
  res_df[order(abs(res_df$p.value)),]
  return(res_df)
}

m <- as.data.frame(dm$t)
rownames(m) <- rownames(dm)
colnames(m) <- "t"

tres <- ttenrich(m=m,genesets=gs_symbols,cores=CORES)

tsig <- subset(tres,fdr<0.05)

nrow(tsig)

## [1] 102

nrow(subset(tsig,t_median>0))

## [1] 60

head(sig[order(-tsig$t_median),],20) %>%
  kbl(caption="hypermethylated pathways") %>%
  kable_paper("hover", full_width = F)

hypermethylated pathways

	n_genes	t_mean	t_median	p.value	logp	fdr
NA	NA	NA	NA	NA	NA	NA
REACTOME_PROCESSING_OF_CAPPED_INTRONLESS_PRE_MRNA	22	0.6435041	0.5868480	0.0001723	3.763765	0.0116054
NA.1	NA	NA	NA	NA	NA	NA
REACTOME_RESPONSE_OF_EIF2AK4_GCN2_TO_AMINO_ACID_DEFICIENCY	73	0.5073936	0.4302307	0.0000004	6.371865	0.0000787
NA.2	NA	NA	NA	NA	NA	NA
REACTOME_NONSENSE_MEDIATED_DECAY_NMD	83	0.5044973	0.4749432	0.0000000	8.141000	0.0000027
REACTOME_TRANSLATION	197	0.3507717	0.4039695	0.0000000	7.950902	0.0000033
REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE	82	0.5311678	0.4434600	0.0000000	8.651174	0.0000017
NA.3	NA	NA	NA	NA	NA	NA
REACTOME_PROCESSING_OF_CAPPED_INTRON_CONTAINING_PRE_MRNA	186	0.2386616	0.2931526	0.0000791	4.101830	0.0065125
NA.4	NA	NA	NA	NA	NA	NA
NA.5	NA	NA	NA	NA	NA	NA
REACTOME_NEGATIVE_REGULATION_OF_ACTIVITY_OF_TFAP2_AP_2_FAMILY_TRANSCRIPTION_FACTORS	10	-1.1117187	-1.1433866	0.0014320	2.844052	0.0471610
NA.6	NA	NA	NA	NA	NA	NA
REACTOME_REGULATION_OF_EXPRESSION_OF_SLITS_AND_ROBOS	121	0.3869456	0.4419840	0.0000009	6.062010	0.0001428
REACTOME_INTRA_GOLGI_AND_RETROGRADE_GOLGI_TO_ER_TRAFFIC	151	0.2553093	0.3731771	0.0009041	3.043774	0.0343566
REACTOME_INFLUENZA_INFECTION	115	0.3835674	0.3955125	0.0000014	5.851508	0.0002086
NA.7	NA	NA	NA	NA	NA	NA
REACTOME_DOWNSTREAM_SIGNALING_EVENTS_OF_B_CELL_RECEPTOR_BCR	56	0.3972437	0.4473579	0.0012624	2.898792	0.0435319
REACTOME_SULFIDE_OXIDATION_TO_SULFATE	5	-0.9406028	-0.9141215	0.0006725	3.172311	0.0284753

nrow(subset(tsig,t_median<0))

## [1] 42

head(sig[order(tsig$t_median),],20) %>%
  kbl(caption="hypomethylated pathways") %>%
  kable_paper("hover", full_width = F)

hypomethylated pathways

	n_genes	t_mean	t_median	p.value	logp	fdr
NA	NA	NA	NA	NA	NA	NA
NA.1	NA	NA	NA	NA	NA	NA
NA.2	NA	NA	NA	NA	NA	NA
NA.3	NA	NA	NA	NA	NA	NA
NA.4	NA	NA	NA	NA	NA	NA
NA.5	NA	NA	NA	NA	NA	NA
REACTOME_GOLGI_TO_ER_RETROGRADE_TRANSPORT	101	0.3191950	0.4481900	0.0009018	3.044908	0.0343566
NA.6	NA	NA	NA	NA	NA	NA
NA.7	NA	NA	NA	NA	NA	NA
NA.8	NA	NA	NA	NA	NA	NA
NA.9	NA	NA	NA	NA	NA	NA
REACTOME_CELL_CELL_COMMUNICATION	92	-0.4780895	-0.3100255	0.0000285	4.544617	0.0030207
NA.10	NA	NA	NA	NA	NA	NA
REACTOME_SIGNALING_BY_GPCR	476	-0.2408081	-0.2447987	0.0000095	5.021140	0.0011763
NA.11	NA	NA	NA	NA	NA	NA
NA.12	NA	NA	NA	NA	NA	NA
NA.13	NA	NA	NA	NA	NA	NA
REACTOME_EUKARYOTIC_TRANSLATION_INITIATION	87	0.5107556	0.4849303	0.0000000	8.163993	0.0000027
NA.14	NA	NA	NA	NA	NA	NA
NA.15	NA	NA	NA	NA	NA	NA

Overlap

Use Euler diagrams to find the similarity. As medians are normally distributed about the zero, one would expect equal numbers of up and down- methylated gene sets. Therefore the t-test based approach could be considered better.

msig_up <- subset(msig,s.dist>0)$set
msig_dn <- subset(msig,s.dist<0)$set

tsig_up <- rownames(subset(tsig,t_median>0))
tsig_dn <- rownames(subset(tsig,t_median<0))

v1 <- list("mitch up"=msig_up, "mitch dn"=msig_dn,
  "ttest up"=tsig_up,"ttest dn"=tsig_dn)

plot(euler(v1),quantities = TRUE)

TODO

Split data set into 2 halves to show consistency.

Session information

For reproducibility.

sessionInfo()

## R version 4.3.1 (2023-06-16)
## 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/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.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: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] ENmix_1.36.03                                      
##  [2] doParallel_1.0.17                                  
##  [3] qqman_0.1.9                                        
##  [4] RCircos_1.2.2                                      
##  [5] beeswarm_0.4.0                                     
##  [6] forestplot_3.1.3                                   
##  [7] abind_1.4-5                                        
##  [8] checkmate_2.2.0                                    
##  [9] reshape2_1.4.4                                     
## [10] gplots_3.1.3                                       
## [11] GEOquery_2.68.0                                    
## [12] RColorBrewer_1.1-3                                 
## [13] topconfects_1.16.0                                 
## [14] DMRcatedata_2.18.0                                 
## [15] ExperimentHub_2.8.1                                
## [16] AnnotationHub_3.8.0                                
## [17] BiocFileCache_2.8.0                                
## [18] dbplyr_2.3.3                                       
## [19] DMRcate_2.14.1                                     
## [20] R.utils_2.12.2                                     
## [21] R.oo_1.25.0                                        
## [22] R.methodsS3_1.8.2                                  
## [23] plyr_1.8.8                                         
## [24] dplyr_1.1.3                                        
## [25] limma_3.56.2                                       
## [26] tictoc_1.2                                         
## [27] IlluminaHumanMethylation450kmanifest_0.4.0         
## [28] eulerr_7.0.0                                       
## [29] mitch_1.12.0                                       
## [30] kableExtra_1.3.4                                   
## [31] missMethyl_1.34.0                                  
## [32] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [33] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [34] minfi_1.46.0                                       
## [35] bumphunter_1.42.0                                  
## [36] locfit_1.5-9.8                                     
## [37] iterators_1.0.14                                   
## [38] foreach_1.5.2                                      
## [39] Biostrings_2.68.1                                  
## [40] XVector_0.40.0                                     
## [41] SummarizedExperiment_1.30.2                        
## [42] Biobase_2.60.0                                     
## [43] MatrixGenerics_1.12.3                              
## [44] matrixStats_1.0.0                                  
## [45] GenomicRanges_1.52.0                               
## [46] GenomeInfoDb_1.36.3                                
## [47] IRanges_2.34.1                                     
## [48] S4Vectors_0.38.1                                   
## [49] BiocGenerics_0.46.0                                
## 
## loaded via a namespace (and not attached):
##   [1] DSS_2.48.0                    ProtGenerics_1.32.0          
##   [3] bitops_1.0-7                  httr_1.4.7                   
##   [5] webshot_0.5.5                 dynamicTreeCut_1.63-1        
##   [7] tools_4.3.1                   doRNG_1.8.6                  
##   [9] backports_1.4.1               utf8_1.2.3                   
##  [11] R6_2.5.1                      HDF5Array_1.28.1             
##  [13] lazyeval_0.2.2                Gviz_1.44.1                  
##  [15] permute_0.9-7                 rhdf5filters_1.12.1          
##  [17] prettyunits_1.1.1             GGally_2.1.2                 
##  [19] gridExtra_2.3                 base64_2.0.1                 
##  [21] preprocessCore_1.62.1         cli_3.6.1                    
##  [23] sass_0.4.7                    readr_2.1.4                  
##  [25] genefilter_1.82.1             askpass_1.2.0                
##  [27] Rsamtools_2.16.0              systemfonts_1.0.4            
##  [29] foreign_0.8-85                siggenes_1.74.0              
##  [31] illuminaio_0.42.0             svglite_2.1.1                
##  [33] dichromat_2.0-0.1             scrime_1.3.5                 
##  [35] BSgenome_1.68.0               readxl_1.4.3                 
##  [37] impute_1.74.1                 rstudioapi_0.15.0            
##  [39] RSQLite_2.3.1                 generics_0.1.3               
##  [41] BiocIO_1.10.0                 gtools_3.9.4                 
##  [43] interp_1.1-4                  Matrix_1.6-1                 
##  [45] fansi_1.0.4                   lifecycle_1.0.3              
##  [47] edgeR_3.42.4                  yaml_2.3.7                   
##  [49] rhdf5_2.44.0                  blob_1.2.4                   
##  [51] promises_1.2.1                crayon_1.5.2                 
##  [53] lattice_0.21-8                echarts4r_0.4.5              
##  [55] GenomicFeatures_1.52.2        annotate_1.78.0              
##  [57] KEGGREST_1.40.0               pillar_1.9.0                 
##  [59] knitr_1.44                    beanplot_1.3.1               
##  [61] rjson_0.2.21                  codetools_0.2-19             
##  [63] glue_1.6.2                    data.table_1.14.8            
##  [65] vctrs_0.6.3                   png_0.1-8                    
##  [67] cellranger_1.1.0              gtable_0.3.4                 
##  [69] assertthat_0.2.1              cachem_1.0.8                 
##  [71] xfun_0.40                     S4Arrays_1.0.6               
##  [73] mime_0.12                     survival_3.5-7               
##  [75] statmod_1.5.0                 interactiveDisplayBase_1.38.0
##  [77] ellipsis_0.3.2                nlme_3.1-163                 
##  [79] bit64_4.0.5                   bsseq_1.36.0                 
##  [81] progress_1.2.2                filelock_1.0.2               
##  [83] bslib_0.5.1                   nor1mix_1.3-0                
##  [85] KernSmooth_2.23-22            rpart_4.1.19                 
##  [87] colorspace_2.1-0              DBI_1.1.3                    
##  [89] Hmisc_5.1-1                   nnet_7.3-19                  
##  [91] tidyselect_1.2.0              bit_4.0.5                    
##  [93] compiler_4.3.1                curl_5.0.2                   
##  [95] rvest_1.0.3                   htmlTable_2.4.1              
##  [97] BiasedUrn_2.0.11              xml2_1.3.5                   
##  [99] RPMM_1.25                     DelayedArray_0.26.7          
## [101] rtracklayer_1.60.1            scales_1.2.1                 
## [103] caTools_1.18.2                quadprog_1.5-8               
## [105] rappdirs_0.3.3                stringr_1.5.0                
## [107] digest_0.6.33                 rmarkdown_2.24               
## [109] htmltools_0.5.6               pkgconfig_2.0.3              
## [111] jpeg_0.1-10                   base64enc_0.1-3              
## [113] sparseMatrixStats_1.12.2      highr_0.10                   
## [115] fastmap_1.1.1                 ensembldb_2.24.0             
## [117] rlang_1.1.1                   htmlwidgets_1.6.2            
## [119] shiny_1.7.5                   DelayedMatrixStats_1.22.6    
## [121] jquerylib_0.1.4               jsonlite_1.8.7               
## [123] BiocParallel_1.34.2           mclust_6.0.0                 
## [125] VariantAnnotation_1.46.0      RCurl_1.98-1.12              
## [127] magrittr_2.0.3                Formula_1.2-5                
## [129] GenomeInfoDbData_1.2.10       Rhdf5lib_1.22.1              
## [131] munsell_0.5.0                 Rcpp_1.0.11                  
## [133] stringi_1.7.12                zlibbioc_1.46.0              
## [135] MASS_7.3-60                   org.Hs.eg.db_3.17.0          
## [137] deldir_1.0-9                  splines_4.3.1                
## [139] multtest_2.56.0               hms_1.1.3                    
## [141] polylabelr_0.2.0              rngtools_1.5.2               
## [143] geneplotter_1.78.0            biomaRt_2.56.1               
## [145] BiocVersion_3.17.1            XML_3.99-0.14                
## [147] evaluate_0.21                 calibrate_1.7.7              
## [149] latticeExtra_0.6-30           biovizBase_1.48.0            
## [151] BiocManager_1.30.22           tzdb_0.4.0                   
## [153] httpuv_1.6.11                 polyclip_1.10-4              
## [155] tidyr_1.3.0                   openssl_2.1.0                
## [157] purrr_1.0.2                   reshape_0.8.9                
## [159] ggplot2_3.4.3                 xtable_1.8-4                 
## [161] restfulr_0.0.15               AnnotationFilter_1.24.0      
## [163] later_1.3.1                   viridisLite_0.4.2            
## [165] tibble_3.2.1                  memoise_2.0.1                
## [167] AnnotationDbi_1.62.2          GenomicAlignments_1.36.0     
## [169] cluster_2.1.4