# knowYourCG

A tool for functional analysis of DNA methylomes

# Quick Start

knowYourCG is a tool for evaluating the enrichment of CpG probes in different methylation feature sets. These features can be categorical (e.g., CpGs located at tissue-specific transcription factors) or continuous (e.g., the local CpG density at a regulatory element). Additionally, the set of CpGs to which the test will be applied can be categorical or continuous as well.

The set of CpGs tested for enrichment is called the query set, and the curated target features are called the database sets. A query set, for example, might be the results of a differential methylation analysis or an epigenome-wide association study. We have curated a variety of database sets that represent different categorical and continuous methylation features such as CpGs associated with chromatin states, technical artifacts, gene association and gene expression correlation, transcription factor binding sites, tissue specific methylation, CpG density, etc.

The following commands prepare the use of knowYourCG. Several database sets are retrieved and caching is performed to enable faster access in future enrichment testing. More information on viewing and accessing available database sets is discussed later on.

    library(knowYourCG)
    library(sesameData)
    sesameDataCache(data_titles=c("genomeInfo.hg38","genomeInfo.mm10",
                      "KYCG.MM285.tissueSignature.20211211",
                      "MM285.tissueSignature","MM285.address",
                      "probeIDSignature","KYCG.MM285.designGroup.20210210",
                      "KYCG.MM285.chromHMM.20210210",
                      "KYCG.MM285.TFBSconsensus.20220116",
                      "KYCG.MM285.HMconsensus.20220116",
                      "KYCG.MM285.chromosome.mm10.20210630"
                      ))

The following example uses a query of CpGs methylated in mouse primordial germ cells (design group PGCMeth). First get the CG list using the following code.

    query <- getDBs("MM285.designGroup")[["PGCMeth"]]
    head(query)

    ## [1] "cg36615889_TC11" "cg36646136_BC21" "cg36647910_BC11" "cg36857173_TC21"
    ## [5] "cg36877289_BC21" "cg36899653_BC21"

Now test the enrichment. By default, KYCG will select all the categorical groups available but we can specify a subset of databases.

    dbs <- c("KYCG.MM285.chromHMM.20210210",
             "KYCG.HM450.TFBSconsensus.20211013",
             "KYCG.MM285.HMconsensus.20220116",
             "KYCG.MM285.tissueSignature.20211211",
             "KYCG.MM285.chromosome.mm10.20210630",
             "KYCG.MM285.designGroup.20210210")
    results_pgc <- testEnrichment(query,databases = dbs,platform="MM285")
    head(results_pgc)

As expected, the PGCMeth group itself appears on the top of the list. But one can also find histone H3K9me3, chromHMM `Het` and transcription factor `Trim28` binding enriched in this CG group.

# Testing Scenarios

There are four testing scenarios depending on the type format of the query set and database sets. They are shown with the respective testing scenario in the table below. `testEnrichment`, `testEnrichmentSEA` are for Fisher’s exact test and Set Enrichment Analysis respectively.

Four knowYourCG Testing Scenarios
| Continuous Database Set | Discrete Database Set
--- | --- | ---
Continuous Query | Correlation-based | Set Enrichment Analysis
Discrete Query | Set Enrichment Analysis | Fisher’s Exact Test

# Enrichment Testing

The main work horse function for testing enrichment of a categorical query against categorical databases is the `testEnrichment` function. This function will perform Fisher’s exact testing of the query against each database set (one-tailed by default, but two-tailed optionally) and reports overlap and enrichment statistics.

> Choice of universe set: Universe set is the set of all probes for a given platform. It can either be passed in as an argument called `universeSet` or the platform name can be passed with argument `platform`. If neither of these are supplied, the universe set will be inferred from the probes in the query.

    library(SummarizedExperiment)

    ## prepare a query
    df <- rowData(sesameDataGet('MM285.tissueSignature'))
    query <- df$Probe_ID[df$branch == "fetal_brain" & df$type == "Hypo"]

    results <- testEnrichment(query, "TFBS", platform="MM285")
    results %>% dplyr::filter(overlap>10) %>% head

    ## prepare another query
    query <- df$Probe_ID[df$branch == "fetal_liver" & df$type == "Hypo"]
    results <- testEnrichment(query, "TFBS", platform="MM285")
    results %>% dplyr::filter(overlap>10) %>%
        dplyr::select(dbname, estimate, test, FDR) %>% head

