---
title: "Find Matches to Canonical siRNA Seeds in Genomic Features with SeedMatchR"
author: "Tareian Cazares and Gulcin Ozer"
vignette: >
  %\VignetteIndexEntry{SeedMatchR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
abstract: >
    On-target gene knockdown using siRNA ideally results from binding 
    fully complementary regions in mRNA transcripts to induce cleavage. 
    Off-target siRNA gene knockdown can occur through several modes, one being 
    a seed-mediated mechanism mimicking miRNA gene regulation. Seed-mediated 
    off-target effects occur when the ~8 nucleotides at the 5’ end of the guide 
    strand, called a seed region, bind the 3’ untranslated regions of mRNA, 
    causing reduced translation. Experiments using siRNA knockdown paired with 
    RNA-seq can be used to detect siRNA sequences with potential off-target 
    effects driven by the seed region. SeedMatchR provides tools for exploring 
    and detecting potential seed-mediated off-target effects of siRNA in RNA-seq 
    experiments. SeedMatchR is designed to extend current differential 
    expression analysis tools, such as DESeq2, by annotating results with 
    predicted seed matches. Using publicly available data, we demonstrate the 
    ability of SeedMatchR to detect cumulative changes in differential gene 
    expression attributed to siRNA seed regions.  
---

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

## Installation

You can install the development version of SeedMatchR from [GitHub](https://github.com/) or the stable build from CRAN.

``` {r eval = FALSE}
# Install from GitHub
install.packages("devtools")
devtools::install_github("tacazares/SeedMatchR")
```

## Quick start example with public siRNA data

```{r include = FALSE}
# Import library
library(SeedMatchR)
library(msa)
library(GenomicFeatures)
```

### Introduction to example dataset

This example uses the siRNA sequence, D1, targeting the Ttr gene in rat liver from the publication:

```
Schlegel MK, Janas MM, Jiang Y, Barry JD, Davis W, Agarwal S, Berman D, Brown CR, Castoreno A, LeBlanc S, Liebow A, Mayo T, Milstein S, Nguyen T, Shulga-Morskaya S, Hyde S, Schofield S, Szeto J, Woods LB, Yilmaz VO, Manoharan M, Egli M, Charissé K, Sepp-Lorenzino L, Haslett P, Fitzgerald K, Jadhav V, Maier MA. From bench to bedside: Improving the clinical safety of GalNAc-siRNA conjugates using seed-pairing destabilization. Nucleic Acids Res. 2022 Jul 8;50(12):6656-6670. doi: 10.1093/nar/gkac539. PMID: 35736224; PMCID: PMC9262600.
```

The guide sequence of interest is 23 bp long and oriented 5' -> 3'. 

```{r}
# siRNA sequence of interest targeting a 23 bp region of the Ttr gene
guide.seq = "UUAUAGAGCAAGAACACUGUUUU"
```


### Required Input Data

The required inputs to `SeedMatchR` are a DESeq2 results data frame in addition to species specific annoation data such as GTF, 2bit DNA, and 3' UTR sequences. 

`SeedMatchR` makes extensive use of `AnnotationDB` objects to help access genomic information in a reproducible manner. The required annotations are: 

*  A character string representing the siRNA RNA sequence. This must be greater than 8 bp.
* `res`: a data frame of DESEQ2 results.
* `GTF`: gene transfer file containing species specific genomic information for gene bodies. This is used to derive the list of 3' UTRs and other features used in the analysis. This is also used to map transcript IDs to gene IDs.
* Feature-specific `DNAStringSet`: A `DNAStringSet` object of sequences for each of the features of interest. The features must be named according to the transcript they were derived from. Examples include those generated by `GenomicFeatures::extractTranscriptSeqs()` paired with functions like `GenomicFeatures::threeUTRsByTranscript()`.

#### Prepare species-specific annotation data

The function `load_species_anno_db()` has built in annotation data for human, rat, and mouse annotations. We can load the species specific annotations using the following approach:

```{r echo = T, results = 'hide', message=FALSE, warning=FALSE, error=FALSE}
# Load the species specific annotation database object
anno.db <- load_species_anno_db("rat")
```

#### Extract features and sequences of interest from annotations

We will use the annotations to derive the features and feature sequences that we want to scan for each gene. 

```{r}
features = get_feature_seqs(anno.db$tx.db, anno.db$dna, feature.type = "3UTR")
```

#### Prepare DESEQ2 Results

`SeedMatchR` assumes that you will be performing your analysis on DESEQ2 results outputs. The first step is to load your DESEQ2 results file as a data frame. 

The test data that is provided with `SeedMatchR` was derived from the 2022 publication by Schlegel et al. The data set represents a DESeq2 analysis performed on rat liver that had been treated with Ttr targeting siRNA. We will use this example to explore seed mediated activity. The data set name is long, so it will be renamed to `res`.

#### Download data (only need to perform once)

We start by downloading the example data set. This function will download three files from the GEO accession [GSE184929](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE184929). These files represent three samples with different siRNA treatments at two dosages. 

```{r echo = T, results = 'hide', message=FALSE, warning=FALSE, error=FALSE}
get_example_data("sirna")
```

#### Load example data

We can load the example data into the environment.

```{r}
sirna.data = load_example_data("sirna")
```

The DESeq2 results are available through the names `Schlegel_2022_Ttr_D1_30mkg`, `Schlegel_2022_Ttr_D4_30mkg` and `Schlegel_2022_Ttr_D1_10mkg`. The data set name is long, so it will be renamed to `res`.

```{r}
res <- sirna.data$Schlegel_2022_Ttr_D1_30mkg
```

The DESeq2 results file is then filtered. The function `filter_deseq()` can be used to filter a results file by log2FoldChange, padj, baseMean, and remove NA entries. 

```{r}
# Dimensions before filtering

dim(res) # [1] 32883    6

# Filter DESeq2 results for SeedMatchR
res = filter_deseq(res, fdr.cutoff=1, fc.cutoff=0, rm.na.log2fc = TRUE)

# Dimensions after filtering
dim(res) # [1] 13582     8
```

### Plot possible seeds

Use the `plot_seeds()` function to visualize the available SeedMatchR options for your input
sequence of interest. The only input to `plot_seeds()` is your input siRNA sequence of interest. 
This function assumes that you are using a RNA input. 

```{r fig.height=5, fig.width=10, out.retina=1}
# Plot the seed sequence options for the siRNA of interest
avail.seed.plot = plot_seeds(guide.seq)

avail.seed.plot
```


### Get the seed sequence of interest

You can extract the seed sequence information from the siRNA input sequence 
using the `get_seed()` function. The inputs to the `get_seed()` function in the 
siRNA sequence of interest and the name of the seed.

```{r}
# Get the seed sequence information for the seed of interest
seed = get_seed(guide.seq, "mer7m8")

seed
```

### Counting seed matches in transcripts

You can perform a seed match for a single seed using the `SeedMatchR()`
function. 

```{r}
res = SeedMatchR(res, 
                 anno.db$gtf, 
                 features$seqs, 
                 guide.seq)

head(res)
```

#### Match multiple seeds

You can perform seed matching for all available seeds using a for loop. The results
will be appended as a new column to the results data frame. 

```{r}
for (seed in c("mer8", "mer6", "mer7A1")){
res <- SeedMatchR(res, 
                  anno.db$gtf, 
                  features$seqs, 
                  guide.seq, 
                  seed.name = seed)
}

head(res)
```

##### Match seeds with mismatches and indels allowed

You can also allow for inexact seed matches in your analysis with the `mismatches` and `indels` arguments. The names can be adjusted to reflect the arguments using the `col.name` argument. 

```{r}
for (indel.bool in c(TRUE, FALSE)){
  for (mm in c(0,1,2)){
    for (seed in c("mer7m8", "mer8", "mer6", "mer7A1")){
      res <- SeedMatchR(res, 
                        anno.db$gtf, 
                        features$seqs, 
                        guide.seq, 
                        seed.name = seed, 
                        col.name = paste0(seed, ".", "mm", mm, "_indel", indel.bool), 
                        mismatches = mm, 
                        indels = indel.bool)
    }
  }
}

head(res)
```

### Comparing the expression profiles of seed targets to background

Many factors that perturb gene expression, like miRNA, show cumulative changes in their targets gene expression. Cumulative changes in the profile of genes expression can be visualized and tested with the emperical distribution function (ecdf) coupled with a statistical test such as the Kolmogorov-Smirnov test. 

`SeedMatchR` provides functions for comparing the log2(Fold Change) of two gene sets. The function `deseq_fc_ecdf` is designed to work directly with a DESeq2 results data frame.

Required Inputs:

- `res`: DESeq2 results data frame
- `gene.lists`: A list of lists containing gene names

```{r fig.height=5, fig.width=10, out.retina=1}
# Gene set 1 
mer7m8.list = res$gene_id[res$mer7m8.mm0_indelFALSE >= 1 & res$mer8.mm0_indelFALSE ==0]

# Gene set 2 
mer8.list = res$gene_id[res$mer8.mm0_indelFALSE >= 1]

background.list = res$gene_id[res$mer7m8.mm0_indelFALSE == 0 & res$mer8.mm0_indelFALSE == 0]

ecdf.results = deseq_fc_ecdf(res, 
                             list("Background" = background.list, "mer8" = mer8.list, "mer7m8" = mer7m8.list),
                             stats.test = "KS", 
                             factor.order = c("Background", "mer8", "mer7m8"), 
                             null.name = "Background",
                             target.name = "mer8", 
                             alternative = "greater")

ecdf.results$plot
```

### Using `SeedMatchR` to explore potential small activating RNA effects

```{r fig.height=5, fig.width=10, out.retina=1}
# Group transcripts by gene
sequences <- transcriptsBy(anno.db$tx.db, by="gene")

# Extract promoter sequences from tx.db object
prom.seq = getPromoterSeq(sequences,
               anno.db$dna,
                 upstream=2000,
                 downstream=100)

# perform a seed search of the promoter sequences. Set tx.id.col to F to use gene annotations
res = SeedMatchR(res, anno.db$gtf, prom.seq@unlistData, guide.seq, tx.id.col = FALSE, col.name = "promoter.mer7m8")

# Find the genes with matches
promoterWseed = res$gene_id[res$promoter.mer7m8 >= 1]

# Generate the background list of genes
background.list = res$gene_id[!(res$gene_id %in% promoterWseed)]

# Plot ecdf results for promoter matches with stats testing
ecdf.results = deseq_fc_ecdf(res, 
                             title = "Ttr D1 30mkg",
                             list("Background" = background.list, 
                                  "Promoter w/ mer7m8" = promoterWseed),
                             stats.test = "KS", 
                             factor.order = c("Background", 
                                              "Promoter w/ mer7m8"), 
                             null.name = "Background",
                             target.name = "Promoter w/ mer7m8",
                             alternative = "less",
                             palette = c("black", "#d35400"))

ecdf.results$plot

```

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