---
title: "Tree-based differential testing using inferential uncertainty for RNA-Seq"
author: "Noor Pratap Singh"
date: "`r format(Sys.Date(), '%m/%d/%Y')`"
output:
  rmarkdown::html_document:
    highlight: tango
    toc: true
    toc_float: true
vignette: |
  %\VignetteIndexEntry{1. Tree-based differential testing using inferential uncertainty for RNA-Seq}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, echo=FALSE}
knitr::opts_chunk$set(tidy=FALSE, cache=FALSE, dev="png",
                      message=FALSE, error=FALSE, warning=FALSE,
		      eval=TRUE, collapse=TRUE, comment = "#>")
```
# Installation

```{r eval=FALSE}
devtools::install_github("NPSDC/mehenDi")
devtools::install_github("NPSDC/brainSimSmall") ## Download the package with a small dataset
suppressPackageStartupMessages(library(mehenDi))
```
We also install the package *brainSimSmall* which is a small RNA-Seq dataset, on which this example is run.

# mehenDi
The *mehenDi* method performs differential testing use the transcripts-tree obtained from *TreeTerminus*.
The transcripts are arranged in a way such that the uncertainty decreases on ascending the tree, with the
inner nodes comprising transcripts groups. *mehenDi* aims to leverage the tree structure to find differential
signal in the tree which might get lost due to uncertainty. The selected nodes can consist of both transcripts
and inner nodes, with inner nodes representing a set of transcripts that exhibit high uncertainty.

# Quick start

We demonstrate *mehenDi* on a small simulated RNA-Seq dataset in the [brainSimSmall](https://github.com/NPSDC/brainSimSmall) package.
We have *Salmon* quantified counts for the samples, with 100 inferential replicates generated using Gibbs Sampling. *Note* Gibbs Sampling
can be replaced by bootstrap sampling. The dataset consists of 12 samples belonging to two-groups, with 6 samples in each.

## Loading data 
We first load the data using [beaveR](https://github.com/NPSDC/beaveR) to create *TreeSummarizedExperiment* object. 
The leaves consist of the transcripts.
```{r eval=TRUE, message=FALSE, warning=FALSE}
dir <- system.file("extdata", package="brainSimSmall")
clustFile <- file.path(dir, "cluster_nwk.txt")
quantDir <- file.path(dir, "salmon_quants")
samples <- as.vector(outer(c(1:6), c(1,2), function(x,y) paste(x,y,sep='_')))
quantFiles <- file.path(quantDir, samples, 'quant.sf')
coldata <- data.frame(files=quantFiles, names=samples, condition=factor(rep(1:2, each=6)))
tse <- beaveR::buildTSE(treeTermFile = clustFile, coldata = coldata)
tse
tree <- TreeSummarizedExperiment::rowTree(tse)
tree
l <- length(tree$tip)
```

## Differential transcript expression
We first perform differential testing on the transcripts using Swish. We can see the transcripts
that are called significant at the 0.01 FDR threshold.
```{r eval=TRUE, message=FALSE, warning=FALSE}
set.seed(1)
yTxps <- tse[1:l,]
yTxps <- fishpond::scaleInfReps(yTxps, quiet=T)
yTxps <- fishpond::labelKeep(yTxps)
yTxps <- fishpond::swish(yTxps, x="condition", quiet=T)
yTxps <- yTxps[match(rownames(yTxps), tree$tip.label),]
table(mcols(yTxps)[["qvalue"]] < 0.1)
```

## Differential gene expression
```{r eval=TRUE, message=FALSE, warning=FALSE}
yAllTxps <- tximeta::tximeta(coldata)
gy <- tximeta::summarizeToGene(yAllTxps)
gy <- fishpond::scaleInfReps(gy, quiet=TRUE)
gy <- fishpond::labelKeep(gy)
gy <- fishpond::swish(gy, x="condition", quiet=TRUE)
dim(gy)
table(mcols(gy)[["qvalue"]] < 0.1)
```

## Running mehenDi
We next run *mehenDi* on the TreeSummarizedExperiment object *tse* created above. The pipeline is similar to Swish, where first
the data is normalized (we used median-ratio scaling in this example). We have used *Swish* for computing the p-values, however we compute them separately for
the leaves and the inner nodes. This is because due to lower uncertainty, the width of the distribution of the 
test-statistic of inner nodes is wider and thus have smaller p-values compared to the inner nodes.

*mehenDi* outputs a set of selected nodes based on a filtering criteria. Once a node
that has a p-value less than a certain threshold (governed by \strong(alpha)) is found, we check for
different criteria. First, nodes should consist of atleast one child node that has a p-value larger than
a threshold, otherwise we are overaggregating. We then check if all the underlying descendant nodes that have a 
confident direction of sign-change have the same direction change between conditions. The confidence is determined
by the \strong(minP) parameters and looks at the inferential replicates. Finally, we also want that the child nodes of a selected node, to have uncertain abundance estimates (governed by \strong(meanInfRV)).

```{r eval=TRUE, message=FALSE, warning=FALSE}
suppressPackageStartupMessages(library(mehenDi))
set.seed(1)
yAll <- beaveR::computeSizeFactors(tse, quiet=TRUE)
yAll <- beaveR::scInfReps(yAll, quiet=TRUE)
yAll <- fishpond::labelKeep(yAll)
yInn <- yAll[(l+1):nrow(yAll),]
yInn <- fishpond::swish(yInn, x="condition", quiet=TRUE)

