Nonlinear Dynamic eQTL analysis on iPSC
Ziang Zhang
2025-02-21
Last updated: 2025-11-09
Checks: 7 0
Knit directory: fashr-paper/
This reproducible R Markdown analysis was created with workflowr (version 1.7.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20251109) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version bb17d9a. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: code/.DS_Store
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: analysis/toy_example.rmd
Untracked: code/00_eQTLs.R
Untracked: code/01_fash.R
Untracked: code/01_fash_uncorrected.R
Untracked: code/02_dyn.R
Untracked: code/03_nonlindyn.R
Untracked: code/04_minlfsr.R
Untracked: code/05_Interact.R
Untracked: code/06_minlfsr_nonlin.R
Untracked: code/07_grid_sensitivity.R
Untracked: code/08_grid_sensitivity.R
Untracked: code/filterVariantPerGene.R
Untracked: data/appendixB/
Untracked: data/dynamic_eQTL_real/
Untracked: output/appendixB/
Untracked: output/dynamic_eQTL_real/
Unstaged changes:
Modified: analysis/about.Rmd
Modified: analysis/index.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown
(analysis/nonlinear_dynamic_eQTL_real.rmd) and HTML
(docs/nonlinear_dynamic_eQTL_real.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | bb17d9a | Ziang Zhang | 2025-11-09 | workflowr::wflow_publish("analysis/nonlinear_dynamic_eQTL_real.rmd") |
knitr::opts_chunk$set(fig.width = 8, fig.height = 6)
library(fashr)
result_dir <- paste0(getwd(), "/output/dynamic_eQTL_real")
data_dir <- paste0(getwd(), "/data/dynamic_eQTL_real")
code_dir <- paste0(getwd(), "/code/dynamic_eQTL_real")
log_prec <- seq(0,10, by = 0.2)
fine_grid <- sort(c(0, exp(-0.5*log_prec)))Obtain the effect size of eQTLs
We use the processed (expression & genotype) data of Strober et.al, 2019 to perform the eQTL analysis.
For the association testing, we use a linear regression model for each gene-variant pair at each time point. Following the practice in Strober et.al, we adjust for the first five PCs.
The code to perform this step can be found in the script
dynamic_eQTL_real/00_eQTLs.R from the code directory.
After this step, we have the effect size of eQTLs for each gene-variant pair at each time point, as well as its standard error.
Fitting FASH
To fit the FASH model on \(\{\beta_i(t_j), s_{ij}\}_{i\in N,j \in [16]}\), we consider fitting two FASH models:
A FASH model based on first order IWP (testing for dynamic eQTLs: \(H_0: \beta_i(t)=c\)).
A FASH model based on second order IWP (testing for nonlinear-dynamic eQTLs: \(H_0: \beta_i(t)=c_1+c_2t\)).
The code to perform this step can be found in the script
dynamic_eQTL_real/01_fash.R from the code directory.
We will directly load the fitted FASH models from the output directory.
load(paste0(result_dir, "/fash_fit2_all.RData"))We will load the datasets from the fitted FASH object:
datasets <- fash_fit2$fash_data$data_list
for (i in 1:length(datasets)) {
datasets[[i]]$SE <- fash_fit2$fash_data$S[[i]]
}
all_genes <- unique(sapply(strsplit(names(datasets), "_"), "[[", 1))In this analysis, we will focus on the FASH(2) model that assumes a second order IWP and tests for nonlinear dynamic eQTLs.
Let’s take a quick overview of the fitted FASH model:
log_prec <- seq(0,10, by = 0.2)
fine_grid <- sort(c(0, exp(-0.5*log_prec)))
fash_fit2 <- fash(Y = "beta", smooth_var = "time", S = "SE", data_list = datasets,
num_basis = 20, order = 2, betaprec = 0,
pred_step = 1, penalty = 10, grid = fine_grid,
num_cores = num_cores, verbose = TRUE)
save(fash_fit2, file = "./results/fash_fit2_all.RData")fash_fit2Fitted fash Object
-------------------
Number of datasets: 1009173
Likelihood: gaussian
Number of PSD grid values: 52 (initial), 9 (non-trivial)
Order of Integrated Wiener Process (IWP): 2As well as the estimated priors:
fash_fit2$prior_weights psd prior_weight
1 0.000000000 9.872378e-01
2 0.006737947 1.209985e-04
3 0.020241911 5.257484e-03
4 0.022370772 2.211188e-03
5 0.055023220 1.113735e-03
6 0.060810063 3.628552e-03
7 0.246596964 3.827515e-04
8 0.449328964 1.965219e-05
9 1.000000000 2.786579e-05Problem with \(\pi_0\) estimation
The original MLE estimated \(\pi_0\) is 0.9872378. This could be under-estimated due to model-misspecification under the alternative hypothesis. To account for this, we will consider the following approaches:
(i): A conservative estimate of \(\pi_0\) based on the BF procedure:
fash_fit2_update <- BF_update(fash_fit2, plot = FALSE)
fash_fit2_update$prior_weights
save(fash_fit2_update, file = paste0(result_dir, "/fash_fit2_update.RData"))The conservative estimate is 0.9992301, which is much more conservative.
(ii): Instead of looking at the FDR which is based on the estimated \(\pi_0\), we can use the minimum local false sign rate (\(\text{min-lfsr}_i\)) to measure significance: \[ \text{min-lfsr}_i = \min_{t} \left\{ \text{lfsr}(W_i(t)) \right\}, \] where \(W_i(t) = \beta_i(t) - \beta_i(0)\).
Detecting Nonlinear dynamic eQTLs
We will use the updated FASH model (2) to detect nonlinear dynamic eQTLs.
alpha <- 0.05
test2 <- fdr_control(fash_fit2_update, alpha = alpha, plot = F)96 datasets are significant at alpha level 0.05. Total datasets tested: 1009173. fash_highlighted2 <- test2$fdr_results$index[test2$fdr_results$FDR <= alpha]How many pairs are detected?
pairs_highlighted2 <- names(datasets)[fash_highlighted2]
length(pairs_highlighted2)[1] 96length(pairs_highlighted2)/length(datasets)[1] 9.51274e-05How many unique genes are detected?
genes_highlighted2 <- unique(sapply(strsplit(pairs_highlighted2, "_"), "[[", 1))
length(genes_highlighted2)[1] 27length(genes_highlighted2)/length(all_genes)[1] 0.004243948Visualize top-ranked pairs:
Some examples of null pairs:
genes_not_highlighted2 <- setdiff(all_genes, genes_highlighted2)
pairs_not_highlighted2 <- setdiff(names(datasets), pairs_highlighted2)
par(mfrow = c(2,2))
for (i in 10:13) {
selected_gene <- (genes_not_highlighted2)[i] #sample(genes_highlighted2, 1)
pairs_of_selected_gene <- grep(selected_gene, pairs_not_highlighted2, value = T)
selected_indices <- which(names(datasets) %in% pairs_of_selected_gene)
selected_index <- selected_indices[which(fash_fit2_update$lfdr[selected_indices] == max(fash_fit2_update$lfdr[selected_indices]))][1]
fitted_result <- predict(fash_fit2_update,
index = selected_index,
smooth_var = seq(0, 15, by = 0.1))
plot(
datasets[[selected_index]]$x,
datasets[[selected_index]]$y,
pch = 20,
col = "black",
xlab = "Time",
ylab = "Effect Est",
main = paste0(names(datasets)[selected_index]),
ylim = c(
min(datasets[[selected_index]]$y - 2 * datasets[[selected_index]]$SE),
max(datasets[[selected_index]]$y + 2 * datasets[[selected_index]]$SE)
)
)
arrows(
datasets[[selected_index]]$x,
datasets[[selected_index]]$y - 2 * datasets[[selected_index]]$SE,
datasets[[selected_index]]$x,
datasets[[selected_index]]$y + 2 * datasets[[selected_index]]$SE,
length = 0.05,
angle = 90,
code = 3,
col = "black"
)
lines(fitted_result$x,
fitted_result$mean,
col = "red",
lwd = 2)
abline(h = mean(datasets[[selected_index]]$y), col = "blue", lty = 2)
polygon(
c(fitted_result$x, rev(fitted_result$x)),
c(fitted_result$lower, rev(fitted_result$upper)),
col = rgb(1, 0, 0, 0.3),
border = NA
)
}
par(mfrow = c(1,1))Comparing with Strober et.al
We will compare the detected dynamic eQTLs with the results from Strober et.al.
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.1 ✔ stringr 1.5.2
✔ ggplot2 4.0.0 ✔ tibble 3.3.0
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.1.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ tidyr::unite() masks ggVennDiagram::unite()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsLet’s take a look at the overlap between the two methods used in Strober et.al and FASH (order 2):
gene_sets <- list(
"Strober (Nonlinear)" = genes_highlighted_strober_nonlinear,
"Strober (Linear)" = genes_highlighted_strober_linear,
"FASH (2)" = genes_highlighted2
)
ggVennDiagram(gene_sets, label = "both") +
scale_fill_gradient(low="grey90",high = "red") +
theme(legend.position = "right") # Move legend to the right
There is a large number of genes only detected by FASH (order 2). Let’s take a look at the 4 pairs that are least significant from FASH:
Let’s also look at the genes that were missed by FASH, but detected by Strober et.al. In this case, we will pick the most significant pair for each gene in FASH:
Classifying nonlinear dynamic eQTLs
Following the definition in Strober et.al, we will classify the detected eQTLs into different categories:
Early: eQTLs with strongest effect during the first three days: \(\max_{t\leq3} |\beta(t)| - \max_{t> 3} |\beta(t)| > 0\).
Late: eQTLs with strongest effect during the last four days: \(\max_{t\geq 12} |\beta(t)| - \max_{t< 12} |\beta(t)| > 0\).
Middle: eQTLs with strongest effect during days 4-11: \(\max_{4\leq t\leq 11} |\beta(t)| - \max_{t> 11 | t< 4} |\beta(t)| > 0\).
Switch: eQTLs with effect sign switch during the time course: ${(t)+,(t)-}-c $ where \(c\) is a threshold that we set to 0.25 (which means with two alleles, the maximal difference of effect size is at least \(\geq 2\times\min\{\max\beta(t)^+,\max\beta(t)^-\}\times2 \geq 2 \times 0.25 \times 2 = 1\)).
We first take a look at the significant pairs detected by FASH (order 2), and classify them based on the false sign rate (lfsr):
smooth_var_refined = seq(0,15, by = 0.1)
functional_early <- function(x){
max(abs(x[smooth_var_refined <= 3])) - max(abs(x[smooth_var_refined > 3]))
}
testing_early_nonlin_dyn <- testing_functional(functional_early,
lfsr_cal = function(x){mean(x <= 0)},
fash = fash_fit2,
indices = fash_highlighted2,
smooth_var = smooth_var_refined)How many pairs and how many unique genes are classified as early dynamic eQTLs?
load(paste0(result_dir, "/classify_nonlin_dyn_eQTLs_early.RData"))
early_indices <- testing_early_nonlin_dyn$indices[testing_early_nonlin_dyn$cfsr <= alpha]
length(early_indices)[1] 0early_genes <- unique(sapply(strsplit(names(datasets)[early_indices], "_"), "[[", 1))
length(early_genes)[1] 0How many pairs are classified as middle dynamic eQTLs?
functional_middle <- function(x){
max(abs(x[smooth_var_refined <= 11 & smooth_var_refined >= 4])) - max(abs(x[smooth_var_refined > 11]), abs(x[smooth_var_refined < 4]))
}
testing_middle_nonlin_dyn <- testing_functional(functional_middle,
lfsr_cal = function(x){mean(x <= 0)},
fash = fash_fit2,
indices = fash_highlighted2,
num_cores = num_cores,
smooth_var = smooth_var_refined)load(paste0(result_dir, "/classify_nonlin_dyn_eQTLs_middle.RData"))
middle_indices <- testing_middle_nonlin_dyn$indices[testing_middle_nonlin_dyn$cfsr <= alpha]
length(middle_indices)[1] 32middle_genes <- unique(sapply(strsplit(names(datasets)[middle_indices], "_"), "[[", 1))
length(middle_genes)[1] 7Take a look at their results:
par(mfrow = c(2,2))
for (i in 1:4) {
selected_index <- sample(middle_indices, 1)
fitted_result <- predict(fash_fit2,
index = selected_index,
smooth_var = seq(0, 15, by = 0.1))
plot(
datasets[[selected_index]]$x,
datasets[[selected_index]]$y,
pch = 20,
col = "black",
xlab = "Time",
ylab = "Effect Est",
main = paste0(names(datasets)[selected_index]),
ylim = c(
min(datasets[[selected_index]]$y - 2 * datasets[[selected_index]]$SE),
max(datasets[[selected_index]]$y + 2 * datasets[[selected_index]]$SE)
)
)
arrows(
datasets[[selected_index]]$x,
datasets[[selected_index]]$y - 2 * datasets[[selected_index]]$SE,
datasets[[selected_index]]$x,
datasets[[selected_index]]$y + 2 * datasets[[selected_index]]$SE,
length = 0.05,
angle = 90,
code = 3,
col = "black"
)
lines(fitted_result$x,
fitted_result$mean,
col = "red",
lwd = 2)
abline(h = 0, lty = 2, col = "blue")
polygon(
c(fitted_result$x, rev(fitted_result$x)),
c(fitted_result$lower, rev(fitted_result$upper)),
col = rgb(1, 0, 0, 0.3),
border = NA
)
}
par(mfrow = c(1,1))How many pairs are classified as late dynamic eQTLs?
functional_late <- function(x){
max(abs(x[smooth_var_refined >= 12])) - max(abs(x[smooth_var_refined < 12]))
}
testing_late_nonlin_dyn <- testing_functional(functional_late,
lfsr_cal = function(x){mean(x <= 0)},
fash = fash_fit2,
indices = fash_highlighted2,
num_cores = num_cores,
smooth_var = smooth_var_refined)load(paste0(result_dir, "/classify_nonlin_dyn_eQTLs_late.RData"))
late_indices <- testing_late_nonlin_dyn$indices[testing_late_nonlin_dyn$cfsr <= alpha]
length(late_indices)[1] 45late_genes <- unique(sapply(strsplit(names(datasets)[late_indices], "_"), "[[", 1))
length(late_genes)[1] 14Let’s take a look at the top-ranked late dynamic eQTLs:
How many pairs and how many unique genes are classified as switch dynamic eQTLs?
switch_threshold <- 0.25
functional_switch <- function(x){
# compute the radius of x, measured by deviation from 0 from below and from above
x_pos <- x[x > 0]
x_neg <- x[x < 0]
if(length(x_pos) == 0 || length(x_neg) == 0){
return(0)
}
min(max(abs(x_pos)), max(abs(x_neg))) - switch_threshold
}
testing_switch_nonlin_dyn <- testing_functional(functional_switch,
lfsr_cal = function(x){mean(x <= 0)},
fash = fash_fit2,
indices = fash_highlighted2,
num_cores = num_cores,
smooth_var = smooth_var_refined)load(paste0(result_dir, "/classify_nonlin_dyn_eQTLs_switch.RData"))
switch_indices <- testing_switch_nonlin_dyn$indices[testing_switch_nonlin_dyn$cfsr <= alpha]
length(switch_indices)[1] 9switch_genes <- unique(sapply(strsplit(names(datasets)[switch_indices], "_"), "[[", 1))
length(switch_genes)[1] 9Let’s take a look at the top-ranked switch dynamic eQTLs:
Gene Set Enrichment Analysis
library(clusterProfiler)
library(tidyverse)
library(msigdbr)
library(org.Hs.eg.db) # Assuming human genes
library(biomaRt)
library(cowplot)
# Retrieve Hallmark gene sets for Homo sapiens
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>%
dplyr::select(gs_name, entrez_gene)
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
## A function to check gene-enrichment
enrich_set <- function(genes_selected, background_gene, q_val_cutoff = 0.05, pvalueCutoff = 0.05) {
genes_converted <- getBM(
filters = "ensembl_gene_id",
attributes = c("ensembl_gene_id", "entrezgene_id"),
values = genes_selected,
mart = mart
)
# Extract Entrez IDs from the converted data
entrez_gene_list <- genes_converted$entrezgene_id
genes_converted_all <- getBM(
filters = "ensembl_gene_id",
attributes = c("ensembl_gene_id", "entrezgene_id"),
values = background_gene,
mart = mart
)
entrez_universe <- as.character(genes_converted_all$entrezgene_id)
entrez_universe <- entrez_universe[!is.na(entrez_universe)]
# Perform enrichment analysis using Hallmark gene sets
enrich_res <- enricher(pAdjustMethod = "BH",
entrez_gene_list,
TERM2GENE = m_t2g,
qvalueCutoff = q_val_cutoff,
pvalueCutoff = pvalueCutoff,
universe = entrez_universe)
enrich_res
}Among all the genes highlighted by FASH:
result <- enrich_set(genes_selected = genes_highlighted2, background_gene = all_genes)
result@result %>%
filter(pvalue < 0.05) %>%
dplyr::select(GeneRatio, BgRatio, pvalue, qvalue)[1] GeneRatio BgRatio pvalue qvalue
<0 rows> (or 0-length row.names)Among the genes highlighted by FASH that are classified as middle dynamic eQTLs:
result <- enrich_set(genes_selected = middle_genes, background_gene = all_genes)
result@result %>%
filter(pvalue < 0.05) %>%
dplyr::select(GeneRatio, BgRatio, pvalue, qvalue)[1] GeneRatio BgRatio pvalue qvalue
<0 rows> (or 0-length row.names)Among the genes highlighted by FASH that are classified as late dynamic eQTLs:
result <- enrich_set(genes_selected = late_genes, background_gene = all_genes)
result@result %>%
filter(pvalue < 0.05) %>%
dplyr::select(GeneRatio, BgRatio, pvalue, qvalue) GeneRatio BgRatio pvalue qvalue
HALLMARK_MYOGENESIS 2/6 67/1837 0.01788046 0.1129292Among the genes highlighted by FASH that are classified as switch dynamic eQTLs:
result <- enrich_set(genes_selected = switch_genes, background_gene = all_genes)
result@result %>%
filter(pvalue < 0.05) %>%
dplyr::select(GeneRatio, BgRatio, pvalue, qvalue) GeneRatio BgRatio pvalue qvalue
HALLMARK_MYOGENESIS 2/6 67/1837 0.01788046 0.1882153
sessionInfo()R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Chicago
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] cowplot_1.2.0 biomaRt_2.64.0 org.Hs.eg.db_3.21.0
[4] AnnotationDbi_1.70.0 IRanges_2.42.0 S4Vectors_0.46.0
[7] Biobase_2.68.0 BiocGenerics_0.54.1 generics_0.1.4
[10] msigdbr_25.1.1 clusterProfiler_4.16.0 lubridate_1.9.4
[13] forcats_1.0.1 stringr_1.5.2 dplyr_1.1.4
[16] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
[19] tibble_3.3.0 ggplot2_4.0.0 tidyverse_2.0.0
[22] ggVennDiagram_1.5.4 fashr_0.1.30 workflowr_1.7.2
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 rstudioapi_0.17.1 jsonlite_2.0.0
[4] magrittr_2.0.4 ggtangle_0.0.8 farver_2.1.2
[7] rmarkdown_2.30 fs_1.6.6 vctrs_0.6.5
[10] memoise_2.0.1 ggtree_3.16.3 mixsqp_0.3-54
[13] progress_1.2.3 htmltools_0.5.8.1 curl_7.0.0
[16] gridGraphics_0.5-1 sass_0.4.10 bslib_0.9.0
[19] httr2_1.2.1 plyr_1.8.9 cachem_1.1.0
[22] TMB_1.9.18 whisker_0.4.1 igraph_2.2.0
[25] lifecycle_1.0.4 pkgconfig_2.0.3 gson_0.1.0
[28] Matrix_1.7-3 R6_2.6.1 fastmap_1.2.0
[31] GenomeInfoDbData_1.2.14 digest_0.6.37 numDeriv_2016.8-1.1
[34] aplot_0.2.9 enrichplot_1.28.4 colorspace_2.1-2
[37] patchwork_1.3.2 ps_1.9.1 rprojroot_2.1.1
[40] irlba_2.3.5.1 RSQLite_2.4.3 filelock_1.0.3
[43] labeling_0.4.3 timechange_0.3.0 httr_1.4.7
[46] compiler_4.5.1 bit64_4.6.0-1 withr_3.0.2
[49] S7_0.2.0 BiocParallel_1.42.2 DBI_1.2.3
[52] R.utils_2.13.0 rappdirs_0.3.3 tools_4.5.1
[55] ape_5.8-1 httpuv_1.6.16 R.oo_1.27.1
[58] glue_1.8.0 callr_3.7.6 nlme_3.1-168
[61] GOSemSim_2.34.0 promises_1.3.3 grid_4.5.1
[64] getPass_0.2-4 reshape2_1.4.4 fgsea_1.34.2
[67] gtable_0.3.6 tzdb_0.5.0 R.methodsS3_1.8.2
[70] data.table_1.17.8 hms_1.1.3 xml2_1.4.0
[73] XVector_0.48.0 ggrepel_0.9.6 pillar_1.11.1
[76] babelgene_22.9 yulab.utils_0.2.1 later_1.4.4
[79] splines_4.5.1 BiocFileCache_2.16.2 treeio_1.32.0
[82] lattice_0.22-7 bit_4.6.0 tidyselect_1.2.1
[85] GO.db_3.21.0 Biostrings_2.76.0 knitr_1.50
[88] git2r_0.36.2 xfun_0.53 LaplacesDemon_16.1.6
[91] stringi_1.8.7 UCSC.utils_1.4.0 lazyeval_0.2.2
[94] ggfun_0.2.0 yaml_2.3.10 evaluate_1.0.5
[97] codetools_0.2-20 qvalue_2.40.0 ggplotify_0.1.3
[100] cli_3.6.5 processx_3.8.6 jquerylib_0.1.4
[103] dichromat_2.0-0.1 Rcpp_1.1.0 GenomeInfoDb_1.44.3
[106] dbplyr_2.5.1 png_0.1-8 parallel_4.5.1
[109] assertthat_0.2.1 blob_1.2.4 prettyunits_1.2.0
[112] DOSE_4.2.0 tidytree_0.4.6 scales_1.4.0
[115] crayon_1.5.3 rlang_1.1.6 fastmatch_1.1-6
[118] KEGGREST_1.48.1