Replicating a published CTD analysis: Yang et al. (2024)

Purpose

ctdR introduces no statistical method: it applies established tests to the chemical-gene interactions of the Comparative Toxicogenomics Database (CTD). What it must show, then, is that it builds the same gene sets as CTD and, given the same CTD release and the same test, returns the same results as a published analysis. Anything else would validate CTD or the test, not the package.

This vignette does that for a recent study, and then reads the same genes the other way round, screening the whole CTD catalogue. The CTD steps take about ten minutes and a 35 MB download, so they are shown with eval = FALSE; the tables are the results of running them, shipped in inst/extdata/yang2024/, and the complete script is inst/scripts/replicate_yang2024.R. These tables hold counts and p-values only, no CTD records.

The study

Yang et al. (2024) asked whether developmental dyslexia genes are targeted by environmental pollutants (doi:10.1186/s12888-024-05952-4, CC BY 4.0). Their design:

  • Genes: dyslexia risk genes from DisGeNET, OMIM and GeneCards; the text says 131, their Table S2 lists 130.
  • Chemicals: 95 chosen in advance from US EPA and NPIC lists, in seven classes (metals, persistent organic pollutants, PAHs, herbicides, insecticides, fungicides, rodenticides).
  • CTD: “2023 data”, interactions of every species, retrieved by uploading the gene symbols to CTD.
  • Test: for each chemical, observed against expected dyslexia genes, hypergeometric test with a fixed universe of 19,000 protein-coding genes, FDR correction within each class.

Uploading the gene symbols and downloading the results is CTD’s Batch Query, so the set sizes and overlaps in the paper’s tables are the output of the Batch Query on CTD of 2023: comparing them with ctdR compares the two tools on the same period of CTD.

This is verification: the chemicals are chosen first, then tested. It is the same pattern as the CTD Batch Query or CTDquerier.

The same CTD release

CTD distributes only its current release and overwrites the file in place. The release of 23 May 2023, the closest to the authors’ “2023 data”, survives as an Internet Archive capture; its SHA-1 matches the digest the archive publishes for that capture.

library(ctdR)
ctd_file <- "CTD_chem_gene_ixns_2023-05-23.csv.gz"
utils::download.file(paste0("https://web.archive.org/web/20230607020439id_/",
    "https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz"), ctd_file,
    mode = "wb")

The authors pooled species, so the import keeps them all. With organism = "all" the symbol-keyed gene sets are forced to pool species by gene symbol, as the CTD website does: a mouse Esr1 interaction counts as ESR1. This is orthology by homonymy, not an orthology mapping.

import_CTD(ctd_file, organism = "all")
ctd_provenance()
#> CTD provenance
#>   Report created: Tue May 23 14:45:41 EDT 2023
#>   Retained:       14,572 chemicals, 1,523,407 chemical-gene pairs
#>   Organism:       all species

Chemicals

The paper’s tables have 95 entries (hexachlorobenzene appears both as a persistent organic pollutant and as a fungicide). Of these, 72 match a CTD name exactly and 22 match a CTD synonym in the chemical vocabulary of January 2023 (for example Chlordane is CTD’s Chlordan, Benzo(j,k)fluorene is fluoranthene; lambda-cyhalothrin maps to cyhalothrin, CTD not distinguishing the isomer). Metam is not in that vocabulary and is left out.

table(rep$Mapping, useNA = "ifany")
#> 
#>                                         CTD synonym 
#>                                                  21 
#> CTD synonym (cyhalothrin, isomer not distinguished) 
#>                                                   1 
#>                                          exact name 
#>                                                  72 
#>                                          not mapped 
#>                                                   1

Does ctdR build the same gene sets?

For each chemical, N is the number of genes it affects and k how many of them are dyslexia genes. Both feed the test directly.

m <- rep[!is.na(rep$ChemicalID) & nzchar(rep$ChemicalID), ]
c(chemicals = nrow(m),
  N_identical = sum(m$N_ctdR == m$N_paper, na.rm = TRUE),
  N_within_2pct = sum(abs(m$N_ctdR / m$N_paper - 1) <= 0.02, na.rm = TRUE),
  k_identical = sum(m$k_ctdR == m$k_paper, na.rm = TRUE))
#>     chemicals   N_identical N_within_2pct   k_identical 
#>            94            48            67            76
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2, 1))
plot(m$N_paper, m$N_ctdR, log = "xy", pch = 19, col = "#1f77b4",
    xlab = "N, Yang et al.", ylab = "N, ctdR", main = "Set size")