pvals <- c(mcols(yTxps)[["pvalue"]], mcols(yInn)[["pvalue"]])
tD <- mehenDi(yAll, x="condition", pvalues = pvals,
                    minP=0.7, mIRVThresh=0.4, alpha=0.1)
length(tD$selNodes) ## Total number of selected nodes
sum(tD$selNodes > l) ## Number of inner nodes that are selected
table(ape::node.depth(tree, 2)[tD[["selNodes"]]]) ## Distribution of the height of the nodes output by mehenDi
```

## Visualizing the node output by mehenDi
We next visualize a node that is output by *mehenDi*. We plot the transcripts in the pileup style corresponding to the gene that is uniquely mapped by the transcripts for that node and the corresponding transcript-tree. We also plot the inferential replicates between the two conditions, for both the transcript with the most significant p-value and the node output by *mehenDi*. While a lot of different libraries and functions are utilized to create this plot, we use `plotGarderner` library to organize all the plots together.


### Transcript tree for a mehenDi node
We first create a tree for a selected *mehenDi* inner node using the library `ggtree`

```{r}
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggtree))
node <- tD$selNodes[tD$selNodes > l][7]
treeSub <- tidytree::tree_subset(tree, node, levels_back = 0)
txShow <- rownames(tse)[phangorn::Descendants(tree,node)[[1]]] ## Transcript names for the transcripts in the tree
txpMin <- treeSub$tip.label[which.min(mcols(yTxps)[treeSub$tip.label, "pvalue"])]
minTInd <- match(txpMin, treeSub$tip)

size <- 5
fs <- 20
nI <- 5

pTree <- ggtree(treeSub) +
        xlim(NA, 80) +
        geom_tiplab(size=6, hjust=-0.1) +
        geom_point2(aes(subset=(node==nI), color="red"),
                size=5, fill='red', show.legend=T) +
        geom_point2(aes(subset=(node==minTInd), color = "black"),
                size=5, fill='black', show.legend=T) +
        geom_cladelab(node = nI, label = "", textcolour="red", barsize=2,
                      barcolour="red",  fontsize=5, offset = 36) +
        theme(legend.position = "bottom", legend.text=element_text(size=fs),
             plot.title=element_text(size=fs, face="bold", hjust=0)) +
        scale_color_manual(name = "", labels=c("Transcript with the lowest p-value", "Selected Node"),
               values=c("black", "red"))
```

### Plots for inferential replicates
We next store the plots for the inferential replicates using `fishpond`'s *plotInfReps* function for both the mehenDi node and its underlying transcript that has the lowest p-value.
```{r}
suppressPackageStartupMessages(library(ggplotify))
cex <- 1.6
legPos <- "topright"
pTxp <- as.grob(function() fishpond::plotInfReps(yAll, txpMin, x = "condition", legend=TRUE,
                              main=txpMin, legendTitle=TRUE, legendCex=cex,	
                             legendPos = legPos))

