---
title: "YY1 ChIA-PET motif analysis (step-by-step)"
author: "Jennifer Hammelman, Konstantin Krismer"
date: "2021-07-12"
output: rmarkdown::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{YY1 ChIA-PET motif analysis (step-by-step)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```r
set.seed(17)
library(spatzie)
```

This vignette describes how to use spatzie to identify pairs of transcription factors whose sequence motifs (that describe their binding sites) are co-enriched in enhancers and promoters that interact with each other. ChIA-PET [@pmid19247990], HiChIP [@pmid27643841] or Hi-C[@pmid19815776] are molecular biology assays commonly used to investigate long-range genomic interactions and the data they generate, once properly processed (BEDPE format), serves as input to spatzie co-enrichment analyses. In this case we used interactions data in BEDPE format based on a ChIA-PET assay targeting the transcription factor YY1.

*Interactions data* in BEDPE format is a tab-separated file, where each line describes one interaction between two *anchors*, i.e., two regions of the genome that are potentially far away from each other.

In order to find motif pairs co-enriched in enhancers and promoters that interact with each other, we first need to annotate all interaction anchors and discard interactions that are not between enhancers and promoters.


# Genome assembly specific configuration

The interactions data was aligned to the mouse genome assembly *mm9* (i.e., the coordinates in the BEDPE file are *mm9* genome coordinates). The annotation package `BSgenome.Mmusculus.UCSC.mm9` contains the genomic sequence, which we need in order to detect transcription factor motifs in the interaction anchor regions. `TxDb.Mmusculus.UCSC.mm9.knownGene` contains promoter annotations.


```r
genome_id <- "BSgenome.Mmusculus.UCSC.mm9"
if (!(genome_id %in% rownames(utils::installed.packages()))) {
  BiocManager::install(genome_id, update = FALSE, ask = FALSE)
}
genome <- BSgenome::getBSgenome(genome_id)


txdb_id <- "TxDb.Mmusculus.UCSC.mm9.knownGene"
if (!(txdb_id %in% rownames(utils::installed.packages()))) {
  BiocManager::install(txdb_id, update = FALSE, ask = FALSE)
}
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene::TxDb.Mmusculus.UCSC.mm9.knownGene

ensembl_data_set <- "mmusculus_gene_ensembl"
gene_symbol <- "mgi_symbol"
```

# Load data

For this vignette we load toy BEDPE example data from a ChIA-PET experiment in murine embryonic stem cells.


```r
yy1_interactions_file <- system.file("extdata/yy1_interactions.bedpe.gz",
                                     package = "spatzie")
yy1_interactions <- GenomicInteractions::makeGenomicInteractionsFromFile(
  yy1_interactions_file,
  type = "bedpe",
  experiment_name = "yy1",
  description = "mESC yy1 chr1")
length(yy1_interactions)
```

```
## [1] 40000
```

# Load mouse promoter information

Unless the interactions data was preprocessed in such a way to include only enhancers in one anchor and only promoters in the other (or any other meaningful grouping), it is imperative to filter out non-enhancer-promoter interactions. To do this, we annotate interaction anchors as either close to a promoter or promoter-distal (i.e., putative enhancer).

Here, promoter regions are loaded, including the surrounding 2500 base pairs.


```r
promoter_ranges <- GenomicFeatures::promoters(txdb,
                                              upstream = 2500,
                                              downstream = 2500,
                                              columns = c("tx_name", "gene_id"))
# trims out-of-bound ranges located on non-circular sequences
promoter_ranges <- GenomicRanges::trim(promoter_ranges)
# remove duplicate promoters from transcript isoforms
promoter_ranges <- BiocGenerics::unique(promoter_ranges)

promoters_df <- as.data.frame(promoter_ranges)
promoters_df$gene_id <- as.character(promoters_df$gene_id)
```

# Annotate interactions with promoters

We use `annotateInteractions` from the `GenomicInteractions` package to annotate the interactions with the promoters obtained in the previous step, thus splitting the interactions into three groups: (1) interactions between promoter regions and regions far away from promoters (i.e., putative enhancers), (2) interactions between two promoter regions, and (3) interactions between two promoter-distal regions.


```r
annotation_features <- list(promoter = promoter_ranges)
GenomicInteractions::annotateInteractions(yy1_interactions, annotation_features)
GenomicInteractions::plotInteractionAnnotations(yy1_interactions)
```

![](figure/is_annotate_interactions-1.png)

# Limit and sort interactions to enhancer:promoter

As we have seen above, the interactions data contains not only interactions between enhancers and promoters. Here all interactions that are not between enhancers and promoters are discarded and enhancer-promoter interactions are sorted in a way that the promoters are always in anchor 1 and the enhancers are always in anchor 2.


```r
distal_promoter_idx <- GenomicInteractions::isInteractionType(
  yy1_interactions, "distal", "promoter")
yy1_pd <- yy1_interactions[distal_promoter_idx]
anchor1 <- GenomicInteractions::anchorOne(yy1_pd)
anchor2 <- GenomicInteractions::anchorTwo(yy1_pd)

promoter_left <- S4Vectors::elementMetadata(anchor1)[, "node.class"] == "promoter"
promoter_right <- S4Vectors::elementMetadata(anchor2)[, "node.class"] == "promoter"
promoter_ranges <- c(anchor1[promoter_left],
                     anchor2[promoter_right])
enhancer_ranges <- c(anchor2[promoter_left],
                     anchor1[promoter_right])
yy1_pd <- GenomicInteractions::GenomicInteractions(promoter_ranges,
                                                   enhancer_ranges)
```

# Scan motifs

For the purpose of this vignette, we use a toy motif database. The HOCOMOCO motif database [@Kulakovskiy2018] is commonly used, but any motif file compatible with `TFBSTools::readJASPARMatrix()` can be used.

`spatzie::scan_motifs()` takes the filtered interactions data and scans the anchors of each interaction with the motifs of the provided motif database, usually a set of transcription factor sequence motifs.

`spatzie::filter_motifs()` operates on the object returned by `spatzie::scan_motifs()` and removes motifs that are present in - in this case - fewer than 40% of the interactions.


```r
motifs_file <- system.file("extdata/motifs_subset.txt.gz",
                           package = "spatzie")
motifs <- TFBSTools::readJASPARMatrix(motifs_file, matrixClass = "PFM")
yy1_pd_interaction <- scan_motifs(yy1_pd, motifs, genome)
yy1_pd_interaction <- filter_motifs(yy1_pd_interaction, 0.4)
```

# Compute motif co-significance

`spatzie::anchor_pair_enrich()` operates on the object returned by `spatzie::filter_motifs()` and calculates the co-enrichment of all motif pairs. It supports three methods to determine co-enrichment, count-based correlation, score-based correlation, and match-based assocation. For details, see the help page (`?spatzie::anchor_pair_enrich`) or the spatzie paper (`citation("spatzie")`).

## Correlation between motif counts


```r
yy1_pd_count_corr <- anchor_pair_enrich(yy1_pd_interaction,
                                        method = "count")
pheatmap::pheatmap(-log2(yy1_pd_count_corr$pair_motif_enrich), fontsize = 6)
```

<img src="figure/is_compute_cosignificance_count-1.png" title="plot of chunk is_compute_cosignificance_count_save" alt="plot of chunk is_compute_cosignificance_count_save" width="95%" style="display: block; margin: auto;" />

## Correlation between max motif scores


```r
yy1_pd_score_corr <- anchor_pair_enrich(yy1_pd_interaction,
                                        method = "score")
pheatmap::pheatmap(-log2(yy1_pd_score_corr$pair_motif_enrich), fontsize = 6)
```

<img src="figure/is_compute_cosignificance_score-1.png" title="plot of chunk is_compute_cosignificance_score_save" alt="plot of chunk is_compute_cosignificance_score_save" width="95%" style="display: block; margin: auto;" />

YY1 binds enhancer and promoter sites providing scaffolding that forms enhancer-promoter interactions in mouse stem cells [@Weintraub2017]. As expected, spatzie identified a statistically significant co-occurrence of YY1 motifs indicating this dependency.

When interpreting spatzie results, keep in mind that motif databases such as HOCOMOCO often include groups of transcription factors with highly similar DNA-binding motifs (in this example YY1 and ZF.5), and the putative co-enrichment of one pair of transcription factor binding sites might be explained by another pair with highly similar motifs.

Please note that the motifs and the interactions data used in this vignette are dummy data used for demonstration purposes only.

# Additional information

Most of the functionality of the spatzie package is also offered through
the website at https://spatzie.mit.edu.

For a more detailed discussion on spatzie, please have a look at
the paper:

**spatzie: An R package for identifying significant transcription factor motif co-enrichment from enhancer-promoter interactions**  
Jennifer Hammelman, Konstantin Krismer, and David K. Gifford  
Nucleic Acids Research, 2022, gkac036; DOI: https://doi.org/10.1093/nar/gkac036  

# Session info

```r
sessionInfo()
```

# References