abline(0, 1, lty = 2)
plot(m$k_paper, m$k_ctdR, pch = 19, col = "#1f77b4",
    xlab = "k, Yang et al.", ylab = "k, ctdR", main = "Dyslexia genes in set")
abline(0, 1, lty = 2)

par(op)

Points on the diagonal are the same sets. The points well below it are parent terms and families (polychlorinated biphenyls, cobalt, nickel, chromium, DDT): the authors counted the chemical together with its descendants in the MeSH hierarchy, as CTD’s web pages do, while ctdR uses the exact ChemicalID. Small differences of a few genes come from the release: the paper reports 14,489 chemicals, the release of May 2023 has 14,572.

The same test on ctdR’s sets

The authors’ p-values are reproduced by the hypergeometric test with a universe of 19,000 and the 130 genes of Table S2. Their “Expected” column, however, is computed with 142 genes, which does not match either number. Applying their test to the sets ctdR builds, on the chemicals whose sets agree (N within 2%, same k, k > 0):

ok <- !is.na(m$PValue_paper) & m$k_paper > 0 &
    abs(m$N_ctdR / m$N_paper - 1) <= 0.02 & m$k_ctdR == m$k_paper
cmp <- m[ok, ]
c(chemicals = nrow(cmp),
  median_abs_log10_ratio = round(median(abs(log10(
      cmp$PValue_sameTest_ctdR / cmp$PValue_paper))), 3),
  same_call_at_0.05 = sum((cmp$PValue_sameTest_ctdR < 0.05) ==
      (cmp$PValue_paper < 0.05)))
#>              chemicals median_abs_log10_ratio      same_call_at_0.05 
#>                 51.000                  0.016                 49.000
show <- cmp[order(cmp$PValue_paper), c("Class", "Chemical", "N_paper",
    "N_ctdR", "k_paper", "PValue_paper", "PValueFDR_paper",
    "PValue_sameTest_ctdR")]
show[, 6:8] <- lapply(show[, 6:8], fmt)
knitr::kable(show, row.names = FALSE, col.names = c("Class", "Chemical",
    "N paper", "N ctdR", "k", "p paper", "FDR paper", "p ctdR sets"))
