Enrichment analysis for 6 Traits, 5 tissues, eQTL + sQTL + stQTL – compute ukbb LD, eqtl, sqtl from predictdb; stQTL from Munro et al

workflowr

• Summary
• Checks
• Past versions

Last updated: 2024-11-18

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.

---

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. 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.

Environment: empty

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.

Seed: set.seed(20231112)

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.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 51219ee

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 51219ee. 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
    Ignored:    analysis/figure/

Unstaged changes:
    Modified:   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.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/multi_group_6traits_15weights_ess_enrichment_genesymbol.Rmd) and HTML (docs/multi_group_6traits_15weights_ess_enrichment_genesymbol.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	51219ee	XSun	2024-11-15	update
html	51219ee	XSun	2024-11-15	update
Rmd	0e6914f	XSun	2024-11-15	update
html	0e6914f	XSun	2024-11-15	update
Rmd	ba41054	XSun	2024-10-31	update
html	ba41054	XSun	2024-10-31	update
Rmd	0fa71fe	XSun	2024-10-29	update
html	0fa71fe	XSun	2024-10-29	update
Rmd	ba30ddd	XSun	2024-10-24	update
html	ba30ddd	XSun	2024-10-24	update
Rmd	2446887	XSun	2024-10-23	update
html	2446887	XSun	2024-10-23	update

---

 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"

pval_threshold <- 0.001

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

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

Model

The model is: glm(PIP ~ gene set membership, family = quasibinomial('log10it')). 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:

1. All genes from ctwas. Here, genes without the credible set & PIP < 0.1 includes only the genes used in ctwas.
2. All genes from the geneset database. Here, genes without the credible set & PIP < 0.1 includes the union of all genes from the GO terms in the geneset database.

p-value calibration

We used permutation testing to assess the significance of associations between combined_pip and GO terms. Permutation testing creates a “null” distribution by shuffling data and recalculating p-values.

1. Initial logistic Regression: For each GO term, a logistic regression is performed using glm(PIP ~ gene set membership, family = quasibinomial('log10it')). The observed p-value (pval_origin) measures the association strength in the actual data.

2. Null Distributions via Permutation:

   • Then we simulate random associations by shuffling the GO term membership (x) multiple times and recalculating the logistic regression p-value for each shuffle.
   • These shuffled p-values represent what would be expected if there were no real association and create a “null” distribution.

3. Calibrated p-values with Increasing Permutations:

   • The initial permutation test uses 1,000 shuffles (n_permutations = 1000) to compute a calibrated p-value (pval_calibrated), comparing the observed p-value against the null distribution.
   • Further Calibration: If the calibrated p-value is significant ( pval_calibrated < 0.05), additional rounds of permutations (100,000 shuffles) are conducted to improve accuracy for small p-values. Each round refines the calibrated p-value by expanding the null distribution, enhancing the robustness of significance estimates.

4. Final Significance: The final p-value reflects the proportion of permutation-derived p-values that are more extreme than the observed one. If only a small number of permutation p-values are smaller, the association is considered statistically significant.

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:

• a is the count where binarized_pip = 0 and geneset_membership = 0.
• b is the count where binarized_pip = 0 and geneset_membership = 1.
• c is the count where binarized_pip = 1 and geneset_membership = 0.
• d is the count where binarized_pip = 1 and geneset_membership = 1.

The 2 different baselines:

1. All genes from ctwas. Here, geneset_membership matrix includes only the genes used in ctwas.
2. All genes from the geneset database. Here, geneset_membership matrix includes the union of all genes from the GO terms in the geneset database.

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)

 lgp_enrichr <- log10(p_enrichr)
 lgp_fet <- log10(p_fet)
 
 df <- data.frame(lgp_enrichr = lgp_enrichr, lgp_fet = lgp_fet)

 # Fit a linear model to calculate the slope
 # fit <- lm(lgp_fet ~ lgp_enrichr)
 # slope <- coef(fit)[2]
 # intercept <- coef(fit)[1]

 ggplot(df, aes(x = lgp_enrichr, y = lgp_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(lgp_enrichr) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_enrichr) * 0.8, y = max(lgp_enrichr) * 0.8, 
            label = "y = x", color = "red", size = 5,hjust = 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 log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-2-1.png

Version	Author	Date
0fa71fe	XSun	2024-10-29
ba30ddd	XSun	2024-10-24

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")

#summary_show <- summary[grep(pattern = "ctwasgene",summary$method),] 

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 GO terms reported by FET and Fractional model (pvalues calibrated) – p-values < 0.001, baseline genes are genes used in ctwas, redundant terms NOT removed

aFib-ebi-a-GCST006414

Comparing Calibrated fractional model and FET

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "aFib-ebi-a-GCST006414" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-4-1.png

Version	Author	Date
0fa71fe	XSun	2024-10-29
2446887	XSun	2024-10-23

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-4-4.png

Version	Author	Date
ba41054	XSun	2024-10-31
0fa71fe	XSun	2024-10-29

Examining ctwas PIP values for genes within GO terms associated with high FET p-values.

GO_Cellular_Component_2023 - Cortical Endoplasmic Reticulum (GO:0032541)

db <- "GO_Cellular_Component_2023"
go <- "Cortical Endoplasmic Reticulum (GO:0032541)"

load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_combinedpip.rdata"))
load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_finemap_annot.rdata"))

genes_notincs <- subset(finemap_gene_res, !gene_name %in% combined_pip_by_type$gene_name)[, "gene_name"]
genes_notincs <- data.frame(gene_name = unique(genes_notincs), combined_pip = 0.5 * min(combined_pip_by_type$combined_pip))
genes_all <- rbind(combined_pip_by_type[, c("gene_name", "combined_pip")], genes_notincs)

load(paste0("/project/xinhe/xsun/data/Go/membership_", db, ".rdata"))
membership_ordered <- membership[match(genes_all$gene_name, rownames(membership)), ]

x <- membership_ordered[,colnames(membership_ordered) == go]
gene_ingo <- names(which(x == 1))

print("ctwas PIP for the genes in GO term")

[1] "ctwas PIP for the genes in GO term"

combined_pip_by_type[combined_pip_by_type$gene_name %in% gene_ingo,]

   gene_name combined_pip  eQTL_pip sQTL_pip stQTL_pip
74     STIM1    0.7376334 0.7376334       NA        NA

GO_Biological_Process_2023 - Cytochrome Complex Assembly (GO:0017004)

db <- "GO_Biological_Process_2023"
go <- "Cytochrome Complex Assembly (GO:0017004)"

load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_combinedpip.rdata"))
load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_finemap_annot.rdata"))

genes_notincs <- subset(finemap_gene_res, !gene_name %in% combined_pip_by_type$gene_name)[, "gene_name"]
genes_notincs <- data.frame(gene_name = unique(genes_notincs), combined_pip = 0.5 * min(combined_pip_by_type$combined_pip))
genes_all <- rbind(combined_pip_by_type[, c("gene_name", "combined_pip")], genes_notincs)

load(paste0("/project/xinhe/xsun/data/Go/membership_", db, ".rdata"))
membership_ordered <- membership[match(genes_all$gene_name, rownames(membership)), ]

x <- membership_ordered[,colnames(membership_ordered) == go]
gene_ingo <- names(which(x == 1))

print("ctwas PIP for the genes in GO term")

[1] "ctwas PIP for the genes in GO term"

combined_pip_by_type[combined_pip_by_type$gene_name %in% gene_ingo,]

   gene_name combined_pip eQTL_pip  sQTL_pip stQTL_pip
77     COX18    0.7264782       NA 0.7264782        NA

GO_Molecular_Function_2023 - Calcium Channel Regulator Activity (GO:0005246)

db <- "GO_Molecular_Function_2023"
go <- "Calcium Channel Regulator Activity (GO:0005246)"

load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_combinedpip.rdata"))
load(paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/", trait, "_finemap_annot.rdata"))

genes_notincs <- subset(finemap_gene_res, !gene_name %in% combined_pip_by_type$gene_name)[, "gene_name"]
genes_notincs <- data.frame(gene_name = unique(genes_notincs), combined_pip = 0.5 * min(combined_pip_by_type$combined_pip))
genes_all <- rbind(combined_pip_by_type[, c("gene_name", "combined_pip")], genes_notincs)

load(paste0("/project/xinhe/xsun/data/Go/membership_", db, ".rdata"))
membership_ordered <- membership[match(genes_all$gene_name, rownames(membership)), ]

x <- membership_ordered[,colnames(membership_ordered) == go]
gene_ingo <- names(which(x == 1))

print("ctwas PIP for the genes in GO term")

[1] "ctwas PIP for the genes in GO term"

combined_pip_by_type[combined_pip_by_type$gene_name %in% gene_ingo,]

    gene_name combined_pip    eQTL_pip  sQTL_pip stQTL_pip
52       SGK1  0.944520556 0.944520556        NA        NA
74      STIM1  0.737633431 0.737633431        NA        NA
79      YWHAE  0.713523514          NA 0.7135235        NA
439     NRXN1  0.001469662 0.001469662        NA        NA

LDL-ukb-d-30780_irnt

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "LDL-ukb-d-30780_irnt" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-8-1.png

Version	Author	Date
ba41054	XSun	2024-10-31

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-8-4.png

Version	Author	Date
ba41054	XSun	2024-10-31

IBD-ebi-a-GCST004131

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "IBD-ebi-a-GCST004131" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-9-1.png

Version	Author	Date
ba41054	XSun	2024-10-31

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-9-4.png

Version	Author	Date
ba41054	XSun	2024-10-31

SBP-ukb-a-360

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "SBP-ukb-a-360" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-10-1.png

Version	Author	Date
ba41054	XSun	2024-10-31

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-10-4.png

Version	Author	Date
ba41054	XSun	2024-10-31

SCZ-ieu-b-5102

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "SCZ-ieu-b-5102" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-11-1.png

Version	Author	Date
ba41054	XSun	2024-10-31

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-11-4.png

Version	Author	Date
ba41054	XSun	2024-10-31

WBC-ieu-b-30

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "WBC-ieu-b-30" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-12-1.png

Version	Author	Date
0e6914f	XSun	2024-11-15

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-12-4.png

Version	Author	Date
0e6914f	XSun	2024-11-15

Comparing the GO terms reported by FET and Fractional model (pvalues calibrated) – p-values < 0.001, baseline genes are all genes in geneset data base, redundant terms NOT removed

aFib-ebi-a-GCST006414

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "aFib-ebi-a-GCST006414" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-13-1.png

Version	Author	Date
0e6914f	XSun	2024-11-15

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-13-4.png

Version	Author	Date
0e6914f	XSun	2024-11-15

LDL-ukb-d-30780_irnt

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "LDL-ukb-d-30780_irnt"

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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)

    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-14-1.png

Version	Author	Date
51219ee	XSun	2024-11-15
0e6914f	XSun	2024-11-15

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"),
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)