pInn <- as.grob(function() fishpond::plotInfReps(yAll, node, x = "condition", legend=TRUE,
                              main=paste0("Node",node), legendTitle=TRUE, legendCex=cex,	
                             legendPos = legPos))
```

### Pileup style plot and putting all the plots together
We use the *plotTranscripts* function in the `plotGardener` package to create the pileup style plot. However, we have to go through a series of steps before we can make use of this function.

#### Extract gene, transcripts and their ranges in txdb for a given node
We use *txdb* and load the reference corresponding to *hg38* and then accordingly extract the information corresponding to the *mehenDi* node.
```{r}
suppressPackageStartupMessages(library(org.Hs.eg.db))
suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene))
suppressPackageStartupMessages(library(plyranges))

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
txpHuman <- transcripts(txdb)

txpHuman <- txpHuman %>%
  mutate(tx_id = as.character(tx_id))

txpHuman <- txpHuman %>%
  mutate(gene_id = mapIds(
           txdb, keys=tx_id,
           column="GENEID", keytype="TXID")
         ) %>%
  filter(!is.na(gene_id))
txpHuman <- txpHuman %>%
  mutate(ensID = mapIds(
           org.Hs.eg.db, keys=gene_id,
           column="ENSEMBL", keytype="ENTREZID")
         ) %>%
  filter(!is.na(ensID))
txpHuman <- txpHuman %>%
                filter(tx_name %in% intersect(tree$tip, mcols(txpHuman)[["tx_name"]]))


genesHuman <- genes(txdb)
genesHuman <- genesHuman %>%
  mutate(ensID = mapIds(org.Hs.eg.db,
                            gene_id, "ENSEMBL", "ENTREZID")) %>%
    filter(!is.na(ensID))

gene <- txpHuman %>%
        filter(tx_name %in% treeSub$tip) %>%
        tibble::as_tibble() %>%
        dplyr::select(ensID) %>%
        unlist %>%
        unique

## We have selected the transcripts only in the positive strand for this gene, since the transcripts corresponding to mehenDi node belong to this strand and plotting transcripts for both strands, would make the figure very big
gDB <- genesHuman %>% 
        filter(ensID == gene, strand == "+") %>%
        tibble::as_tibble()
```


#### Putting it all together
```{r, fig.width=17, fig.height=17}
suppressPackageStartupMessages(library(plotgardener))
suppressPackageStartupMessages(library(grid))

pageCreate(width = 17, height = 17, showGuides = FALSE, default.units = "inches")

par <- pgParams(
          chrom = as.character(gDB[["seqnames"]]), 
          chromstart = gDB[["start"]], chromend = gDB[["end"]],
          assembly = "hg38", just = c("left", "bottom"), 
        default.units = "inches", fontsize = fs-7
        )
    
hilite <- data.frame(transcript=txShow, 
                     color=rep("red", length(txShow))
                             )
parText <- pgParams(
    just = c("left"), fontsize = fs
)#, fontface = "bold")

plotText(label="A) Transcripts mapping to gene PDE4DIP", x = 0.5, y = 0.5,
        params = parText)

plotTranscripts(
    params = par, x = 1.5, y = 9.8, width = 5, height = 9, 
    spaceHeight=1, boxHeight=0.1,
#     labels = "transcript", fill="blue",
    transcriptHighlights = hilite
)

plotGenomeLabel(
  params = par, x = 0.5, y = 10.1, length=6,
  just = c("left", "top"), fontsize=fs
)


plotText(label="B) Tree representing the transcripts covered by gene\n PDE4DIP", 
         x = 8, y = 0.5, params = parText)
plotGG(
    plot = pTree,
    x = 7.7, y = 2, width = 8, height=6
)

plotText(label="C) Scaled inferential replicates for the \ntranscript ENST00000369356.8", 
         x = 0.5, y = 11, params = parText)

plotGG(
    plot = pTxp,
    x = 0.5, y = 11.5, width = 6, height=6
)

plotText(label="D) Scaled inferential replicates for the selected node", 
         x = 7.8, y = 11, params = parText)

plotGG(
    plot = pInn,
    x = 8.3, y = 11.5, width = 6, height=6
)

```

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