Class Chemical N paper N ctdR k p paper FDR paper p ctdR sets
Herbicides Cyanazine 12 12 4 1.48e-07 1.18e-06 9.93e-07
POPs Dieldrin 1494 1494 25 2.1e-05 0.000202 2.43e-05
PAHs Phenanthrene 154 155 7 7.21e-05 0.001 9.59e-05
Insecticides Lambda-cyhalothrin 210 210 8 8.04e-05 0.000241 9.99e-05
PAHs Anthracene 21 21 3 0.000242 0.002 0.00038
Metals Mercury 635 633 13 0.00037 0.002 0.000406
Fungicides Cycloheximide 501 503 11 0.001 0.003 0.000686
Herbicides Glyphosate 1831 1842 25 0.001 0.002 0.000655
Herbicides Metolachlor 72 72 4 0.001 0.003 0.0015
POPs Toxaphene 80 80 4 0.002 0.006 0.00221
Rodenticides Brodifacoum 660 660 12 0.002 0.005 0.00193
Metals Zinc 2678 2677 30 0.004 0.008 0.00384
Herbicides Trifluralin 18 18 2 0.005 0.01 0.00662
Herbicides Pendimethalin 21 21 2 0.007 0.01 0.00896
PAHs Dibenzo(a, h)pyrene 2 2 1 0.007 0.037 0.0136
Fungicides Benomyl 62 62 3 0.008 0.029 0.00882
Metals Arsenic 4977 4977 46 0.012 0.022 0.0127
Fungicides Ferbam 3 3 1 0.014 0.038 0.0204
POPs Chlordane 215 214 5 0.015 0.037 0.0159
PAHs 7,12-Dimethylbenz(a)anthracene 1304 1305 16 0.016 0.063 0.0165
Metals Barium 81 81 3 0.017 0.027 0.0181
Fungicides Thiabendazole 337 336 6 0.027 0.06 0.0282
PAHs Benzo(j, k)fluorene 956 956 12 0.03 0.08 0.0309
PAHs Indeno(1,2,3-cd)pyrene 255 255 5 0.03 0.08 0.031
Metals Iron 458 458 7 0.037 0.054 0.0382
Fungicides Mancozeb 121 121 3 0.048 0.072 0.0503
PAHs Benzo(g, h,i)perylene 290 290 5 0.048 0.107 0.0493
PAHs Pyrene 300 299 5 0.054 0.107 0.0548
Herbicides 2,4-D 153 152 3 0.085 0.086 0.0864
Herbicides Dicamba 14 14 1 0.086 0.086 0.0917
Rodenticides Bromadiolone 14 14 1 0.086 0.115 0.0917
POPs Polychlorinated dibenzofurans 75 75 2 0.09 0.121 0.0934
POPs Endrin 86 86 2 0.114 0.137 0.117
Insecticides Carbofuran 120 120 2 0.196 0.293 0.198
POPs Aldrin 225 225 3 0.198 0.216 0.2
Metals Aluminum 1937 1932 16 0.249 0.266 0.246
PAHs Benzo(k)fluoranthene 753 752 7 0.256 0.411 0.256
PAHs Benzo(a)anthracene 1413 1413 12 0.259 0.411 0.26
PAHs Dibenzo(a, l)pyrene 394 393 4 0.282 0.411 0.282
PAHs Benzo(a)phenanthrene 954 953 8 0.328 0.438 0.328
Fungicides Chlorothalonil 93 92 1 0.472 0.519 0.469
PAHs Benzo(b)fluoranthene 1684 1683 12 0.487 0.599 0.485
Rodenticides Warfarin 748 748 5 0.585 0.585 0.584
PAHs Dibenzo(a, h)anthracene 916 916 6 0.603 0.689 0.602
Fungicides Ziram 157 157 1 0.662 0.662 0.661
PAHs 3-Methylcholanthrene 1305 1304 8 0.678 0.724 0.676
Metals Vanadium 182 182 1 0.716 0.716 0.715
Insecticides Methyl parathion 101 101 5 0.8 0.997 0.00065
PAHs 1-Nitropyrene 263 264 1 0.839 0.839 0.839
Insecticides Terbufos 733 733 1 0.994 0.997 0.994
Insecticides Aldicarb 810 810 1 0.997 0.997 0.997

The two calls that differ at 0.05:

  • Mancozeb: 0.048 in the paper, 0.050 here, on the threshold.
  • Methyl parathion: the paper reports p = 0.800 with 5 dyslexia genes out of 101 against 0.75 expected; its own numbers give 6.5e-4, so the published value is a transcription error.

The same genes, in discovery mode

ctdR asks a different question: not whether chosen chemicals are enriched, but which chemicals, among all those in CTD, are. The same 130 genes were screened against every chemical of the May 2023 release with ctdR‘s ORA. The universe is restricted to human protein-coding genes, the closest to the authors’ 19,000 that ctdR can apply (it narrows the background, it does not widen it), and the Benjamini-Hochberg correction runs over every chemical tested.

entrez <- AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db,
    genes$GeneSymbol, "ENTREZID", "SYMBOL")
coding <- unique(AnnotationDbi::select(org.Hs.eg.db::org.Hs.eg.db,
    keys = "protein-coding", keytype = "GENETYPE", columns = "SYMBOL")$SYMBOL)
discovery <- enrichment_CTD(
    data.frame(EntrezID = unname(entrez[!is.na(entrez)]), pvalue = 0.01),
    method = "ORA", universe = coding, minGSSize = 2)
knitr::kable(summ, col.names = c("Analysis", "Chemicals tested",
    "Adjusted p < 0.05", "Unadjusted p < 0.05"))
Analysis Chemicals tested Adjusted p < 0.05 Unadjusted p < 0.05
main 10776 640 1871
all_human 10783 1039 2295
min10 4540 658 1136

main is the analysis above; all_human uses every human gene symbol as the universe; min10 requires at least ten genes per set.

The 30 chemicals at the top, with unadjusted and adjusted p-values:

t30 <- top
t30$PValue <- fmt(t30$PValue)
t30$PValueAdjusted <- fmt(t30$PValueAdjusted)
knitr::kable(t30[, c("Rank", "ChemicalName", "Count", "GeneSetSize",
    "FoldEnrichment", "PValue", "PValueAdjusted", "InYangList")],
    col.names = c("Rank", "Chemical", "k", "N", "Fold", "p",
        "Adjusted p", "In Yang's list"))
