---
title: "Epiregulon tutorial with MultiAssayExperiment"
author: "Xiaosai Yao"
email: "yao.xiaosai@gene.com"
output:
  BiocStyle::html_document:
    toc: true
    number_section: true
    self_contained: true
    titlecaps: true
package: "epiregulon"
vignette: >
  %\VignetteIndexEntry{Epiregulon tutorial with MultiAssayExperiment}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE
)
```

# Introduction

Gene regulatory networks model the underlying gene regulation hierarchies that drive gene expression and cell states. The main functions of this package are to construct gene regulatory networks and infer transcription factor (TF) activity at the single cell level by integrating scATAC-seq and scRNA-seq data and incorporating of public bulk TF ChIP-seq data.

There are three related packages: `r BiocStyle::Biocpkg("epiregulon")`, `r BiocStyle::Biocpkg("epiregulon.extra")` and `epiregulon.archr`, the two of which are available through Bioconductor and the last of which is only available through [github](https://github.com/xiaosaiyao/epiregulon.archr). The basic `epiregulon` package takes in gene expression and chromatin accessibility matrices in the form of `SingleCellExperiment` objects, constructs gene regulatory networks ("regulons") and outputs the activity of transcription factors at the single cell level. The `epiregulon.extra` package provides a suite of tools for enrichment analysis of target genes, visualization of target genes and transcription factor activity, and network analysis which can be run on the `epiregulon` output. If the users would like to start from `ArchR` projects instead of `SingleCellExperiment` objects, they may choose to use `epiregulon.archr` package, which allows for seamless integration with the [ArchR](https://www.archrproject.com/) package, and continue with the rest of the workflow offered in `epiregulon.extra`.


For full documentation of the `epiregulon` package, please refer to the [epiregulon book](https://xiaosaiyao.github.io/epiregulon.book/).

# Installation

```{r, results='hide', message=FALSE, eval=FALSE}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
 
BiocManager::install("epiregulon")

```

Alternatively, you could install from github 

```{r setup, message=FALSE, eval=FALSE}
devtools::install_github(repo ='xiaosaiyao/epiregulon')
```

Load package

```{r}
library(epiregulon)
```

# Data preparation

Prior to using `epiregulon`, single cell preprocessing needs to performed by user's favorite methods. The following components are required: <br>
1. Peak matrix from scATAC-seq containing the chromatin accessibility information <br>
2. Gene expression matrix from either paired or unpaired scRNA-seq. RNA-seq integration needs to be performed for unpaired dataset. <br>
3. Dimensionality reduction matrix from either single modality dataset or joint scRNA-seq and scATAC-seq <br>

This tutorial demonstrates the basic functions of `epiregulon`, using the reprogram-seq dataset which can be downloaded from the `r Biocpkg("scMultiome")` package. In this example, prostate cancer cells (LNCaP) were infected in separate wells with viruses encoding 4 transcription factors (NKX2-1, GATA6, FOXA1 and FOXA2) and a positive control (mNeonGreen) before pooling. The identity of the infected transcription factors was tracked through cell hashing (available in the field `hash_assignment` of the `colData`) and serves as the ground truth.

```{r}
# load the MAE object
library(scMultiome)
library(beachmat.hdf5)

mae <- scMultiome::reprogramSeq()

# extract peak matrix
PeakMatrix <- mae[["PeakMatrix"]]

# extract expression matrix
GeneExpressionMatrix <- mae[["GeneExpressionMatrix"]]
rownames(GeneExpressionMatrix) <- rowData(GeneExpressionMatrix)$name

# define the order of hash_assigment
GeneExpressionMatrix$hash_assignment <- 
  factor(as.character(GeneExpressionMatrix$hash_assignment),
         levels = c("HTO10_GATA6_UTR", "HTO2_GATA6_v2", "HTO8_NKX2.1_UTR", "HTO3_NKX2.1_v2",
                    "HTO1_FOXA2_v2", "HTO4_mFOXA1_v2", "HTO6_hFOXA1_UTR", "HTO5_NeonG_v2"))
# extract dimensional reduction matrix
reducedDimMatrix <- reducedDim(mae[['TileMatrix500']], "LSI_ATAC")

# transfer UMAP_combined from TileMatrix to GeneExpressionMatrix
reducedDim(GeneExpressionMatrix, "UMAP_Combined") <- 
  reducedDim(mae[['TileMatrix500']], "UMAP_Combined")
```
Visualize singleCellExperiment by UMAP

```{r}

scater::plotReducedDim(GeneExpressionMatrix, 
                       dimred = "UMAP_Combined", 
                       text_by = "Clusters", 
                       colour_by = "Clusters")