df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9,
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8,
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-14-4.png

Version	Author	Date
51219ee	XSun	2024-11-15
0e6914f	XSun	2024-11-15

IBD-ebi-a-GCST004131

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "IBD-ebi-a-GCST004131"

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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)

    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-15-1.png

Version	Author	Date
51219ee	XSun	2024-11-15
0e6914f	XSun	2024-11-15

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"),
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)

df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9,
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8,
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-15-4.png

Version	Author	Date
51219ee	XSun	2024-11-15
0e6914f	XSun	2024-11-15

SBP-ukb-a-360

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "SBP-ukb-a-360" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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
   )

Past versions of unnamed-chunk-16-1.png

Version	Author	Date
51219ee	XSun	2024-11-15

all_fractional <- all_fractional[,c("trait","db","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Past versions of unnamed-chunk-16-4.png

Version	Author	Date
51219ee	XSun	2024-11-15

SCZ-ieu-b-5102

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "SCZ-ieu-b-5102" 
 
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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)
 
df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9, 
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8, 
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

WBC-ieu-b-30

all_fractional <- c()
all_fet <- c()
supporting_genes <- c()

trait <- "WBC-ieu-b-30"

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)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)

    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

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)
#supporting_genes$id <- paste0(supporting_genes$trait,"-",supporting_genes$db,"-",supporting_genes$GO)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fet_pass <- all_fet[as.numeric(all_fet$pvalue) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fet_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fet_pass$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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fet, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO.x","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pvalue","fdr","Overlap","Genes")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fet","fdr_fet","Overlap_fet","Overlapped_Genes_fet")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"),
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fet) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fet<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fet) < pval_threshold,]
unique_fet <- unique_fet[complete.cases(unique_fet),]
DT::datatable(unique_fet,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for FET'),options = list(pageLength = 10) )

p_fractional_calibrated <- as.numeric(unique_fractional$pvalue_calibrated_fractional)
p_fet <- as.numeric(unique_fractional$pvalue_fet)

lgp_fractional_calibrated <- log10(p_fractional_calibrated)
lgp_fet <- log10(p_fet)

df <- data.frame(lgp_fractional_calibrated = lgp_fractional_calibrated, lgp_fet = lgp_fet)

 ggplot(df, aes(x = lgp_fractional_calibrated, y = lgp_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(lgp_fractional_calibrated) * 0.8, y = max(lgp_fet) * 0.9,
   #          label = paste0("Slope: ", round(slope, 3)),
   #          color = "blue") +  # Slope text
   annotate("text", x = max(lgp_fractional_calibrated) * 0.8, y = max(lgp_fractional_calibrated) * 0.8,
            label = "y = x", color = "red", size = 5, hjust = 1, vjust = -0.5) +  # y = x text near the line
   ggtitle("Comparison of p-values between \n fractional_calibrated and FET \n for unique fractional GO terms") +  # Add title
   xlab("fractional_calibrated log10(p)") +  # x-axis label
   ylab("FET log10(p)") +  # y-axis label
   theme_minimal()

Comparing the GO terms reported by FGSEA and Fractional model (pvalues calibrated) – p-values < 0.001, baseline genes are all genes in geneset data base, redundant terms NOT removed

aFib-ebi-a-GCST006414

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "aFib-ebi-a-GCST006414" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

LDL-ukb-d-30780_irnt

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "LDL-ukb-d-30780_irnt" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

IBD-ebi-a-GCST004131

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "IBD-ebi-a-GCST004131" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

SBP-ukb-a-360

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "SBP-ukb-a-360" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

SCZ-ieu-b-5102

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "SCZ-ieu-b-5102" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

WBC-ieu-b-30

all_fractional <- c()
all_fgsea <- c()
supporting_genes <- c()

trait <- "WBC-ieu-b-30" 
 
for (db in dbs) {

    file_fgsea <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fgsea_blgeneset_",trait,"_",db,".rdata")
    load(file_fgsea)
    all_fgsea <- rbind(all_fgsea,fgseaRes)

    file_fractional <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_calibrated_blgeneset_",trait,"_",db,".rdata")
    load(file_fractional)
    summary$trait <- trait
    summary$db <- db
    all_fractional <- rbind(all_fractional,summary)
    
    file_supporting <- paste0("/project/xinhe/xsun/multi_group_ctwas/11.multi_group_1008/postprocess/enrichment_fractional_supporting_genes_",trait,"_",db,".rdata")
    load(file_supporting)
    supporting_genes <- rbind(supporting_genes,gene_supp_db)

}

all_fractional$id <- paste0(all_fractional$trait,"-",all_fractional$db,"-",all_fractional$GO)
all_fgsea$id <- paste0(all_fgsea$trait,"-",all_fgsea$db,"-",all_fgsea$pathway)

fractional_pass <- all_fractional[as.numeric(all_fractional$pvalue_calibrated) < pval_threshold, ]
fractional_pass <-  fractional_pass[complete.cases( fractional_pass$fdr_origin),]
fgsea_pass <- all_fgsea[as.numeric(all_fgsea$pval) < pval_threshold, ]

venn.plot <- draw.pairwise.venn(
     area1 = nrow(fractional_pass),          # Size of Group A
     area2 = nrow(fgsea_pass),          # Size of Group B
     cross.area = sum(fractional_pass$id %in%  fgsea_pass$id),     # Overlap between Group A and Group B
     category = c("Fractional", "fgsea"),  # 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_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","id")]

merged_two <- merge(all_fractional,all_fgsea, by = "id")
merged_two <- merged_two[,c("trait.x","db.x","GO","pvalue_origin","fdr_origin","pvalue_calibrated","fdr_calibrated","pval","padj","leadingEdge")]
colnames(merged_two) <- c("trait","db","GO","pvalue_origin_fractional","fdr_origin_fractional","pvalue_calibrated_fractional","fdr_calibrated_fractional","pvalue_fgsea","fdr_fgsea","leadingEdge_fgsea")
merged_two <- merge(merged_two, supporting_genes[, c("trait", "db", "GO", "supporting_genes")],
                     by = c("trait", "db", "GO"), 
                     all.x = TRUE)
colnames(merged_two)[which(colnames(merged_two) == "supporting_genes")] <- "supporting_genes(pip>0.5)_fractional"

unique_fractional <- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) < pval_threshold & as.numeric(merged_two$pvalue_fgsea) > pval_threshold,]
unique_fractional <- unique_fractional[complete.cases(unique_fractional$trait),]
unique_fractional <- unique_fractional[order(as.numeric(unique_fractional$pvalue_calibrated_fractional)),]

DT::datatable(unique_fractional,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fractional model'),options = list(pageLength = 10) )

unique_fgsea<- merged_two[as.numeric(merged_two$pvalue_calibrated_fractional) > pval_threshold & as.numeric(merged_two$pvalue_fgsea) < pval_threshold,]
unique_fgsea <- unique_fgsea[complete.cases(unique_fgsea$trait),]
DT::datatable(unique_fgsea,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Unique GO terms for fgsea'),options = list(pageLength = 10) )

Session information

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          colorspace_2.0-3     vctrs_0.6.5         
 [7] generics_0.1.2       htmltools_0.5.2      yaml_2.3.5          
[10] utf8_1.2.2           rlang_1.1.2          jquerylib_0.1.4     
[13] later_1.3.0          pillar_1.9.0         glue_1.6.2          
[16] withr_2.5.0          lambda.r_1.2.4       lifecycle_1.0.4     
[19] stringr_1.5.1        munsell_0.5.0        gtable_0.3.0        
[22] workflowr_1.7.0      htmlwidgets_1.5.4    evaluate_0.15       
[25] labeling_0.4.2       knitr_1.39           fastmap_1.1.0       
[28] crosstalk_1.2.0      httpuv_1.6.5         fansi_1.0.3         
[31] highr_0.9            Rcpp_1.0.12          promises_1.2.0.1    
[34] scales_1.3.0         DT_0.22              formatR_1.12        
[37] jsonlite_1.8.0       farver_2.1.0         fs_1.5.2            
[40] digest_0.6.29        stringi_1.7.6        rprojroot_2.0.3     
[43] cli_3.6.1            tools_4.2.0          magrittr_2.0.3      
[46] sass_0.4.1           tibble_3.2.1         futile.options_1.0.1
[49] whisker_0.4          pkgconfig_2.0.3      rmarkdown_2.25      
[52] rstudioapi_0.13      R6_2.5.1             git2r_0.30.1        
[55] compiler_4.2.0