Rank Chemical k N Fold p Adjusted p In Yang’s list
1 Oxidopamine 15 199 11.30 4.82e-12 5.19e-08 FALSE
2 Levodopa 11 113 14.59 2.46e-10 1.33e-06 FALSE
3 Ethanol 81 6936 1.75 3.7e-10 1.33e-06 FALSE
4 Clozapine 22 751 4.39 4.8e-09 1.13e-05 FALSE
5 1-bromopropane 8 59 20.32 5.25e-09 1.13e-05 FALSE
6 Methamphetamine 57 4339 1.97 3.03e-08 5.45e-05 FALSE
7 Rotenone 55 4193 1.97 6.81e-08 0.000105 FALSE
8 Aluminum Chloride 10 160 9.37 1.2e-07 0.000161 FALSE
9 1-Methyl-4-phenyl-1,2,3,6-tetrahydropyridine 14 361 5.81 1.36e-07 0.000163 FALSE
10 Melatonin 14 368 5.70 1.72e-07 0.000171 FALSE
11 delta-hexachlorocyclohexane 5 20 37.47 1.75e-07 0.000171 FALSE
12 manganese chloride 24 1101 3.27 2.46e-07 0.000182 FALSE
13 Permethrin 28 1454 2.89 2.67e-07 0.000182 TRUE
14 CL 218872 3 3 149.89 2.9e-07 0.000182 FALSE
15 3-N-methylspiperone 3 3 149.89 2.9e-07 0.000182 FALSE
16 (1-(2-(bis(4-fluorophenyl)methoxy)ethyl)-4-(3-hydroxy-3-phenylpropyl) piperazinyl decanoate) 3 3 149.89 2.9e-07 0.000182 FALSE
17 Caffeine 50 3774 1.99 2.99e-07 0.000182 FALSE
18 2-(2-amino-3-methoxyphenyl)-4H-1-benzopyran-4-one 15 448 5.02 3.14e-07 0.000182 FALSE
19 methoxyacetic acid 12 275 6.54 3.21e-07 0.000182 FALSE
20 Vorinostat 40 2673 2.24 3.7e-07 0.000199 FALSE
21 deguelin 20 808 3.71 3.98e-07 0.000204 FALSE
22 Ammonium Chloride 56 4549 1.85 4.48e-07 0.000219 FALSE
23 gamma-Aminobutyric Acid 8 106 11.31 5.49e-07 0.000257 FALSE
24 liquiritigenin 5 25 29.98 5.83e-07 0.000259 FALSE
25 cyanazine 4 11 54.51 6.02e-07 0.000259 TRUE
26 arsenite 54 4353 1.86 6.6e-07 0.000262 FALSE
27 Methylnitronitrosoguanidine 12 295 6.10 6.77e-07 0.000262 FALSE
28 nonylphenol 16 544 4.41 6.99e-07 0.000262 FALSE
29 DEET 23 1082 3.19 7.06e-07 0.000262 FALSE
30 Dopamine 11 244 6.76 7.39e-07 0.000265 FALSE

Most of the chemicals at the top were not in the authors’ list, and most act on the nervous system: oxidopamine, levodopa, clozapine, methamphetamine, MPTP, melatonin, as well as the pesticide rotenone and the solvent 1-bromopropane. This fits the gene list, which contains dopaminergic and glutamatergic genes (DRD2, COMT, GRIN2B, GABRA1), and it is what a list of pollutants fixed in advance cannot show. Manganese chloride and aluminum chloride appear on their own: the MeSH hierarchy again, seen from the other side.

Where the 35 chemicals the authors found significant fall in the full ranking:

p35 <- pos[pos$Analysis == "main", c("Class", "Chemical", "Rank", "Count",
    "GeneSetSize", "PValue", "PValueAdjusted")]
p35$PValue <- fmt(p35$PValue)
p35$PValueAdjusted <- fmt(p35$PValueAdjusted)
knitr::kable(p35, row.names = FALSE, col.names = c("Class", "Chemical",
    "Rank", "k", "N", "p", "Adjusted p"))