The output of each test contains multiple variables including: the estimate (fold enrichment), p-value, overlap statistics, type of test, as well as the name of the database set and the database group. By default, the results are sorted by -log10 of the of p-value and the fold enrichment.

The `nQ` and `nD` columns identify the length of the query set and the database set, respectively. Often, it’s important to examine the extent of overlap between the two sets, so that metric is reported as well in the `overlap` column.

# Database Sets

The success of enrichment testing depends on the availability of biologically-relevant databases. To reflect the biological meaning of databases and facilitate selective testing, we have organized our database sets into different groups. Each group contains one or multiple databases. Here is how to find the names of pre-built database groups:

    listDBGroups("MM285")

The `listDBGroups()` function returns a data frame containing information of these databases. The Title column is the accession key one needs for the `testEnrichment` function. With the accessions, one can either directly use them in the `testEnrichment` function or explicitly call the `getDBs()` function to retrieve databases themselves. Caching these databases on the local machine is important, for two reasons: it limits the number of requests sent to the Bioconductor server, and secondly it limits the amount of time the user needs to wait when re-downloading database sets. For this reason, one should run `sesameDataCache()` before loading in any database sets. This will take some time to download all of the database sets but this only needs to be done once per installation. During the analysis the database sets can be identified using these accessions. knowYourCG also does some guessing when a unique substring is given. For example, the string “MM285.designGroup” retrieves the “KYCG.MM285.designGroup.20210210” database. Let’s look at the database group which we had used as the query (query and database are reciprocal) in our first example:

    dbs <- getDBs("MM285.design")

    ## Selected the following database groups:

    ## 1. KYCG.MM285.designGroup.20210210

In total, 32 datasets have been loaded for this group. We can get the “PGCMeth” as an element of the list:

    str(dbs[["PGCMeth"]])

    ##  chr [1:474] "cg36615889_TC11" "cg36646136_BC21" "cg36647910_BC11" ...
    ##  - attr(*, "group")= chr "KYCG.MM285.designGroup.20210210"
    ##  - attr(*, "dbname")= chr "PGCMeth"

On subsequent runs of the `getDBs()` function, the database loading can be faster thanks to the sesameData in-memory caching, if the corresponding database has been loaded.

# Query Set(s)

A query set represents probes of interest. It may either be in the form of a character vector where the values correspond to probe IDs or a named numeric vector where the names correspond to probe IDs. The query and database definition is rather arbitrary. One can regard a database as a query and turn a query into a database, like in our first example. In real world scenario, query can come from differential methylation testing, unsupervised clustering, correlation with a phenotypic trait, and many others. For example, we could consider CpGs that show tissue-specific methylation as the query. We are getting the B-cell-specific hypomethylation.

    df <- rowData(sesameDataGet('MM285.tissueSignature'))
    query <- df$Probe_ID[df$branch == "B_cell"]
    head(query)

    ## [1] "cg32668003_TC11" "cg45118317_TC11" "cg37563895_TC11" "cg46105105_BC11"
    ## [5] "cg47206675_TC21" "cg38855216_TC21"

This query set represents hypomethylated probes in Mouse B-cells from the MM285 platform. This specific query set has 168 probes.

# Gene Enrichment

A special case of set enrichment is to test whether CpGs are associated with specific genes. Automating the enrichment test process only works when the number of database sets is small. This is important when targeting all genes as there are tens of thousands of genes on each platform. By testing only those genes that overlap with the query set, we can greatly reduce the number of tests. For this reason, the gene enrichment analysis is a special case of these enrichment tests. We can perform this analysis using the `buildGeneDBs()` function.

    query <- names(sesameData_getProbesByGene("Dnmt3a", "MM285"))
    results <- testEnrichment(query,
        buildGeneDBs(query, max_distance=100000, platform="MM285"),
        platform="MM285")
    main_stats <- c("dbname","estimate","gene_name","FDR", "nQ", "nD", "overlap")
    results[,main_stats]

