Enrichment analysis for 6 Traits, 5 tissues, eQTL + sQTL + stQTL – compute ukbb LD, eqtl, sqtl from predictdb; stQTL from Munro et al
XSun
2024-10-23
Last updated: 2024-10-23
Checks: 6 1
Knit directory: multigroup_ctwas_analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown is untracked by Git. To know which version of the R
Markdown file created these results, you’ll want to first commit it to
the Git repo. If you’re still working on the analysis, you can ignore
this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
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(20231112) 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 2db1bdc. 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: .Rhistory
Untracked files:
Untracked: analysis/multi_group_6traits_15weights_ess_enrichment_genesymbol.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.
There are no past versions. Publish this analysis with
wflow_publish() to start tracking its development.
library(tidyr)
library(dplyr)
library(VennDiagram)
library(ggplot2)
traits <- c("LDL-ukb-d-30780_irnt","SBP-ukb-a-360","WBC-ieu-b-30","aFib-ebi-a-GCST006414","SCZ-ieu-b-5102","IBD-ebi-a-GCST004131")
dbs <- c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023")
trait<- "LDL-ukb-d-30780_irnt"
db <- "GO_Biological_Process_2023"Methods
We do enrichment analysis for the genes with PIP > 0.8 here: https://sq-96.github.io/multigroup_ctwas_analysis/multi_group_6traits_15weights_ess.html
The gene set membership was downloaded here: https://maayanlab.cloud/Enrichr/#libraries
Background genes
For Fractional model and Fisher exact test, we selected 2 kind of backgroud genes
- All genes used in ctwas
- All genes in the selected geneset database.
For enrichR, the background genes are not modifiable. The background genes are all genes in the selected geneset database
Remove the redundant
Using EnrichR package.
This package was used in our earlier ctwas paper.
- It takes a list of genes(genes with PIP > 0.8) as input and returns the enriched GO terms with adjusted p-values.
Fractional model
The model is:
glm(PIP ~ gene set membership, family = quasibinomial('logit')).
We do this regression for one gene set at a time.
The PIP vector contains:
- all genes within the credible set: we use their actual PIPs
- genes without the credible set & PIP < 0.1: we set the PIPs
as
0.5*min(gene pip within credible set)
The 2 different baselines:
- All genes from ctwas. Here,
genes without the credible set & PIP < 0.1includes only the genes used in ctwas. - All genes from the geneset database. Here,
genes without the credible set & PIP < 0.1includes the union of all genes from the GO terms in the geneset database.
Fisher exact test
We assign 1 to the genes with PIP > 0.5/0.8 & in cs and 0 for others. We name this vector as binarized_PIP. We test the association between the binarized_PIP and geneset_membership.
The testing matrix is:
| geneset_membership | 0 | 1 |
|---|---|---|
| binarized_pip 0 | a | b |
| binarized_pip 1 | c | d |
Where:
ais the count wherebinarized_pip = 0andgeneset_membership = 0.bis the count wherebinarized_pip = 0andgeneset_membership = 1.cis the count wherebinarized_pip = 1andgeneset_membership = 0.dis the count wherebinarized_pip = 1andgeneset_membership = 1.
The 2 different baselines:
- All genes from ctwas. Here,
geneset_membershipmatrix includes only the genes used in ctwas. - All genes from the geneset database. Here,
geneset_membershipmatrix includes the union of all genes from the GO terms in the geneset database.
Summary for the number of Go terms with p-values < 0.001
load("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_summary_for_all_redundant_p.rdata")
DT::datatable(summary,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Number of enriched GO terms under different settings'),options = list(pageLength = 20) )
Comparing the p-values from Enrichr and Fisher exact test – baseline genes are all genes from gene sets
p_enrichr <- c()
p_fet <- c()
#compare_diff <- c()
for (trait in traits) {
for (db in dbs) {
file_enrichr <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_redundant_",trait,"_",db,".rdata")
if(file.exists(file_enrichr)) {
load(file_enrichr)
load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fisher_blgeneset_pip08_",trait,"_",db,".rdata"))
merged <- merge(db_enrichment, summary, by.x = "Term", by.y = "GO")
p_enrichr <- c(p_enrichr, merged$P.value)
p_fet <- c(p_fet, merged$pvalue)
#compare_diff <- rbind(compare_diff, merged)
}
}
}
p_enrichr <- as.numeric(p_enrichr)
p_fet <- as.numeric(p_fet)
df <- data.frame(p_enrichr = p_enrichr, p_fet = p_fet)
# Fit a linear model to calculate the slope
fit <- lm(p_fet ~ p_enrichr)
slope <- coef(fit)[2]
intercept <- coef(fit)[1]
ggplot(df, aes(x = p_enrichr, y = p_fet)) +
geom_point() + # Add points
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # y = x line
geom_smooth(method = "lm", se = FALSE, color = "blue") + # Best-fit line
annotate("text", x = max(p_enrichr) * 0.8, y = max(p_fet) * 0.9,
label = paste0("Slope: ", round(slope, 3)),
color = "blue") + # Slope text
annotate("text", x = max(p_enrichr) * 0.8, y = max(p_enrichr) * 0.8,
label = "y = x", color = "red", size = 5, hjust = -0.1, vjust = -0.5) + # y = x text near the line
ggtitle("Comparison of p-values between enrichr and FET, baseline -- all genes from gene sets") + # Add title
xlab("Enrichr p-values") + # x-axis label
ylab("FET p-values") + # y-axis label
theme_minimal()
Comparing the GO terms reported by FET and Fractional model – p-values < 0.001, baseline genes are genes used in ctwas, redundant terms NOT removed
pval <- 0.001
all_fractional <- c()
all_fet <- c()
for (trait in traits) {
for (db in dbs) {
file_fet <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fisher_blctwas_pip08_",trait,"_",db,".rdata")
load(file_fet)
all_fet <- rbind(all_fet,summary[as.numeric(summary$pvalue) < pval,])
file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_",trait,"_",db,".rdata")
load(file_fractional)
summary$trait <- trait
summary$db <- db
all_fractional <- rbind(all_fractional,summary[as.numeric(summary$pvalue) < pval,])
}
}
all_fractional <- all_fractional[complete.cases(all_fractional),]
all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fet$id <- paste0(all_fet$trait,"-",all_fet$db,"-",all_fet$GO)
venn.plot <- draw.pairwise.venn(
area1 = nrow(all_fractional), # Size of Group A
area2 = nrow(all_fet), # Size of Group B
cross.area = sum(all_fractional$id %in% all_fet$id), # Overlap between Group A and Group B
category = c("Fractional", "FET"), # Labels for the groups
fill = c("red", "blue"), # Colors for the groups
lty = "blank", # Line type for the circles
cex = 2, # Font size for the numbers
cat.cex = 2 # Font size for the labels
)
all_fractional <- all_fractional[,c("trait","db","GO","pvalue","fdr","id")]
DT::datatable(all_fractional[!all_fractional$id %in% all_fet$id,],caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) ) DT::datatable(all_fet[!all_fet$id %in% all_fractional$id,],caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )Comparing the GO terms reported by FET and Fractional model – p-values < 0.001, baseline genes are all genes from gene sets, redundant terms NOT removed
pval <- 0.001
all_fractional <- c()
all_fet <- c()
for (trait in traits) {
for (db in dbs) {
file_fet <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fisher_blgeneset_pip08_",trait,"_",db,".rdata")
load(file_fet)
all_fet <- rbind(all_fet,summary[as.numeric(summary$pvalue) < pval,])
file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_blgeneset_",trait,"_",db,".rdata")
load(file_fractional)
summary$trait <- trait
summary$db <- db
all_fractional <- rbind(all_fractional,summary[as.numeric(summary$pvalue) < pval,])
}
}
all_fractional <- all_fractional[complete.cases(all_fractional),]
all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fet$id <- paste0(all_fet$trait,"-",all_fet$db,"-",all_fet$GO)
venn.plot <- draw.pairwise.venn(
area1 = nrow(all_fractional), # Size of Group A
area2 = nrow(all_fet), # Size of Group B
cross.area = sum(all_fractional$id %in% all_fet$id), # Overlap between Group A and Group B
category = c("Fractional", "FET"), # Labels for the groups
fill = c("red", "blue"), # Colors for the groups
lty = "blank", # Line type for the circles
cex = 2, # Font size for the numbers
cat.cex = 2 # Font size for the labels
)
all_fractional <- all_fractional[,c("trait","db","GO","pvalue","fdr","id")]
DT::datatable(all_fractional[!all_fractional$id %in% all_fet$id,],caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) ) DT::datatable(all_fet[!all_fet$id %in% all_fractional$id,],caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )
sessionInfo()R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /software/openblas-0.3.13-el7-x86_64/lib/libopenblas_haswellp-r0.3.13.so
locale:
[1] C
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggplot2_3.5.1 VennDiagram_1.7.3 futile.logger_1.4.3
[4] dplyr_1.1.4 tidyr_1.3.0
loaded via a namespace (and not attached):
[1] tidyselect_1.2.0 xfun_0.41 bslib_0.3.1
[4] purrr_1.0.2 lattice_0.20-45 splines_4.2.0
[7] colorspace_2.0-3 vctrs_0.6.5 generics_0.1.2
[10] htmltools_0.5.2 yaml_2.3.5 mgcv_1.8-40
[13] utf8_1.2.2 rlang_1.1.2 jquerylib_0.1.4
[16] later_1.3.0 pillar_1.9.0 glue_1.6.2
[19] withr_2.5.0 lambda.r_1.2.4 lifecycle_1.0.4
[22] stringr_1.5.1 munsell_0.5.0 gtable_0.3.0
[25] workflowr_1.7.0 htmlwidgets_1.5.4 evaluate_0.15
[28] labeling_0.4.2 knitr_1.39 fastmap_1.1.0
[31] httpuv_1.6.5 crosstalk_1.2.0 fansi_1.0.3
[34] highr_0.9 Rcpp_1.0.12 promises_1.2.0.1
[37] scales_1.3.0 DT_0.22 formatR_1.12
[40] jsonlite_1.8.0 farver_2.1.0 fs_1.5.2
[43] digest_0.6.29 stringi_1.7.6 rprojroot_2.0.3
[46] cli_3.6.1 tools_4.2.0 magrittr_2.0.3
[49] sass_0.4.1 tibble_3.2.1 futile.options_1.0.1
[52] pkgconfig_2.0.3 Matrix_1.5-3 rmarkdown_2.25
[55] rstudioapi_0.13 R6_2.5.1 nlme_3.1-157
[58] git2r_0.30.1 compiler_4.2.0