Class Chemical Rank k N p Adjusted p
Insecticides Permethrin 13 28 1454 2.67e-07 0.000182
Herbicides Cyanazine 25 4 11 6.02e-07 0.000259
POPs Dieldrin 51 25 1321 1.81e-06 0.000383
POPs Dichlorodiphenyl trichloroethane 53 16 596 2.29e-06 0.000463
Herbicides Atrazine 74 74 7383 6.92e-06 0.00101
Insecticides Lambda-cyhalothrin 90 8 165 1.51e-05 0.0018
POPs Heptachlor 96 7 122 1.77e-05 0.00195
Metals Copper 107 53 4743 2.35e-05 0.00236
Metals Lead 125 37 2872 3.91e-05 0.00292
PAHs Phenanthrene 127 7 139 4.1e-05 0.00292
Herbicides Glyphosate 185 25 1642 7.49e-05 0.00434
Metals Mercury 207 13 571 0.000116 0.00592
Metals Manganese 250 12 510 0.000159 0.00685
Metals Zinc 277 30 2325 0.000257 0.00953
PAHs Anthracene 315 3 21 0.000353 0.0119
Fungicides Cycloheximide 324 11 483 0.000396 0.0132
Herbicides Metolachlor 377 4 59 0.000648 0.0174
Fungicides Maneb 409 9 360 0.000701 0.0185
Rodenticides Brodifacoum 494 12 644 0.00126 0.0276
POPs Toxaphene 505 4 74 0.00152 0.0318
Metals Arsenic 532 46 4618 0.00169 0.0342
Herbicides Trifluralin 566 2 10 0.00192 0.0364
Insecticides Chlorpyrifos 578 54 5746 0.00212 0.0394
Metals Selenium 680 23 1907 0.00366 0.058
Insecticides Dimethoate 685 3 47 0.0038 0.0584
Herbicides Pendimethalin 756 2 16 0.00498 0.0705
Fungicides Benomyl 775 3 53 0.00533 0.0741
POPs Chlordane 895 5 190 0.00897 0.108
Metals Barium 1122 3 78 0.0153 0.147
Insecticides Cyfluthrin 1179 3 86 0.0199 0.167
Fungicides Ferbam 1197 1 3 0.0199 0.167
Metals Cobalt 1345 6 335 0.025 0.195
Metals Chromium 2293 4 241 0.078 0.367
Rodenticides Cholecalciferol 2419 4 255 0.0914 0.407
PAHs Dibenzo(a, h)pyrene NA NA NA 
sapply(split(pos, pos$Analysis), function(d) c(
    tested = sum(!is.na(d$Rank)),
    adjusted_below_0.05 = sum(d$PValueAdjusted < 0.05, na.rm = TRUE),
    unadjusted_below_0.05 = sum(d$PValue < 0.05, na.rm = TRUE)))
#>                       all_human main min10
#> tested                       34   34    33
#> adjusted_below_0.05          28   23    27
#> unadjusted_below_0.05        34   32    31

Verification and discovery

The two readings answer different questions, and each has a cost.

  • Verification (Yang et al., CTD Batch Query, CTDquerier) tests chemicals chosen in advance. Few tests mean a mild correction and more power for the chemicals one already suspects; the choice itself, though, is a selection the correction does not see, and chemicals outside the list cannot appear.
  • Discovery (ctdR) tests every chemical. The correction is harsher, so some suspected chemicals lose significance (11 of the authors’ 34 tested in the main analysis), but nothing is excluded in advance and the multiplicity is accounted for in full.

ctdR does not replace verification; it comes before it, producing the hypotheses that a targeted study then tests. Two cautions when reading a discovery ranking: very small sets reach significance with two or three genes (read k and N beside the p-value), and very large sets (ethanol, caffeine) reach it by size, because they are the most studied chemicals.

What differs from the original analysis

Difference Effect
CTD release: May 2023 against the authors’ early 2023 a few genes per chemical
MeSH descendants counted by the authors, not by ctdR larger sets in the paper for parent terms
Gene symbols: HGNC symbols of today’s org.Hs.eg.db a renamed gene changes name (for example BVES is now POPDC1)
Universe: 19,000 fixed in the paper; ctdR narrows its background different p-values in discovery mode, same sets
130 genes in Table S2, 131 in the text, 142 in the Expected column p-values reproduced with 130

Yang et al. and ctdR, chemical by chemical

Each of the 95 entries of the paper’s tables, with the authors’ call next to ctdR’s call in discovery, both at 0.05, twice: without correction (the paper’s p-value against ctdR’s raw p-value) and with it (the paper’s FDR within each class against ctdR’s Benjamini-Hochberg p-value over every chemical tested). The table is step 7 of inst/scripts/replicate_yang2024.R.