```

# Quick start

## Retrieve bulk TF ChIP-seq binding sites 

First, we retrieve a `GRangesList` object containing the binding sites of all the transcription factors and co-regulators. These binding sites are derived from bulk ChIP-seq data in the [ChIP-Atlas](https://chip-atlas.org/) and [ENCODE](https://www.encodeproject.org/) databases. For the same transcription factor, multiple ChIP-seq files from different cell lines or tissues are merged. For further information on how these peaks are derived, please refer to `?epiregulon::getTFMotifInfo`. Currently, human genomes hg19 and hg38 and mouse mm10 are supported. 

```{r getTFMotifInfo}
grl <- getTFMotifInfo(genome = "hg38")
grl
```

We recommend the use of ChIP-seq data over motif for estimating TF occupancy. However, if the user would like to start from motifs, it is possible to switch to the `motif` mode. The user can provide the peak regions as a `GRanges` object and the location of the motifs will be annotated based on Cisbp from chromVARmotifs (human_pwms_v2 and mouse_pwms_v2, version 0.2) 
```{r eval=FALSE}
grl.motif <- getTFMotifInfo(genome = "hg38",  
                            mode = "motif", 
                            peaks = rowRanges(PeakMatrix))
```
## Link ATAC-seq peaks to target genes

Next, we try to link ATAC-seq peaks to their putative target genes. We assign a peak to a gene within a size window (default ±250kb) if the chromatin accessibility of the peak and expression of the target genes are highly correlated (default threshold 0.5). To compute correlations, we first create cell aggregates by performing k-means clustering on the reduced dimensionality matrix. Then we aggregate the counts of the gene expression and peak matrix and average across the number of cells. Correlations are computed on the averaged gene expression and chromatin accessibility.

If cluster labels are provided, peak-to-gene correlations are computed on all the cells and for each cluster. Peak-to-gene links are retained as long as any of the correlations pass the threshold; the longer list of peak-to-gene links capture both inter- and intra-cluster variations.

We also provide a second method of target gene assignment based on the nearest expressed genes to each ATAC-seq peak. To use this function, the user can specify `assignment_method = "nearest"`. With this option, we will still calculate the correlation and apply a correlation cutoff. If no correlation filter is needed, the user can set `cor_cutoff` to 0.
```{r calculateP2G}
set.seed(1010)
p2g <- calculateP2G(peakMatrix = PeakMatrix, 
                    expMatrix = GeneExpressionMatrix, 
                    reducedDim = reducedDimMatrix,
                    exp_assay = "normalizedCounts",
                    peak_assay = "counts")

p2g

```


## Add TF motif binding to peaks

The next step is to add the TF binding information by overlapping regions of the peak matrix with the bulk chip-seq database. The output is a data frame object with three columns:

1. `idxATAC` - index of the peak in the peak matrix
2. `idxTF` - index in the gene expression matrix corresponding to the transcription factor
3. `tf` - name of the transcription factor


```{r addTFMotifInfo}

overlap <- addTFMotifInfo(grl = grl, p2g = p2g, 
                          peakMatrix = PeakMatrix)
head(overlap)
```

## Generate regulons

A DataFrame, representing the inferred regulons, is then generated. The DataFrame consists of ten columns:

1. `idxATAC` - index of the peak in the peak matrix
2. `chr` - chromosome number
3. `start` - start position of the peak
4. `end` - end position of the peak
5. `idxRNA` - index in the gene expression matrix corresponding to the target gene
6. `target` - name of the target gene
7. `distance` - distance between the transcription start site of the target gene and the middle of the peak
8. `idxTF` - index in the gene expression matrix corresponding to the transcription factor
9. `tf` - name of the transcription factor
10. `corr` - correlation between target gene expression and the chromatin accessibility at the peak. if cluster labels are provided, 
this field is a matrix with columns names corresponding to correlation across all cells and for each of the clusters.

```{r, warning=FALSE, getRegulon}

regulon <- getRegulon(p2g = p2g, overlap = overlap, 
                      aggregate = FALSE)
regulon

```

## Network pruning (highly recommended)

Epiregulon prunes the network by performing tests of independence on the observed number of cells jointly expressing transcription factor (TF), regulatory element (RE) and target gene (TG) vs the expected number of cells if TF/RE and TG are independently expressed. The users can choose between two tests, the binomial test and the chi-square test. In the binomial test, the expected probability is P(TF, RE) * P(TG), and the number of trials is the total number of cells, and the observed successes is the number of cells jointly expressing all three elements. In the chi-square test, the expected probability for having all 3 elements active is also P(TF, RE) * P(TG) and the probability otherwise is 1- P(TF, RE) * P(TG). The observed cell count for the category of "active TF" is the number of cells jointly expressing all three elements, and the cell count for the inactive category is n - n_triple.

We calculate cluster-specific p-values if users supply cluster labels. This is useful if we are interested in cluster-specific networks. The pruned regulons can then be used to visualize differential networks for transcription factors of interest. 


```{r network pruning, results = "hide", message = FALSE}

