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.
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:
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.
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
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
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 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:
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
The two readings answer different questions, and each has a cost.
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.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.
| 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 |
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 |
ChemicalID in ctdR), so the two analyses test different sets.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.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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: America/New_York
#> 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.1 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 tinytex_0.61
#> [25] purrr_1.2.2 bit_4.6.0
#> [27] xfun_0.62 cachem_1.1.0
#> [29] jsonlite_2.0.0 blob_1.3.0
#> [31] rhdf5filters_1.25.5 DelayedArray_0.39.8
#> [33] Rhdf5lib_2.1.0 BiocParallel_1.47.0
#> [35] parallel_4.6.1 R6_2.6.1
#> [37] bslib_0.12.0 RColorBrewer_1.1-3
#> [39] limma_3.99.4 jquerylib_0.1.4
#> [41] Rcpp_1.1.2 bookdown_0.48
#> [43] knitr_1.52 GSVA_2.7.28
#> [45] readr_2.2.0 BiocBaseUtils_1.15.1
#> [47] Matrix_1.7-6 tidyselect_1.2.1
#> [49] dichromat_2.0-1 abind_1.4-8
#> [51] yaml_2.3.12 codetools_0.2-20
#> [53] curl_8.1.0 lattice_0.23-1
#> [55] tibble_3.3.1 withr_3.0.3
#> [57] KEGGREST_1.53.6 S7_0.2.2
#> [59] evaluate_1.0.5 BiocFileCache_3.3.2
#> [61] Biostrings_2.81.9 pillar_1.11.1
#> [63] BiocManager_1.30.27 filelock_1.0.3
#> [65] vroom_1.7.1 hms_1.1.4
#> [67] ggplot2_4.0.3 sparseMatrixStats_1.25.0
#> [69] scales_1.4.0 xtable_1.8-8
#> [71] glue_1.8.1 tools_4.6.1
#> [73] BiocIO_1.23.3 data.table_1.18.6.1
#> [75] fgsea_1.39.2 annotate_1.91.0
#> [77] XML_3.99-0.25 fastmatch_1.1-8
#> [79] cowplot_1.2.0 rhdf5_2.57.20
#> [81] grid_4.6.1 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