outcome_levels <- c("significant in both", "significant only in Yang et al.",
    "significant only in ctdR", "significant in neither", "not tested by ctdR")
cbind(
    unadjusted = table(factor(yang_vs_ctdR$Outcome_unadjusted, outcome_levels)),
    adjusted = table(factor(yang_vs_ctdR$Outcome_adjusted, outcome_levels)))
#>                                 unadjusted adjusted
#> significant in both                     40       23
#> significant only in Yang et al.          2       11
#> significant only in ctdR                 7        2
#> significant in neither                  42       55
#> not tested by ctdR                       4        4

Without correction the two analyses agree on most entries; the correction moves entries from “both” to “only in Yang et al.” because the correction over every chemical in CTD is harsher than the one within each class of a list chosen in advance.

The entries on which the calls disagree, with or without correction, or that ctdR did not test:

agree <- c("significant in both", "significant in neither")
disagree <- yang_vs_ctdR[!(yang_vs_ctdR$Outcome_unadjusted %in% agree &
    yang_vs_ctdR$Outcome_adjusted %in% agree), ]
disagree <- disagree[, c("Chemical", "k_paper", "N_paper", "PValue_paper",
    "PValue_discovery", "PValueFDR_paper", "PValueAdjusted_discovery",
    "Outcome_unadjusted", "Outcome_adjusted")]
disagree[, 4:7] <- lapply(disagree[, 4:7], fmt)
knitr::kable(disagree, row.names = FALSE, col.names = c("Chemical", "k", "N",
    "p paper", "p ctdR", "FDR paper", "Adjusted p ctdR", "Unadjusted",
    "Adjusted"))
Chemical k N p paper p ctdR FDR paper Adjusted p ctdR Unadjusted Adjusted
Cholecalciferol 32 2830 0.002 0.0914 0.005 0.407 significant only in Yang et al. significant only in Yang et al.
Selenium 24 1978 0.004 0.00366 0.008 0.058 significant in both significant only in Yang et al.
Cobalt 48 4849 0.003 0.025 0.008 0.195 significant in both significant only in Yang et al.
Pendimethalin 2 21 0.007 0.00498 0.01 0.0705 significant in both significant only in Yang et al.
Chromium 26 2352 0.008 0.078 0.017 0.367 significant only in Yang et al. significant only in Yang et al.
Dimethoate 3 62 0.008 0.0038 0.018 0.0584 significant in both significant only in Yang et al.
Barium 3 81 0.017 0.0153 0.027 0.147 significant in both significant only in Yang et al.
Benomyl 3 62 0.008 0.00533 0.029 0.0741 significant in both significant only in Yang et al.
Chlordane 5 215 0.015 0.00897 0.037 0.108 significant in both significant only in Yang et al.
Dibenzo(a, h)pyrene 1 2 0.007 0.037 not tested by ctdR not tested by ctdR
Cyfluthrin 3 88 0.021 0.0199 0.038 0.167 significant in both significant only in Yang et al.
Ferbam 1 3 0.014 0.0199 0.038 0.167 significant in both significant only in Yang et al.
Cadmium 45 5266 0.05 0.000643 0.063 0.0175 significant only in ctdR significant only in ctdR
Triadimefon 7 492 0.051 0.0205 0.072 0.172 significant only in ctdR significant in neither
Hexachlorobenzene 4 224 0.066 0.0458 0.081 0.268 significant only in ctdR significant in neither
Polychlorinated biphenyls 72 9137 0.057 0.0339 0.098 0.227 significant only in ctdR significant in neither
Hexachlorobenzene 4 224 0.066 0.0458 0.099 0.268 significant only in ctdR significant in neither
Pyrene 5 300 0.054 0.00916 0.107 0.11 significant only in ctdR significant in neither
Methyl parathion 5 101 0.8 0.000164 0.997 0.007 significant only in ctdR significant only in ctdR
Metam 0 105 not tested by ctdR not tested by ctdR
Imazethapyr 0 1 not tested by ctdR not tested by ctdR
Diphacinone 0 1 not tested by ctdR not tested by ctdR
  • Significant only in Yang et al., without correction: chromium and cholecalciferol are parent terms whose sets the authors enlarged with their MeSH descendants (2,352 and 2,830 genes in the paper, 258 and 264 for the exact ChemicalID in ctdR), so the two analyses test different sets.
  • Significant only in ctdR, without correction: the paper’s p-values are just above 0.05 (0.050 to 0.066), except methyl parathion, the transcription error discussed above.
  • Not tested by ctdR: Metam has no CTD identifier; the other three have sets of one or two genes in the paper and do not reach ctdR’s minimum of two genes in discovery.