As expected, we recover our targeted gene (Dnmt3a).

Gene enrichment testing can easily be included with default or user specified database sets by setting include_genes=TRUE:

    query <- names(sesameData_getProbesByGene("Dnmt3a", "MM285"))
    dbs <- c("KYCG.MM285.chromHMM.20210210","KYCG.HM450.TFBSconsensus.20211013",
             "KYCG.MM285.chromosome.mm10.20210630")
    results <- testEnrichment(query,databases=dbs,
                              platform="MM285",include_genes=TRUE)
    main_stats <- c("dbname","estimate","gene_name","FDR", "nQ", "nD", "overlap")
    results[,main_stats] %>%
        head()

# GO/Pathway Enrichment

One can get all the genes associated with a probe set and test the Gene Ontology of the probe-associated genes using the `testGO()` function, which internally utilizes g:Profiler2 for the enrichment analysis:

    library(gprofiler2)
    df <- rowData(sesameDataGet('MM285.tissueSignature'))
    query <- df$Probe_ID[df$branch == "fetal_liver" & df$type == "Hypo"]
    res <- testGO(query, platform="MM285",organism = "mmusculus")
    head(res$result)

# Genomic Proximity Testing

Sometimes it may be of interest whether a query set of probes share close genomic proximity. Co-localization may suggest co-regulation or co-occupancy in the same regulatory element. KYCG can test for genomic proximity using the `testProbeProximity()`function. Poisson statistics for the expected # of co-localized hits from the given query size (lambda) and the actual co-localized CpG pairs along with the p value are returned:

    df <- rowData(sesameDataGet('MM285.tissueSignature'))
    probes <- df$Probe_ID[df$branch == "fetal_liver" & df$type == "Hypo"]
    res <- testProbeProximity(probeIDs=probes)
    head(res)

    ## $Stats
    ##    nQ Hits Lambda        P.val
    ## 1 194    4   0.08 2.554721e-08
    ##
    ## $Clusters
    ##   seqnames     start       end distance
    ## 1     chr1 165770666 165770667       11
    ## 2     chr1 165770677 165770678   377829
    ## 3     chr5  75601915  75601916       29
    ## 4     chr5  75601944  75601945 73617660
    ## 5     chr9 110235046 110235047       26
    ## 6     chr9 110235072 110235073       NA
    ## 7    chr11  32245638  32245639       95
    ## 8    chr11  32245733  32245734 63088309

# Set Enrichment Analysis

The query may be a named continuous vector. In that case, either a gene enrichment score will be calculated (if the database is discrete) or a Spearman correlation will be calculated (if the database is continuous as well). The three other cases are shown below using biologically relevant examples.

To display this functionality, let’s load two numeric database sets individually. One is a database set for CpG density and the other is a database set corresponding to the distance of the nearest transcriptional start site (TSS) to each probe.

    query <- getDBs("KYCG.MM285.designGroup")[["TSS"]]

    sesameDataCache(data_titles = c("KYCG.MM285.seqContextN.20210630"))
    res <- testEnrichmentSEA(query, "MM285.seqContextN")
    main_stats <- c("dbname", "test", "estimate", "FDR", "nQ", "nD", "overlap")
    res[,main_stats]

The estimate here is enrichment score.

> NOTE: Negative enrichment score suggests enrichment of the categorical database with the higher values (in the numerical database). Positive enrichment score represent enrichment with the smaller values. As expected, the designed TSS CpGs are significantly enriched in smaller TSS distance and higher CpG density.

Alternatively one can test the enrichment of a continuous query with discrete databases. Here we will use the methylation level from a sample as the query and test it against the chromHMM chromatin states.

    library(sesame)
    sesameDataCache(data_titles = c("MM285.1.SigDF"))
    beta_values <- getBetas(sesameDataGet("MM285.1.SigDF"))
    res <- testEnrichmentSEA(beta_values, "MM285.chromHMM")
    main_stats <- c("dbname", "test", "estimate", "FDR", "nQ", "nD", "overlap")
    res[,main_stats]