pruned.regulon <- pruneRegulon(expMatrix = GeneExpressionMatrix,
                               exp_assay = "normalizedCounts",
                               peakMatrix = PeakMatrix,
                               peak_assay = "counts",
                               test = "chi.sq",
                               regulon = regulon[regulon$tf %in% 
                                         c("NKX2-1","GATA6","FOXA1","FOXA2", "AR"),],
                               clusters = GeneExpressionMatrix$Clusters,
                               prune_value = "pval",
                               regulon_cutoff = 0.05
                               )

pruned.regulon
```

## Add Weights

While the `pruneRegulon` function provides statistics on the joint occurrence of TF-RE-TG, we would like to further estimate the strength of regulation. Biologically, this can be interpreted as the magnitude of gene expression changes induced by transcription factor activity. Epiregulon estimates the regulatory potential using one of the three measures: 1) correlation between TF and target gene expression, 2) mutual information between the TF and target gene expression and 3) Wilcoxon test statistics of target gene expression in cells jointly expressing all 3 elements vs cells that do not.

Two of the measures (correlation and Wilcoxon statistics) give both the magnitude and directionality of changes whereas mutual information is always positive. The correlation and mutual information statistics are computed on pseudobulks aggregated by user-supplied cluster labels, whereas the Wilcoxon method groups cells into two categories: 1) the active category of cells jointly expressing TF, RE and TG and 2) the inactive category consisting of the remaining cells. 

We calculate cluster-specific weights if users supply cluster labels.

```{r addWeights, results = "hide", warning = FALSE, message = FALSE}

regulon.w <- addWeights(regulon = pruned.regulon,
                        expMatrix  = GeneExpressionMatrix,
                        exp_assay  = "normalizedCounts",
                        peakMatrix = PeakMatrix,
                        peak_assay = "counts",
                        clusters = GeneExpressionMatrix$Clusters,
                        method = "wilcox")

regulon.w

```

## (Optional) Annotate with TF motifs

So far the gene regulatory network was constructed from TF ChIP-seq exclusively. Some users would prefer to further annotate regulatory elements with the presence of motifs. We provide an option to annotate peaks with motifs from the Cisbp database. If no motifs are present for this particular factor (as in the case of co-factors or chromatin modifiers), we return NA. If motifs are available for a factor and the RE contains a motif, we return 1. If motifs are available and the RE does not contain a motif, we return 0. The users can also provide their own motif annotation through the `pwms` argument.


```{r addMotifScore}

regulon.w.motif <- addMotifScore(regulon = regulon.w,
                                 peaks = rowRanges(PeakMatrix),
                                 species = "human",
                                 genome = "hg38")

# if desired, set weight to 0 if no motif is found
regulon.w.motif$weight[regulon.w.motif$motif == 0] <- 0

regulon.w.motif
```
## (Optional) Annotate with log fold changes

It is sometimes helpful to filter the regulons based on gene expression changes between two conditions. The `addLogFC` function is a wrapper around `scran::findMarkers` and adds extra columns of log changes to the regulon DataFrame. The users can specify the reference condition in `logFC_ref` and conditions for comparison in
`logFC_condition`. If these are not provided, log fold changes are calculated for every condition in the cluster labels against the rest of the conditions.
 
```{r add logFC}
regulon.w <- addLogFC(regulon = regulon.w,
                      clusters = GeneExpressionMatrix$hash_assignment,
                      expMatrix  = GeneExpressionMatrix,
                      assay.type  = "normalizedCounts",
                      pval.type = "any",
                      logFC_condition = c("HTO10_GATA6_UTR", 
                                          "HTO2_GATA6_v2", 
                                          "HTO8_NKX2.1_UTR"),
                      logFC_ref = "HTO5_NeonG_v2")

regulon.w
```
## Calculate TF activity 

Finally, the activities for a specific TF in each cell are computed by averaging expressions of target genes linked to the TF weighted by the test statistics of choice, chosen from either correlation, mutual information or the Wilcoxon test statistics. 
$$y=\frac{1}{n}\sum_{i=1}^{n} x_i * weights_i$$
where $y$ is the activity of a TF for a cell,
$n$ is the total number of targets for a TF,
$x_i$ is the log count expression of target $i$ where $i$ in {1,2,...,n} and
$weights_i$ is the weight of TF - target $i$

```{r calculateActivity}
score.combine <- calculateActivity(expMatrix = GeneExpressionMatrix,
                                   regulon = regulon.w, 
                                   mode = "weight", 
                                   method = "weightedMean", 
                                   exp_assay = "normalizedCounts",
                                   normalize = FALSE)

score.combine[1:5,1:5]
```

# Session Info

```{r}
sessionInfo()
```