The full table, inst/extdata/yang2024/yang2024_comparison.tsv, also holds the authors’ test applied to ctdR’s sets for every entry.

Session information

sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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.32.so;  LAPACK version 3.12.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] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [3] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [5] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [7] BiocGenerics_0.59.12        generics_0.1.4             
#>  [9] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [11] ctdR_0.99.12                BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] DBI_1.3.0                   httr2_1.3.0                
#>   [3] GSEABase_1.75.0             rlang_1.3.0                
#>   [5] magrittr_2.0.5              otel_0.2.0                 
#>   [7] compiler_4.6.2              RSQLite_3.53.3             
#>   [9] DelayedMatrixStats_1.35.0   png_0.1-9                  
#>  [11] vctrs_0.7.3                 memuse_4.2-3               
#>  [13] SpatialExperiment_1.23.0    pkgconfig_2.0.3            
#>  [15] crayon_1.5.3                fastmap_1.2.0              
#>  [17] magick_2.9.1                dbplyr_2.6.0               
#>  [19] XVector_0.53.0              labeling_0.4.3             
#>  [21] rmarkdown_2.32              tzdb_0.5.0                 
#>  [23] graph_1.91.0                purrr_1.2.2                
#>  [25] bit_4.6.0                   xfun_0.62                  
#>  [27] cachem_1.1.0                jsonlite_2.0.0             
#>  [29] blob_1.3.0                  rhdf5filters_1.25.5        
#>  [31] DelayedArray_0.39.8         Rhdf5lib_2.1.0             
#>  [33] BiocParallel_1.47.0         parallel_4.6.2             
#>  [35] R6_2.6.1                    bslib_0.12.0               
#>  [37] RColorBrewer_1.1-3          limma_3.99.4               
#>  [39] jquerylib_0.1.4             Rcpp_1.1.2                 
#>  [41] knitr_1.52                  GSVA_2.7.28                
#>  [43] readr_2.2.0                 BiocBaseUtils_1.15.1       
#>  [45] Matrix_1.7-6                tidyselect_1.2.1           
#>  [47] abind_1.4-8                 yaml_2.3.12                
#>  [49] codetools_0.2-20            curl_8.1.0                 
#>  [51] lattice_0.23-1              tibble_3.3.1               
#>  [53] withr_3.0.3                 KEGGREST_1.53.6            
#>  [55] S7_0.2.2                    evaluate_1.0.5             
#>  [57] BiocFileCache_3.3.1         Biostrings_2.81.9          
#>  [59] pillar_1.11.1               BiocManager_1.30.27        
#>  [61] filelock_1.0.3              vroom_1.7.1                
#>  [63] hms_1.1.4                   ggplot2_4.0.3              
#>  [65] sparseMatrixStats_1.25.0    scales_1.4.0               
#>  [67] xtable_1.8-8                glue_1.8.1                 
#>  [69] maketools_1.3.2             tools_4.6.2                
#>  [71] BiocIO_1.23.3               sys_3.4.3                  
#>  [73] data.table_1.18.6.1         fgsea_1.39.2               
#>  [75] annotate_1.91.0             buildtools_1.0.0           
#>  [77] XML_3.99-0.25               fastmatch_1.1-8            
#>  [79] cowplot_1.2.0               rhdf5_2.57.20              
#>  [81] grid_4.6.2                  SingleCellExperiment_1.35.2
#>  [83] AnnotationDbi_1.75.2        HDF5Array_1.41.3           
#>  [85] cli_3.6.6                   S4Arrays_1.13.2            
#>  [87] dplyr_1.2.1                 gtable_0.3.6               
#>  [89] sass_0.4.10                 digest_0.6.39              
#>  [91] SparseArray_1.13.4          rjson_0.2.23               
#>  [93] org.Hs.eg.db_3.24.4         farver_2.1.2               
#>  [95] memoise_2.0.1               htmltools_0.5.9            
#>  [97] lifecycle_1.0.5             h5mread_1.5.3              
#>  [99] httr_1.4.9                  statmod_1.5.2              
#> [101] bit64_4.8.6