In neurogenomics/MultiEWCE: Multi-Scale Target Explorer

library(MSTExplorer)

Introduction

The MSTExplorer package is an extension of the EWCE package. It is designed to run expression weighted celltype enrichment (EWCE) on multiple gene lists in parallel. The results are then stored both as separate .rds files, one for each individual EWCE analysis, as well as a in a single dataframe containing all the results.

This package is useful in cases where you have a large number of related, but separate, gene lists. In this vignette we will use an example from the Human Phenotype Ontology (HPO). The HPO contains over 9000 clinically relevant phenotypes annotated with lists of genes that have been found to be associated with the particular phenotype.

Loading Phenotype Associated Gene Lists from the HPO

The MSTExplorer package requires the gene data to be in a particular format. It must be a data.frame that includes one column of gene list names, and another column of genes. For example:

| hpo_name | Gene | | ---------------- | ------ | | "Abnormal heart" | gene X | | "Abnormal heart" | gene Y | | "Poor vision" | gene Z | | "Poor vision" | gene Y | | "Poor vision" | gene W | | "Short stature" | gene V | etc...

Now we will get a dataset like this from the HPO.

gene_data <- HPOExplorer::load_phenotype_to_genes()
knitr::kable(head(gene_data))

In this example our gene list names column is called Phenotype and our column of genes is called Gene. However, different column names can be specified to the MSTExplorer package.

Setting up input arguments for the gen_results function

# Loading CTD file 
ctd <- load_example_ctd()
list_names <- unique(gene_data$hpo_id)[seq(10)] 
reps <- 10 # in practice would use more reps
cores <- 1 # in practice would use more cores
save_dir <- file.path(tempdir())
save_dir_tmp <- file.path(save_dir,"results")

ctd

The ctd (cell type data) file contains the single cell RNA sequence data that is required for EWCE. for further information about generating a ctd please see the EWCE documentation. In this example we will use a CTD of human gene expression data, generated from the Descartes Human Cell Atlas. Replace this with your own CTD file.

gene_data

Gene data is the dataframe containing gene list names and genes, in this case we have already loaded it and assigned it to the variable gene_data.

list_names

This is a character vector containing all the gene list names. This can be obtained from your gene_data as follows. To save time in this example analysis we will only use the first 10 gene lists ([1:10])

bg

This is a character vector of genes to be used as the background genes. See EWCE package docs for more details on background genes.

Column names

list_name_column is the name of the column in gene_data that contains the gene list names and gene_column contains the genes.

Processing results arguments

The save_dir argument is the path to the directory where the individual EWCE results will be saved.

The force_new argument can be set to TRUE or FALSE and states if you want to redo and overwrite analysis of gene lists that have already been saved to the save_dir. Setting this to FALSE is useful in cases where you stopped an analysis midway and would like to carry on from where you left off.

Number of cores

The cores argument is the number of cores you would like to run in parallel. This is dependent on what is available to you on your computer. In this case we will just run it on one core, no parallelism.

EWCE::bootstrap_enrichment_test arguments

The gen_results function calls the EWCE::bootstrap_enrichment_test function. Here we set the input parameters related to this.

reps is the number of bootstrap reps to run, for this tutorial we will only do 10 to save time, but typically you would want to do closer to 100,000.

Run analysis

Now we have set up all our desired inputs, we can run the analysis.

out <- MSTExplorer::gen_results(ctd = ctd,
                              gene_data = gene_data,
                              list_names = list_names, 
                              list_name_column = "hpo_id",
                              reps = reps,
                              cores = cores,
                              save_dir = save_dir, 
                              force_new = TRUE,
                              save_dir_tmp = save_dir_tmp)
results <- out$results

Visualise the results

Just as an example, we will create a plot showing the number of significant enrichments per phenotype in the all_results data.frame. We will use q <= 0.05 as the significance threshold.

library(ggplot2)
library(data.table)
#### Aggregate results ####
n_signif <- results[q<=0.05 & !is.na(q),
                        list(sig_enrichments = .N,
                             mean_effect=mean(effect)), 
                        by="hpo_id"]
#### Plot #### 
plot1 <-  ggplot(n_signif, aes(x = stringr::str_wrap(hpo_id,width = 10),
                               y = sig_enrichments,
                               fill = mean_effect)) +
    geom_col() +
    labs(x="Phenotype",y="Enrichments (n)") + 
    theme_bw() 

methods::show(plot1)

Other functions in MSTExplorer

merge_results

If you have a results directory of individual EWCE results but do not have the merged dataframe of all results, you can call the merge_results function manually. The save_dir argument is the path to your results directory and the list_name_column argument is the name of the column containing gene list names. In this case we used "Phenotype" as this column name when we generated the results.

all_results_2 <- MSTExplorer::merge_results(save_dir = save_dir_tmp)

get_gene_list

This function gets a character vector of genes associated with a particular gene list name.

phenotypes <- c("Scoliosis")
gene_set <- HPOExplorer::get_gene_lists(phenotypes = phenotypes,
                                        phenotype_to_genes = gene_data)

cat(paste(length(unique(gene_set$gene_symbol)),
          "genes associated with",shQuote(phenotypes),":",
          paste(unique(gene_set$gene_symbol)[seq(5)],collapse = ", ")))

get_unfinished_list_names

This function is used to find which gene lists you have not yet analysed

all_phenotypes <- unique(gene_data$hpo_id)
unfinished <- MSTExplorer::get_unfinished_list_names(list_names = all_phenotypes,
                                                   save_dir_tmp = save_dir_tmp)

methods::show(paste(length(unfinished),"/",length(all_phenotypes),
                    "gene lists not yet analysed"))

Run disease-level enrichment tests

So far, we've iterated over gene list grouped by phenotypes. But we can also do this at the level of diseases (which are composed of combinations of phenotypes).

gene_data <- HPOExplorer::load_phenotype_to_genes("genes_to_phenotype.txt")
#### Filter only to those with >=4 genes ####
gene_counts <- gene_data[,list(genes=length(unique(gene_symbol))),
                         by="disease_id"][genes>=4]
list_names <- unique(gene_counts$disease_id)[seq(5)]

out <- MSTExplorer::gen_results(ctd = ctd,
                              gene_data = gene_data,
                              list_name_column = "disease_id",
                              list_names = list_names,
                              annotLevel = 1,
                              force_new = TRUE,
                              reps = 10)
results <- out$results

Full analysis

Run the following code the replicate the main analysis in the study described here.

gene_data <- HPOExplorer::load_phenotype_to_genes()
gene_data[,n_gene:=length(unique(gene_symbol)),by="hpo_id"]
gene_data <- gene_data[n_gene>=4,]
ctd <- load_example_ctd("ctd_DescartesHuman.rds")

all_results <- MSTExplorer::gen_results(ctd = ctd,
                                      list_name_column = "hpo_id", 
                                      gene_data = gene_data,  
                                      annotLevel = 2,
                                      reps = 100000,
                                      cores = 10)

Session info

utils::sessionInfo()

---

neurogenomics/MultiEWCE documentation built on Sept. 28, 2024, 2:27 a.m.