--- title: "GeneSCARAB" author: - name: Marcos Ramos-González affiliation: - &ibvf "Institute for Plant Biochemistry and Photosynthesis, University of Seville - Spanish National Research Council, Spain" - &ccia "Department of Computer Science and Artificial Intelligence, University of Seville, Spain" - name: Emma Serrano-Pérez affiliation: *ibvf - name: Víctor Ramos-González affiliation: *ccia - name: Mercedes García-González affiliation: *ibvf - name: Francisco J. Romero-Campero affiliation: - *ibvf - *ccia email: mramos5@us.es package: GeneSCARAB abstract: > Enrichment analyses are the gold standard for interpreting omics data, allowing for the identification of overrepresented biological processes or sets of elements (genes, metabolites, etc.). However, in chronobiology, a field where the goal is to determine at what time of day each process occurs, these analyses can be misleading due to assumptions of linearity and background size that do not hold true in most studies. GeneSCARAB addresses this by directly analyzing the cohesion and circular uniformity of the phases of the elements in each set, and expands the range of possible analyses by offering multimodal analysis, timing comparisons across various conditions, analysis of the contribution of individual elements, and a variety of publication-ready plots. In terms of input, GeneSCARAB accepts results from popular phase-determination packages such as RAIN, CircaCompare, or Metacycle, facilitating its integration into established pipelines output: BiocStyle::html_document: toc: true toc_depth: 2 toc_float: true vignette: > %\VignetteIndexEntry{GeneSCARAB} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r knit, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction To date, although many recent studies have focused on proposing best practices for circular analysis (Humphreys & Ruxton, 2017; Landler et al., 2018, 2020, 2021), their adoption in omics-based chronobiology has been scarce. Often, the goal of these studies is to determine at what time of day each biological process occurs, analyzing large volumes of data. The simplest analysis typically involves estimating the time of maximum expression (or abundance) in the cycle for each element, referred as phase of the element. This allows for the grouping of each element in bins based on discrete intervals for the different phases. In turn, enrichment analysis using one of the standard bioinformatics ontologies can be applied to each bin, based on software such as Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005). These assume monotonic data and do not account for distances in a circular space, incorrectly measuring those distances and leading to potential errors when the goal is to determine when a process occurs. Similarly, once the bins are created, the total number of elements associated with each is highly variable, so the enrichments may be biased as the background changes in each bin. This is common, for example, in circadian transcriptomic studies, in which there exist high and low transcription time windows, resulting in highly unbalanced bins. GeneSCARAB (Gene Set Circular Analysis and Rhythms Ascertainment in Bioprocesses) is designed for analyzing the cohesion and clustering toward specific points in the cycle of biological sets (genes, metabolites, ...) for chronobiology studies following the proposed best practices in circular statistics. It employs rigorous tests to ensure maximum statistical power in the analyses and functions to accurately plot the results. Also, GeneSCARAB seamlessly integrates into common pipelines for this kind of studies, using as input standard outputs from rhythm detection software, and standard Bioconductor annotation packages to adapt the analysis to any species of interest. ## Installation GeneSCARAB can be installed as follows: ```{r install, echo=FALSE, message=FALSE, eval=FALSE} #This package is currently under review for inclusion in Bioconductor. #You can install the development version from GitHub (ramosgonzmarc/GeneSCARAB) if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes") remotes::install_github("ramosgonzmarc/GeneSCARAB") ``` ## Quick start The first step of this tutorial consist on loading the package and the associated data. The annotation package used in this tutorial can be loaded using the `load_example_annot` function. ```{r setup, echo=FALSE, message=FALSE} library(GeneSCARAB) org_Otaurireduced_eg_db <- load_example_annot() data("circa_table_genescarab") total_phases_table_ld <- data.frame(names = rownames( circa_table_genescarab ), phase = as.numeric( circa_table_genescarab[["ld.peak.time.hours"]] )) total_phases_table_sd <- data.frame(names = rownames( circa_table_genescarab ), phase = as.numeric( circa_table_genescarab[["sd.peak.time.hours"]] )) ``` The basic pipeline implemented in GeneSCARAB consists of taking an annotation package, a table with the estimated phases for each gene, and a set of GO terms of interest, and the software will determine the parameters of the circular distributions for each process, as well as significance measures to identify rhythmic processes with respect to one or more diel time points. Circular uniformity testing is performed via Rayleigh (Rayleigh, 1880), Kuiper (Kuiper, 1960), Hermans-Rasson (HERMANS & RASSON, 1985; Landler et al., 2019b) and Rao's spacing (Rao, 1976; Landler et al., 2019a) tests. For this, we will use data from the paper "Multiomics integration unveils photoperiodic plasticity in the molecular rhythms of marine phytoplankton" ([doi.org/10.1093/plcell/koaf033] (https://doi.org/10.1093/plcell/koaf033)). We will first generate a subset of GOs using the annotation package. ```{r gos_sel, echo=FALSE, message=FALSE} functional_data <- AnnotationDbi::select(org_Otaurireduced_eg_db, keys = AnnotationDbi::keys(org_Otaurireduced_eg_db, keytype = "GID"), columns = c("GID", "GO") ) complete_gos <- unique(functional_data$GO) complete_gos <- complete_gos[!is.na(complete_gos)] set.seed(2345) subset_gos <- complete_gos[sample( 1:length(complete_gos), 50, replace = FALSE )] ``` Then we generate a GO terms list with genes per GO term, including ancestors of selected GO terms. ```{r gos_list, echo=FALSE, message=FALSE} go.list.test <- create_gene_list_go( go_vector = subset_gos, org.package = org_Otaurireduced_eg_db, go_column = "GO", id_column = "GID" ) ``` We can use this list along with the total phase table to create a list of gene phases per GO for the short-day (SD) condition, and clean it by removing GOs with no associated rhythmic genes. ```{r phases_list, echo=FALSE, message=FALSE} phases.list.sd <- gene_list_to_phases( go.list.test, total_phases_table_sd ) phases.list.sd.clean <- phases.list.sd[which( sapply(phases.list.sd, nrow) != 0 )] ``` Based on this list, the `complete_circular_table` function allows the user to calculate a table showing the circular distribution of the phases of the elements in each set (genes by GO term in this case) in terms of quartiles, mean, variance, and magnitude (rho), as well as the simple and BH-adjusted p-values of the Rayleigh, Kuiper, Hermans-Rasson, and Rao’s spacing tests against circular uniformity. This function accepts the list of input phases as key arguments, as well as three types of additional arguments. The function automatically runs only the Rayleigh and Kuiper tests to save computation time in situations where HR and Rao are not necessary. Thus, `force.hr.th` indicates the p-value threshold reached by Rayleigh or Kuiper below which the user does not consider it necessary to run HR (because uniformity would have been already rejected), while `hr.on.large.sets.th` refers to the maximum size of the sets for which we want to apply HR. GeneSCARAB runs the tests sequentially, so these same arguments are used to control the execution of the Rao test, the last of the three. HR and Rao require permutations, so `iter.hr`and `iter.rao` controls the number of iterations for each run. ```{r circa_table, warning=FALSE} go.circa.table.sd <- complete_circular_table( phase.list = phases.list.sd.clean, hr.on.large.sets.th = 200, force.hr.th = 0.04, rao.on.large.sets.th = 200, force.rao.th = 0.04, iter.rao = 999, iter.hr = 999 ) ``` Some software returns discrete phases rather than continuous ones. In this case, the `complete_circular_table_grouped` function works the same as `complete_circular_table`, with the difference that it incorporates tiebreaking in the tests used as described in Landler et al., 2020) and also allows the user to control the execution of a permutation-based Kuiper test. Now we can explore the results. It is recommended to filter by n \> 2 (or more), as HR may give significant results for gene sets with a small number of associated rhythmic genes (even 1). Additionally, since the tests performed when using HR and Rao are typically only for negative results in Rayleigh and Kuiper, the expected distribution of p-values will favor high p-values, which can lead to an increase in false negatives when BH-adjusting them. Therefore, it may be more informative to filter by p-value in these cases. In any case, this decision is left to the user, and tests can be forced to be performed always, by setting p-value thresholds to 0 and size thresholds to a large number. ```{r rhythmic_sd} go.circa.table.sd <- subset(go.circa.table.sd, n > 2) hr_boolean_sd <- go.circa.table.sd$hr_p_value < 0.05 hr_boolean_sd[is.na(hr_boolean_sd)] <- FALSE rao_boolean_sd <- go.circa.table.sd$rao_p_value < 0.05 rao_boolean_sd[is.na(rao_boolean_sd)] <- FALSE sd_rhythmic <- rownames(go.circa.table.sd)[ go.circa.table.sd$rayleigh_p_value_adj < 0.05 | go.circa.table.sd$kuiper_p_value_adj < 0.05 | hr_boolean_sd | rao_boolean_sd ] sd_rhythmic <- unique(sd_rhythmic[!is.na(sd_rhythmic)]) length(sd_rhythmic) ``` We can plot these results with 3 different visualizations. Let's plot a sample of 8 biological processes. First a circular boxplot. This function exports a png containing the plot. MetBrewer's palettes are used for coloring. ```{r sd_boxplot, fig.cap = "Circular boxplot showing the phase distribution of genes associated to some example GO terms."} ex_plots <- c( "GO:0016684", "GO:0006220", "GO:0012507", "GO:0042592", "GO:0016780", "GO:0012507", "GO:0008169" ) plot_table_7 <- go.circa.table.sd[ex_plots, ] circular_boxplot(plot_table_7, color.palette = "Tam" ) ``` GeneSCARAB also offers a circular dotplot visualization. ```{r sd_dotplot, fig.cap = "Circular dotplot showing the phase distribution of genes associated to some example GO terms."} plot_phase_list_7 <- phases.list.sd.clean[ex_plots] circular_dotplot(plot_phase_list_7, color.palette = "Tam" ) ``` And also circular histograms. `nbins` control the number of breaks. ```{r sd_hist, fig.cap = "Circular histogram showing the phase distribution of genes associated to some example GO terms."} circular_histogram(plot_phase_list_7, color.palette = "Tam", nbins = 48) ``` Sets that reject uniformity based on HR or Rao but not by Rayleigh or Kuiper are probably multimodal. These kind of GOs may be further analyzed using `multimodal_analysis` function. Taking advantage of `r CRANpkg("OptCirClust")` (Debnath & Song, 2021), it takes a phase table from an individual set, which has already been determined to be non-uniform, and allows the user to identify its modality, identifying processes with 24, 12, and 8-hour periods. It also clusters the elements (genes) into each of the corresponding modes and recalculates their associated circular statistics tables. ```{r multi_analysis} multi_table <- go.circa.table.sd[which( go.circa.table.sd$hr_p_value < 0.05 ), ] multi_go <- rownames(multi_table)[1] phase_table_multi_go <- phases.list.sd.clean[[multi_go]] multimodal_list <- multimodal_analysis(phase_table_multi_go) ``` This GO term presents significant bimodality, with circular medians in ZT4 and ZT15 and high rhos We can also use GeneSCARAB to plot these clusters. When creating the list, remember that it must be a named list. ```{r multi_dotplot, fig.cap = "Circular dotplot showing two distinct clusters of genes belonging to a bimodal gene set."} circular_dotplot( phases.list = list( cluster1 = multimodal_list$cluster_1, cluster2 = multimodal_list$cluster_2 ), color.palette = "Tam", tr.height = 0.15 ) ``` ## Other ways of creating sets #### KEGG Just as in the previous case for GO terms, GeneSCARAB allows you to automatically generate lists of genes associated with specific KEGG pathways based on an annotation package and a set of KO terms to include (or calculate it for all those appearing in the annotation package) using the `create_gene_list_kegg` function instead of `create_gene_list_go`. In addition, it allows you to filter pathways based on whether they belong to bacteria, archaea, plants, animals, fungi, or protists. ```{r kegg_set, echo=FALSE, message=FALSE} ko.list.test <- create_gene_list_kegg( ko_vector = "K10666", org.package = org_Otaurireduced_eg_db, ko_column = "KO", id_column = "GID", ko_prefix = "map", species = "plants" ) ``` #### Custom gene sets GeneSCARAB also allows the user to use custom gene lists, for example, based on associations found in other databases or previous analyses. All you need to do is define a list structure in which each element corresponds to a vector of genes (or elements). ```{r custom_set} setA <- unique(functional_data$GID)[1:50] setB <- unique(functional_data$GID)[100:120] setC <- unique(functional_data$GID)[400:450] custom.list.test <- list(setA = setA, setB = setB, setC = setC) ``` And, as in the previous cases, the gene_list_to_phases function creates the phase table for the rhythmic genes contained in each of the custom gene sets. ```{r phases_custom} phases.list.custom <- gene_list_to_phases( custom.list.test, total_phases_table_ld ) ``` ## Comparison of phase distributions of sets under two conditions Different biological processes or gene sets may exhibit significantly different circular distributions under different experimental conditions. GeneSCARAB allows for a direct statistical comparison using the `test_two_dist` function. This is useful, for example, for determining advancement or delay in the diel expression of different processes for different conditions or genotypes. For example, we will compare data collected under short-day (SD) and long-day (LD) conditions. For this, we first generate the circular statistics table for LD condition. ```{r rhythmic_ld, warning=FALSE} phases.list.ld <- gene_list_to_phases( go.list.test, total_phases_table_ld ) phases.list.ld.clean <- phases.list.ld[which( sapply(phases.list.ld, nrow) != 0 )] go.circa.table.ld <- complete_circular_table( phase.list = phases.list.ld.clean, hr.on.large.sets.th = 200, force.hr.th = 0.04, rao.on.large.sets.th = 200, force.rao.th = 0.04, iter.rao = 999, iter.hr = 999 ) go.circa.table.ld <- subset(go.circa.table.ld, n > 2) hr_boolean_ld <- go.circa.table.ld$hr_p_value < 0.05 hr_boolean_ld[is.na(hr_boolean_ld)] <- FALSE rao_boolean_ld <- go.circa.table.ld$rao_p_value < 0.05 rao_boolean_ld[is.na(rao_boolean_ld)] <- FALSE ld_rhythmic <- rownames(go.circa.table.ld)[ go.circa.table.ld$rayleigh_p_value_adj < 0.05 | go.circa.table.ld$kuiper_p_value_adj < 0.05 | hr_boolean_ld | rao_boolean_ld ] ld_rhythmic <- unique(ld_rhythmic[!is.na(ld_rhythmic)]) length(ld_rhythmic) ``` After this, we compare gene phase distributions under both conditions for the sets identified as non-uniform under both LD and SD. `test_two_dist` renders a table with the comparison set-by-set, offering a p-value and adjusted p-value for each comparison, based on the result of MANOVA over the sines and cosines of the phases (Landler et al., 2021, 2022). ```{r dist_compare} rhythmic_gos <- intersect(ld_rhythmic, sd_rhythmic) diff_distributed_gos <- test_two_dist( phases.list.ld.clean[rhythmic_gos], phases.list.sd.clean[rhythmic_gos] ) sum(diff_distributed_gos$p_value_adj < 0.05) ``` And we can plot one of these. ```{r compare_boxplot, fig.cap = "Circular boxplot showing differences in phase distribution of genes associated to GO:0006334 gene set due to photoperiod of entrainment."} go_circa_plot_diff <- rbind( go.circa.table.ld["GO:0006334", ], go.circa.table.sd["GO:0006334", ] ) rownames(go_circa_plot_diff) <- c("LD", "SD") circular_boxplot(go_circa_plot_diff, color.palette = "Austria") ``` ## Identification of sets that do not follow the experimental phase distribution GeneSCARAB also compares a set’s phases against the overall phase distribution of the entire study. Biologically, it relates to sets of elements that display specific and tightly scheduled temporal placement, suggesting that precise temporal control is likely functionally important for these sets. This comparison is performed by the `test_against_gen_dist` function, which also relies on MANOVA test on the sine and cosine values of the phases of the elements in the set and of the entire set phases (Landler et al., 2021, 2022). ```{r diff_dist} not_experimentally_distributed_gos <- test_against_gen_dist( phases.list.ld.clean[rhythmic_gos], total_phases_table_ld ) sum(not_experimentally_distributed_gos$gen_dist_p_value_adj < 0.05) ``` ## Contribution of each individual element to the complete phase distribution of the set To determine the relative importance of each element to its set's statistics, GeneSCARAB also measures the individual contribution of each element to the mean and circular dispersion of the complete group using leave-one-out and ranks them. ```{r gene_contribution} processes_for_gene_contribution <- rhythmic_gos[100:105] gene_contributions <- gene_contribution_to_set( circa_result = go.circa.table.ld[processes_for_gene_contribution, ], phase_list = phases.list.ld.clean ) ``` Finally, we can identify the genes that have the greatest impact on the mean and variance and visualize these genes with respect to the complete distribution. ```{r ranked_contribution, fig.cap = "Circular dotplot showing the gene whose removal exhibits the highest impact on the mean of the gene set."} ranked_mean <- sapply(gene_contributions, function(x) rownames(x)[1]) ranked_mean ranked_rho <- sapply(gene_contributions, function(x) { names( which.min(x[, "Rho_rank"]) ) }) ranked_rho gene_of_interest_mean <- subset( phases.list.ld.clean$`GO:0016773`, names == ranked_mean[1] ) circular_dotplot( phases.list = list( ranked_gene = gene_of_interest_mean, complete_go = phases.list.ld.clean$`GO:0016773` ), color.palette = "Austria", tr.height = 0.15 ) ``` ## References Debnath, T., & Song, M. (2021). Fast Optimal Circular Clustering and Applications on Round Genomes. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 18(6), 2061–2071. HERMANS, M., & RASSON, J. P. (1985). A new Sobolev test for uniformity on the circle. Biometrika, 72(3), 698–702. Kuiper, N. H. (1960). Tests concerning random points on a circle. Nederl. Akad. Wetensch. Proc. Ser. A, 63(1), 38–47 Landler, L., Ruxton, G. D., & Malkemper, E. P. (2018). Circular data in biology: advice for effectively implementing statistical procedures. Behavioral Ecology and Sociobiology, 72(8), 128. Landler, L., Ruxton, G. D., & Malkemper, E. P. (2019a). Circular statistics meets practical limitations: a simulation-based Rao’s spacing test for non-continuous data. Movement Ecology, 7(1), 15. Landler, L., Ruxton, G. D., & Malkemper, E. P. (2019b). The Hermans–Rasson test as a powerful alternative to the Rayleigh test for circular statistics in biology. BMC Ecology, 19(1), 30. Landler, L., Ruxton, G. D., & Malkemper, E. P. (2020). Grouped circular data in biology: advice for effectively implementing statistical procedures. Behavioral Ecology and Sociobiology, 74(8), 100. Landler, L., Ruxton, G. D., & Malkemper, E. P. (2021). Advice on comparing two independent samples of circular data in biology. Scientific Reports, 11(1), 20337. Landler, L., Ruxton, G. D., & Malkemper, E. P. (2022). The multivariate analysis of variance as a powerful approach for circular data. Movement Ecology, 10(1), 21. Rao, J. S. (1976). Some tests based on arc-lengths for the circle. Sankhyā: The Indian Journal of Statistics, Series B, 329–338. Rayleigh, Lord. (1880). XII. On the resultant of a large number of vibrations of the same pitch and of arbitrary phase. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 10(60), 73–78. ```{r session_info} sessionInfo() ```