As expected, chromatin states `Tss`, `Enh` has negative enrichment score, meaning these databases are associated with small values of the query (DNA methylation level). On the contrary, `Het` and `Quies` states are associated with high methylation level.

# Session Info

    sessionInfo()

    ## R version 4.4.0 beta (2024-04-15 r86425 ucrt)
    ## Platform: x86_64-w64-mingw32/x64
    ## Running under: Windows Server 2022 x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## time zone: America/New_York
    ## tzcode source: internal
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] sesame_1.22.0               gprofiler2_0.2.3
    ##  [3] SummarizedExperiment_1.34.0 Biobase_2.64.0
    ##  [5] GenomicRanges_1.56.0        GenomeInfoDb_1.40.0
    ##  [7] IRanges_2.38.0              S4Vectors_0.42.0
    ##  [9] MatrixGenerics_1.16.0       matrixStats_1.3.0
    ## [11] knitr_1.46                  sesameData_1.21.10
    ## [13] ExperimentHub_2.12.0        AnnotationHub_3.12.0
    ## [15] BiocFileCache_2.12.0        dbplyr_2.5.0
    ## [17] BiocGenerics_0.50.0         knowYourCG_1.0.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] DBI_1.2.2               bitops_1.0-7            rlang_1.1.3
    ##  [4] magrittr_2.0.3          compiler_4.4.0          RSQLite_2.3.6
    ##  [7] png_0.1-8               vctrs_0.6.5             reshape2_1.4.4
    ## [10] stringr_1.5.1           pkgconfig_2.0.3         crayon_1.5.2
    ## [13] fastmap_1.1.1           XVector_0.44.0          utf8_1.2.4
    ## [16] rmarkdown_2.26          tzdb_0.4.0              UCSC.utils_1.0.0
    ## [19] preprocessCore_1.66.0   purrr_1.0.2             bit_4.0.5
    ## [22] xfun_0.43               zlibbioc_1.50.0         cachem_1.0.8
    ## [25] jsonlite_1.8.8          blob_1.2.4              DelayedArray_0.30.0
    ## [28] BiocParallel_1.38.0     parallel_4.4.0          R6_2.5.1
    ## [31] bslib_0.7.0             stringi_1.8.3           RColorBrewer_1.1-3
    ## [34] jquerylib_0.1.4         Rcpp_1.0.12             wheatmap_0.2.0
    ## [37] readr_2.1.5             Matrix_1.7-0            tidyselect_1.2.1
    ## [40] abind_1.4-5             yaml_2.3.8              codetools_0.2-20
    ## [43] curl_5.2.1              lattice_0.22-6          tibble_3.2.1
    ## [46] plyr_1.8.9              withr_3.0.0             KEGGREST_1.44.0
    ## [49] evaluate_0.23           Biostrings_2.72.0       pillar_1.9.0
    ## [52] BiocManager_1.30.22     filelock_1.0.3          plotly_4.10.4
    ## [55] generics_0.1.3          RCurl_1.98-1.14         BiocVersion_3.19.1
    ## [58] hms_1.1.3               ggplot2_3.5.1           munsell_0.5.1
    ## [61] scales_1.3.0            glue_1.7.0              lazyeval_0.2.2
    ## [64] tools_4.4.0             data.table_1.15.4       grid_4.4.0
    ## [67] tidyr_1.3.1             AnnotationDbi_1.66.0    colorspace_2.1-0
    ## [70] GenomeInfoDbData_1.2.12 cli_3.6.2               rappdirs_0.3.3
    ## [73] fansi_1.0.6             S4Arrays_1.4.0          viridisLite_0.4.2
    ## [76] dplyr_1.1.4             gtable_0.3.5            sass_0.4.9
    ## [79] digest_0.6.35           SparseArray_1.4.0       htmlwidgets_1.6.4
    ## [82] memoise_2.0.1           htmltools_0.5.8.1       lifecycle_1.0.4
    ## [85] httr_1.4.7              bit64_4.0.5