---
title: "Swish: differential expression accounting for inferential uncertainty"
date: "`r format(Sys.Date(), '%m/%d/%Y')`"
output:
  rmarkdown::html_document:
    highlight: tango
    toc: true
    toc_float: true
bibliography: library.bib
vignette: |
  %\VignetteIndexEntry{1. Swish: DE analysis accounting for inferential uncertainty}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- run this document with rmarkdown::render("swish.Rmd") -->

```{r setup, echo=FALSE, results="hide"}
knitr::opts_chunk$set(tidy=FALSE, cache=FALSE, dev="png",
                      message=FALSE, error=FALSE, warning=FALSE)
```

# The Swish method

The *Swish* method for differential expression analysis of bulk or
single-cell RNA-seq data using inferential replicate counts is
described in the following reference: 
@swish [doi: 10.1093/nar/gkz622](https://doi.org/10.1093/nar/gkz622).

We note that *Swish* extends and builds on another method, *SAMseq*
[@samseq], implemented in the *samr* package,
by taking into account inferential uncertainty, and allowing to
control for batch effects and matched samples. Additionally, *Swish*
has methods for testing changes in effect size across secondary
covariates, which we refer to as "interactions". *Swish* also
implements correlation tests as in the original *SAMseq* package.
`swish` calls functions from the *qvalue* [@qvalue] or *samr* package
for calculation of local FDR and q-value. This vignette gives an
example of differential analysis of matched samples, and an
interaction test for matched samples, to see if a condition effect
changes in magnitude across two groups of samples.

# Quick start

The following lines of code will perform a basic transcript-level
`swish` two group analysis of bulk RNA-seq. For more details, read
on. There is a special section below for two-group analysis
of [scRNA-seq](#alevin_scrna-seq).

```{r eval=FALSE}
# 'coldata.csv': sample information table
coldata <- read.csv("coldata.csv")
library(tximeta)
y <- tximeta(coldata) # reads in counts and inf reps
library(fishpond)
y <- scaleInfReps(y) # scales counts
y <- labelKeep(y) # labels features to keep
set.seed(1)
y <- swish(y, x="condition") # simplest Swish case
```

**Note:** Inferential replicates, either from Gibbs sampling or
bootstrapping of reads, are required for the *swish* method. 
When running Salmon, you can set `--numGibbsSamples 30` to generate
Gibbs samples, or `--numBootstraps 30` to generate bootstrap samples
(we typically recommend 20-30 inferential replicates).

The results can be found in `mcols(y)`. For example, one can calculate
the number of genes passing a 5% FDR threshold:

```{r eval=FALSE}
table(mcols(y)$qvalue < .05)
```

One can at any point remove the genes that didn't pass the expression
filter with the following line of code (can be run before or after
`swish`). These genes are ignored by `swish`, and so will have `NA` in
the results columns in `mcols(y)`.

```{r eval=FALSE}
y <- y[mcols(y)$keep,]
```

A gene-level analysis looks identical to a transcript-level analysis,
only the input data changes. Examples follow.

Lastly, what is the structure of the output of `tximeta` [@tximeta], which is
used in `swish`? See the section below, *Structure of tximeta
output / swish input*.

## Macrophage stimulation experiment

We begin the *fishpond* vignette by loading data from a Bioconductor
Experiment Data package, *macrophage*. The package contains RNA-seq
quantification from 24 RNA-seq samples, which are a subset of the
RNA-seq samples generated and analyzed by @alasoo - 
[doi: 10.1038/s41588-018-0046-7](https://doi.org/10.1038/s41588-018-0046-7).

The experiment involved treatment of macrophage cell lines from a number
of human donors with IFN gamma, *Salmonella* infection, or both
treatments combined. In the beginning of this vignette, we will focus
on comparing the IFN gamma stimulated cell lines with the control cell
lines, accounting for the paired nature of the data (cells from the
same donor). Later in the vignette we will analyze differences in the
*Salmonella* infection response by IFN gamma treatment status --
whether the cells are primed for immune response.

We load the package, and point to the `extdata` directory. For a
typical analysis, the user would just point `dir` to the location on
the machine or cluster where the transcript quantifications are stored
(e.g. the `quant.sf` files).

```{r}
library(macrophage)
dir <- system.file("extdata", package="macrophage")
```

The data was quantified using *Salmon* [@salmon] 0.12.0 against the
Gencode v29 human reference transcripts [@gencode]. For more details
and all code used for quantification, refer to the
[macrophage](https://bioconductor.org/packages/macrophage) 
package vignette. 

Importantly, `--numGibbsSamples 20` was used to generate 20
inferential replicates with *Salmon*'s Gibbs sampling procedure.
We also 
recommend to use `--gcBias` when running *Salmon* to protect against
common sample-specific biases present in RNA-seq data.

# Data import

## Read in the column data from CSV

We start by reading in a CSV with the *column data*, that is,
information about the samples, which are represented as columns of
the *SummarizedExperiment* object we will construct containing the
counts of reads per gene or transcript.

```{r}
coldata <- read.csv(file.path(dir, "coldata.csv"))
head(coldata)
```

We will subset to certain columns of interest, and re-name them for
later.

```{r}
coldata <- coldata[,c(1,2,3,5)]
names(coldata) <- c("names","id","line","condition")
```

## Add a column pointing to your files

`coldata` needs to have a column `files` which specifies the path to
the quantification files. In this case, we've gzipped the
quantification files, so we point to the `quant.sf.gz` file. We make
sure that all the files exist in the location we specified.

```{r}
coldata$files <- file.path(dir, "quants", coldata$names, "quant.sf.gz")
all(file.exists(coldata$files))
```

## Read in quants with tximeta

We will read in quantification data for some of the samples. First we
load the *SummarizedExperiment* package. We will store out data and
the output of the statistical method in a *SummarizedExperiment*
object. We use the *tximeta* [@tximeta] package to read in the data:

```{r}
suppressPackageStartupMessages(library(SummarizedExperiment))
```

```{r include=FALSE}
# This hidden code chunk is only needed for Bioc build machines,
# so that 'fishpond' will build regardless of whether
# the machine can connect to ftp.ebi.ac.uk.
# Using linkedTxomes to point to a GTF that lives in the macrophage pkg.
# The chunk can be skipped if you have internet connection,
# as tximeta will automatically ID the transcriptome and DL the GTF.
library(tximeta)
makeLinkedTxome(
  indexDir=file.path(dir, "gencode.v29_salmon_0.12.0"),
  source="myGENCODE",
  organism="Homo sapiens",
  release="29",
  genome="GRCh38",
  fasta="ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.transcripts.fa.gz",
  gtf=file.path(dir, "gencode.v29.annotation.gtf.gz"), # local version
  write=FALSE
)
# 'myGENCODE' is used to avoid tximeta() attempting to
# pull down the TxDb from AnnotationHub, alterantively
# once can specify useHub=FALSE.
```

The *macrophage* package contains a number of samples under different 
conditions. In this vignette, we will  demonstrate a two group 
comparison, and so we first subset our `coldata` to the `"naive"` 
and `"IFNg"` groups. Log fold changes will be made comparing 
`IFNg` to `naive` (the reference level).

```{r}
coldata <- coldata[coldata$condition %in% c("naive","IFNg"),]
coldata$condition <- factor(coldata$condition, 
                            levels=c("naive","IFNg"))
```

We load in the quantification data with `tximeta`:

```{r}
library(tximeta)
se <- tximeta(coldata)
```

We can see that all the assays have been loaded:

```{r}
assayNames(se)
```

`tximeta` loads transcript-level data, although it can later be
summarized to the gene levels:

```{r}
head(rownames(se))
```

```{r include=FALSE}
# hidden chunk to trick tximeta into knowing this is GENCODE,
# this is only needed because above we labelled the txome myGENCODE
metadata(se)$txomeInfo$source <- "GENCODE"
```

We will rename our *SummarizedExperiment* `y` for the statistical
analysis. For speed of the vignette, we subset to the transcripts on
chromosome 4.

```{r}
y <- se
y <- y[seqnames(y) == "chr4",]
```

**Note on factor levels**: 
The `swish` function compares expression level across factors such
that log2 fold changes are reported as the non-reference level over
the reference level.  By default, R will choose a *reference level*
for factors based on alphabetical order, unless `levels` are
explicitly set. It is recommended to set the factors levels, as in
the above code chunk, with the reference level coming first in the
character vector, so that log2 fold changes correspond to the
comparison of interest.

# Differential transcript expression

## Running Swish at the transcript level

Running `swish` has three steps: scaling the inferential replicates,
labeling the rows with sufficient counts for running differential
expression, and then calculating the statistics. As `swish` makes use
of pseudo-random number generation in breaking ties and in calculating
permutations, to obtain identical results, one needs to set a random
seed before running `swish()`, as we do below.

```{r results="hide", message=FALSE}
library(fishpond)
y <- scaleInfReps(y)
y <- labelKeep(y)
y <- y[mcols(y)$keep,]
set.seed(1)
y <- swish(y, x="condition", pair="line")
```

**Note:** the simple paired test is the slowest of all the other
designs provided for Swish
(see [further details](#further-details)), because it requires
recomputing the signed ranks at each permutation. Other designs avoid
recomputing ranks per permutation, or don't use ranks for the test
statistic. You can set `fast=1` for simple paired designs, 
which implements a one-sample z-score as the test statistic,
and permutations to provide p-values and q-values. 
These permutations are much faster than those that use the signed 
rank test, for example with one core, the signed rank test is 12x 
slower than the one-sample z-score method when run across all
transcripts, and they nevertheless detect nearly the same sets.

```{r}
y_fast <- swish(y, x="condition", pair="line", fast=1)
table(original=mcols(y)$qvalue < .1,
      fast=mcols(y_fast)$qvalue < .1)
```

The default number of permutations for computing p-values is
`nperms=100`. However, for paired datasets as this one, you may have
fewer maximum permutations. In this case, there are 64 possible ways
to switch the condition labels for six pairs of samples. We can set
the `nperms` manually to 64, or if we had just used the default value,
`swish` would set `nperms` to the maximum value possible and notify
the user that it had done so.

A note about `labelKeep`: by default we keep features with `minN=3`
samples with a minimal count of 10. For scRNA-seq data with
de-duplicated UMI counts, we recommend to lower the count, e.g. a
count of 3, across a higher number of `minN` cells, depending on the
number of cells being compared. You can also set `x="condition"` when
running `labelKeep` which will use the condition variable to set
`minN`.

The results are stored in `mcols(y)`. We will show below how to pull
out the top up- and down-regulated transcripts.

We can see how many transcripts are in a 5% FDR set:

```{r}
table(mcols(y)$qvalue < .05)
```

A number of features may have the same `pvalue` and `qvalue` due to
the procedure for assessing significance. We can additionally rank
features by their effect size, to break ties in the `qvalue`:

```{r}
most.sig <- with(mcols(y),
                 order(qvalue, -abs(log2FC)))
mcols(y)[head(most.sig),c("log2FC","qvalue")]
```

The top 6 genes by `qvalue` and effect size in this case all have
positive LFC, although we have ranked by the largest in absolute
value, so large negative values will also appear in this list.

## Plotting results

We can check the distribution of p-values. This looks as expected for
a comparison where we expect many transcripts will be affected by the
treatment (IFNg stimulation of macrophage cells). There is a flat
component and then an enrichment of transcripts with p-values near 0.

```{r}
hist(mcols(y)$pvalue, col="grey")
```

Of the transcripts in this set, which have the most extreme log2
fold change? Note that often many transcripts will share the same
q-value, so it's valuable to look at the log2 fold change as well (see
further note below on q-value computation). The log2 fold
change computed by `swish` is the median over inferential replicates,
and uses a pseudo-count of 5 on the scaled counts, to stabilize
the variance on the fold change from division by small counts. 
See the note above for interpretation of log2 fold changes with
respect to the levels of `x`.

With the following code chunk, we construct two vectors that give the
significant genes with the lowest (most negative) and highest (most
positive) log2 fold changes.

```{r}
with(mcols(y),
     table(sig=qvalue < .05, sign.lfc=sign(log2FC))
     )
sig <- mcols(y)$qvalue < .05
lo <- order(mcols(y)$log2FC * sig)
hi <- order(-mcols(y)$log2FC * sig)
```

Here we print a small table with just the calculated statistics for
the large positive log fold change transcripts (up-regulation):

```{r}
top_up <- mcols(y)[head(hi),]
names(top_up)
cols <- c("log10mean","log2FC","pvalue","qvalue")
print(as.data.frame(top_up)[,cols], digits=3)
```

Likewise for the largest negative log fold change transcripts
(down-regulation): 

```{r}
top_down <- mcols(y)[head(lo),]
print(as.data.frame(top_down)[,cols], digits=3)
```

We can plot the scaled counts for the inferential replicates, and also
group the samples by a covariate, in this case the cell line. The
analysis was paired, so the statistic assessed if the change within
pairs was consistent. Here we plot the 100th top up-regulated
transcript. `plotInfReps` also has functionality for plotting
uncertainty in single cell data, as will be covered in a later
section. 

```{r}
plotInfReps(y, idx=hi[100], x="condition", cov="line")
```

We can make an 
[MA plot](https://en.wikipedia.org/wiki/MA_plot), 
where the transcripts in our FDR set are colored:

```{r}
plotMASwish(y, alpha=.05)
```

Using the `addIds` function from *tximeta*, we can easily add gene
symbols. By specifying `gene=TRUE`, this will use the gene ID to match
to gene symbols for all of the transcripts.

```{r}
library(org.Hs.eg.db)
y <- addIds(y, "SYMBOL", gene=TRUE)
```

We can then add gene symbols to our MA plot:

```{r}
plotMASwish(y, alpha=.05, xlim=c(.5,5.5))
with(
  subset(mcols(y), qvalue < .05 & abs(log2FC) > 4),
     text(log10mean, log2FC, SYMBOL,
          col="blue", pos=4, cex=.7)
)
```

# Differential gene expression

## Running Swish at the gene level

We can also run swish at the gene level. First we summarize all of the
data to the gene level, using the `summarizeToGene` function from
*tximeta*. Again, we rename the object for statistical analysis, and
then we subset to the genes on chromosome 4 for the demonstration.

```{r}
gse <- summarizeToGene(se)
gy <- gse
gy <- gy[seqnames(gy) == "chr4",]
```

Next we can run the same steps as before. Again we set a random seed
in order to be able to reproduce exact results in the future:

```{r results="hide", message=FALSE}
gy <- scaleInfReps(gy)
gy <- labelKeep(gy)
gy <- gy[mcols(gy)$keep,]
set.seed(1)
gy <- swish(gy, x="condition", pair="line")
```

As before, the number of genes in a 1% FDR set:

```{r}
table(mcols(gy)$qvalue < .05)
```

## Plotting gene results

The histogram of p-values:

```{r}
hist(mcols(y)$pvalue, col="grey")
```

As before, finding the genes with the most extreme log2 fold change:

```{r}
with(mcols(gy),
     table(sig=qvalue < .05, sign.lfc=sign(log2FC))
     )     
sig <- mcols(gy)$qvalue < .05
glo <- order(mcols(gy)$log2FC * sig)
ghi <- order(-mcols(gy)$log2FC * sig)
```

```{r}
gtop_up <- mcols(gy)[head(ghi),]
print(as.data.frame(gtop_up)[,cols], digits=3)
gtop_down <- mcols(gy)[head(glo),]
print(as.data.frame(gtop_down)[,cols], digits=3)
```

We can plot a particular one of these genes:

```{r}
plotInfReps(gy, idx=ghi[100], x="condition", cov="line")
```

As expected, the highly up-regulated genes are involved in immune
response. Many genes encoding guanylate-binding proteins (GBP) are
up-regulated, and these proteins are induced by interferon,
produced in response to infection by pathogenic microbes.

We can make an MA plot, where the genes in our FDR set are colored:

```{r}
plotMASwish(gy, alpha=.05)
```

Again, using the `addIds` function from *tximeta*, we can easily add
gene symbols to our gene-level expression analysis:

```{r}
library(org.Hs.eg.db)
gy <- addIds(gy, "SYMBOL", gene=TRUE)
```

We can then add gene symbols to our MA plot:

```{r}
plotMASwish(gy, alpha=.05, xlim=c(.5,5.5))
with(
  subset(mcols(gy), qvalue < .05 & abs(log2FC) > 3),
     text(log10mean, log2FC, SYMBOL,
          col="blue", pos=4, cex=.7)
)
```

# Differential transcript usage

We have added a new function `isoformProportions` which can be run
after `scaleInfReps` (and optionally after removing genes via 
`labelKeep` and subsetting the SummarizedExperiment). This function,
`isoformProportions` will create a new assay `isoProp` from the
scaledTPM counts, containing isoform proportions per gene. The same
procedure will also be applied to all the inferential replicates. Note
that after `isoformProportions` the transcripts from single isoform
genes will be removed, and the transcripts will be re-ordered by gene
(alphabetically by gene).

Following this function, running `swish` will be equivalent to a test
of differential transcript usage, taking account of the uncertainty in
transcript abundances, as it will look for transcripts where the
isoform proportions change across condition.

```{r}
# run on the transcript-level dataset
iso <- isoformProportions(y)
iso <- swish(iso, x="condition", pair="line")
```

# Interaction designs

We also provide in `swish` methods for testing if a condition effect
varies *across a secondary covariate*, using matched samples for
condition, or un-matched samples, which we refer to as "interactions"
in the software.

If matched samples are available, we compute the log2 fold change for
each pair of samples across condition in the same covariate group, and
then we use a Wilcoxon rank sum statistic for comparing the log2 fold
changes across the secondary covariate. For permutation significance,
the secondary covariate labels of the pairs are permuted. For
un-matched samples, multiple random "pseudo-pairs" of samples across
condition within the two covariate groups are chosen, and the
statistic computed as above, averaging over the random
pseudo-pairings. The motivation for the above permutation schemes is
to ensure the following condition, that "under the null hypothesis,
the likelihood of the data is invariant under these permutations"
[@anderson], where our null hypothesis specifically involves the
interaction between condition and the secondary covariate.

For the macrophage dataset we have been working with [@alasoo], we
have a 2x2 experimental design, with IFN gamma stimulation,
*Salmonella* infection, and both treatments, as well as control
samples. We have these four conditions across 6 cell lines from 6
donors (a subset of all the RNA-seq samples available). So we can use
the first method described above, where the cell line is used to match
samples across condition. 

## Condition and secondary covariates

We begin the interaction analysis by now loading in the full set of
samples, as before we only loaded the samples for the two group
analysis.
We then define two new factors indicating IFNg status and *Salmonella*
status.

```{r}
coldata <- read.csv(file.path(dir, "coldata.csv"))
coldata <- coldata[,c(1,2,3,5)]
names(coldata) <- c("names","id","line","condition")
coldata$files <- file.path(dir, "quants", coldata$names, "quant.sf.gz")
se <- tximeta(coldata)
```

```{r}
se$ifng <- factor(ifelse(
  grepl("IFNg",se$condition),
  "treated","control"))
se$salmonella <- factor(ifelse(
  grepl("SL1344",se$condition),
  "infected","control"))
with(colData(se),
     table(ifng, salmonella))
```

We will work with the chromosome 1 transcripts for demonstration:

```{r}
y2 <- se
y2 <- y2[seqnames(y2) == "chr1",]
```

## Create and check paired samples

Our implementation of the interaction design for matched samples takes
into account matched samples within the `x` condition, which we will
specify to be the *Salmonella* infection status. We will specify the
secondary covariate `cov` to be the IFN gamma treatment. We will look
for transcripts where the infection response changes based on IFN
gamma treatment.

To perform the analysis, we create a new variable `pair` which will
record which samples are related within groups defined by IFN gamma
treatment status.

```{r}
y2$pair <- factor(y2$line)
levels(y2$pair) <- LETTERS[1:6] # simplify names
y2$pair <- as.character(y2$pair)
y2$pair[y2$ifng == "control"]
y2$pair[y2$ifng == "treated"]
y2$pair[y2$ifng == "treated"] <- rep(LETTERS[7:12], each=2)
y2$pair <- factor(y2$pair)
head(table(y2$pair, y2$salmonella))
```

## Swish for interaction effects

We now perform `swish` analysis, specifying the *Salmonella* infection
as our main condition, the IFN gamma treatment as the secondary
covariate, and providing the pairing within IFN gamma treatment
groups. We specify `interaction=TRUE` to test for differences in
infection response across IFN gamma treatment group.

```{r results="hide", message=FALSE}
y2 <- scaleInfReps(y2)
y2 <- labelKeep(y2)
y2 <- y2[mcols(y2)$keep,]
set.seed(1)
y2 <- swish(y2, x="salmonella", cov="ifng", pair="pair", interaction=TRUE)
```

## Plotting interaction results

In this case, we appear to have fewer non-null p-values from first
impression of the p-value histogram:

```{r}
hist(mcols(y2)$pvalue, col="grey")
```

The MA plot shows significant transcripts on either side of
`log2FC=0`. Note that the log2 fold change reported is the
*difference* between the log2 fold change in the IFN gamma treated and
IFN gamma control group. So positive `log2FC` in this plot indicates
that the effect is higher with IGN gamma treatment than in absence of
the treatment.

```{r}
plotMASwish(y2, alpha=.05)
```

We can plot some of the transcripts with high log fold change
*difference* across IFN gamma treatment group, and which belong to the
less than 5% nominal FDR group:

```{r}
idx <- with(mcols(y2), which(qvalue < .05 & log2FC > 5))
plotInfReps(y2, idx[1], x="ifng", cov="salmonella")
plotInfReps(y2, idx[2], x="ifng", cov="salmonella")
```

# Allelic expression analysis

We have a separate vignette for allelic expression analysis, either
navigate to the 
[allelic vignette](https://thelovelab.github.io/fishpond/articles/allelic.html)
on the *fishpond* website, or the following from the R prompt:

```{r eval=FALSE}
vignette(topic="allelic", package="fishpond")
```

The allelic analysis methods are described in @wu2022 
[doi: 10.1101/2022.08.12.503785](https://doi.org/10.1101/2022.08.12.503785).

# Correlation test

*Swish* now has methods to compute correlations
(`"spearman"` or `"pearson"`) of a continuous variable `x` with the
log counts, and then assess the significance of those correlations
using inferential replicates and permutation tests.
Additionally, one can compute correlations of a covariate `cov` with
log fold changes across a pairing variable (where `x` labels the
condition for the log fold change and `pair` labels the pairs).
The Spearman-based test statistic follows the procedure from the
original method *SAMseq* that *Swish* is based upon [@samseq],
e.g. the methods used in that package with a quantitative response
type. See `?swish` and the `cor` argument for more details on the
use of correlations within *Swish*.

For examples of the correlation tests, see the dynamic AI examples in
the `"allelic"` fishpond vignette, or there are also some simulated
examples in the `test_correlation.R` testing script located in
fishpond's `tests/testthat` directory.

# alevin scRNA-seq

The *alevin* [@alevin] and *tximport* / *tximeta* maintainers have
created an efficient format for storing and importing the sparse 
scRNA-seq estimated gene counts, inferential mean and
variance data, and optionally the inferential replicate counts. `tximeta`
will automatically import these matrices if *alevin* was run using
`--numCellBootstraps` (in order to generate inferential variance) and
additionally `--dumpFeatures` (in order to dump the inferential
replicates themselves, see below on thoughts on avoiding this step
though). 

The storage format for counts, mean, variance, and inferential
replicates, involves writing one cell at a time, storing the locations
of the non-zero counts, and then the non-zero counts. The matrices are
loaded into R sparely using the *Matrix* package. The storage format is
efficient, for example, the estimated counts for the 900 mouse neuron
dataset from 10x Genomics takes up 4.2 Mb, the mean/variance matrices
take up 8.6 Mb each, and the inferential replicates takes up 72 Mb (20
bootstrap inferential replicates). Hence avoiding writing and
importing the inferential replicates themselves, and only using the
mean and variance, is desirable.

The `swish` and *alevin* authors have developed a workflow that uses
inferential mean and variance to generate pseudo-inferential
replicates in place of the actual inferential replicates. Storing and
importing only the inferential mean and variance is much more
efficient (dramatically faster load time and less space on disk and in
memory). The faster workflow would then skip `--dumpFeatures` when
running *alevin*, or subsequently use `dropInfReps=TRUE` to not load
the inferential replicates into R. 

**Plotting:** To demonstrate how the inferential mean and variance
look across real scRNA-seq data, we load the neurons dataset and
plot the inferential replicate data across cells. First we read in the
counts, in this case using `dropInfReps=TRUE` because the directory
includes the actual inferential replicates, not only the mean and
variance.  We set `skipMeta=TRUE`, although in general you would want
`tximeta` to add the gene range information and other metadata if
working with human, mouse, or fruit fly.

```{r eval=FALSE}
dir <- system.file("extdata", package="tximportData")
files <- file.path(dir,"alevin/neurons_900_v014/alevin/quants_mat.gz")
neurons <- tximeta(files, type="alevin",
                   skipMeta=TRUE, # just for vignette
                   dropInfReps=TRUE,
                   alevinArgs=list(filterBarcodes=TRUE))
```

We can easily make a *SingleCellExperiment* object [@Amezquita2020],
and plot counts across cluster (here a randomly assigned cluster
label). For more details on working with *SingleCellExperiment*
objects, consult the following online book: 
[Orchestrating Single-Cell Analysis with Bioconductor](https://osca.bioconductor.org/)
[@Amezquita2020].

```{r eval=FALSE}
library(SingleCellExperiment)
sce <- as(neurons, "SingleCellExperiment")
sce$cluster <- factor(paste0("cl",sample(1:6,ncol(sce),TRUE)))
par(mfrow=c(2,1), mar=c(2,4,2,1))
plotInfReps(sce, "ENSMUSG00000072235.6", x="cluster",
            legend=TRUE)
plotInfReps(sce, "ENSMUSG00000072235.6", x="cluster",
            reorder=FALSE)
```

`plotInfReps` has a number of options and convenience arguments for
visualizing single cell data. One can:

* add a `legend`,
* `reorder` the cells by expression value within cluster,
* apply size factors to the counts (`applySF`) (size factor scaling
will use the values in `sizeFactors(sce)` or equivalently
`mcols(sce)$sizeFactor`),
* use a column of `mcols(sce)` to set the `main` title,
  e.g. `mainCol="SYMBOL"`,
* specify `x` as a numeric covariate (e.g. pseudotime), and use `cov`
  to distinguish groups (e.g. lineage). Points will then be colored by
  `cov` instead of by discrete `x`.

See `?plotInfReps` for more description of arguments.

Note that the figures scale to some degree by the number of cells; for
example with $n \in [1,150)$ or $[150,400)$, more visual elements
per cell will be included: 

```{r echo=FALSE, eval=FALSE}
par(mfrow=c(2,1), mar=c(2,4,2,1))
plotInfReps(sce[,1:50], "ENSMUSG00000072235.6", x="cluster")
plotInfReps(sce[,1:150], "ENSMUSG00000072235.6", x="cluster")
```

**Advice on Swish testing:** `swish` can be run on *alevin* counts
imported with `tximeta`, but there are a few extra steps
suggested. The following does not use evaluated code chunks, but
provides suggestions for how to tailor `swish` to single-cell
datasets. 

**Filter genes:** we recommend to filter
genes as the first step, to reduce the size of the data before losing
sparsity on the count matrices (conversion of data to ranks loses data
sparsity inside the `swish()` function). One can run `labelKeep`
therefore before `scaleInfReps`. E.g., to remove genes for which there
are not 10 cells with a count of 3 or more:

```{r eval=FALSE}
y <- labelKeep(y, minCount=3, minN=10)
y <- y[mcols(y)$keep,] # subset genes
```

**Subset cells:** One should also subset to the groups of cells of
interest for differential testing, in order to take up the least
amount of memory when the sparse matrices are converted to dense
matrices. Note that `swish` only allows for differential testing of
two groups (although it allows for blocking factors and interaction
tests).

**(Slower) With inferential replicates:** After one has filtered both
genes and cells down to the set that are of interest for differential
testing, one can run the following commands, to make the sparse
matrices into dense ones, scale the cells, and perform *Swish*
differential expression, however read on for faster suggestions.

```{r eval=FALSE}
assays(y) <- lapply(assays(y), as.matrix) # make dense matrices
y <- scaleInfReps(y, lengthCorrect=FALSE, sfFun=sfFun)
y <- swish(y, x="condition")
```

**Size factor function:** Note that `scaleInfReps` has an argument
`sfFun` which allows the user to provide their own size factor
calculation function. We have found that `calculateSumFactors`
[@Lun2016] in the *scran* package [@scran] works well for computing
size factors. `sfFun` should be specified as a function that returns a
vector of size factors, or a numeric vector of the size factors
themselves. 

**(Faster) Without inferential replicates:** The following workflow
can be used in the case that `assayNames(y)` only contains counts,
mean, and variance, which is much faster by avoiding writing/importing
inferential replicates. We first generate pseudo-inferential
replicates from inferential mean and variance matrices before running
`scaleInfReps` from the code chunk above. The generation of
pseudo-inferential replicates is described in @compression.

The following function can be used just before `scaleInfReps`:

```{r, eval=FALSE}
y <- makeInfReps(y, 20)
```

There is alternatively a scheme for splitting the operation of
generating (dense) inferential replicate matrices across multiple
jobs, and running `swish` across batches of genes at a time. 
This job-splitting procedure is also described and benchmarked in
@compression.  This helps to reduce the total memory used, in the case
that the counts, mean, and variance matrices are too large to be made
dense altogether. This scheme involves 1) splitting the object into
smaller pieces, written out to `.rds` files, 2) running `swish` as a
distributed job, 3) reading the `.csv` output back into R. The
following code chunk would start with a 'y' with sparse matrix assays,
and without ever running `scaleInfReps` (it is run within the
distributed job, and with `lengthCorrect=FALSE` by default).

```{r, eval=FALSE}
library(SingleCellExperiment)
y <- as(y, "SingleCellExperiment")
# then, after filtering genes and cells...
# compute and store sizeFactors calculated over all genes
y <- scran::computeSumFactors(y)
# split swish objects into 8 RDS files:
splitSwish(y, nsplits=8, prefix="swish", snakefile="Snakefile")
# now, run snakemake from command line
# after finished, results back into R:
y <- addStatsFromCSV(y, "summary.csv")
```

The `splitSwish` function will write out a `Snakefile` that can be
used with [snakemake](https://snakemake.readthedocs.io/en/stable/) in
order to run distributed `swish` jobs in an easily customized
workflow. Then the `addStatsFromCSV` will read in and attach the
results to the original object. This final alternative avoids
generating dense matrices until they have been split into `nsplits`
pieces, and so can be used to reduce the memory requirements for the
individual jobs. If one is new to running `snakemake`, it is
recommended to first run with the flags `-np` as a "dry-run" to see
the operations that will be performed. The `swish` command can be
customized in the `swish` rule in the `Snakefile`, e.g. to control for
batches or test for interactions.

# Further details

## Analysis types supported by Swish

There are currently seven types of analysis supported by `swish`:

* Two group analysis
* Two groups with two or more batches
* Two group paired or matched samples
* Two condition x two group paired samples, interaction test
* Two condition x two group samples, not paired, interaction test
* Correlation test of expression with continuous x
* Correlation test of LFC for paired samples with continuous covariate

This vignette demonstrated the third in this list, but the others
can be run by either not specifying any additional covariates, or by
specifying a batch variable with the argument `cov` instead of `pair`.
The two interaction tests can be run by specifying `interaction=TRUE`
and providing `x`, `cov`, and optionally `pair`. The last two tests
can be run using the `cor` argument 
(see [Correlation test](#correlation-test) section in this vignette).

## Accounting for continuous variables

We have two recommended approaches for using `swish` in combination
with estimated batch effects, e.g. factors of unwanted variation
estimated by *RUVSeq* or surrogate variables estimated by *SVA*
(without loss of generality, we will call these batch factors).
First, examine a plot of the estimated batch factors, e.g. a stripplot
(1 factor) or scatterplot of pairs (2 or more factors). 

1. If the samples seem to fall into discrete clusters, one can run
   `kmeans` to assign a discrete estimated batch to each cluster. If
   the number of samples per condition group per batch is 2 or more,
   then one can use the `cluster` output of the `kmeans` function as
   input to the `cov` argument of `swish`. This will then perform a
   stratified Mann-Whitney Wilcoxon test.
2. If the samples do not fall into discrete clusters, one can use the
   same approach that we use to correct the samples for sequencing
   depth variation: direct scaling of the estimated counts across the
   inferential replicates. Why this works: usually we do not scale
   counts in statistical inference pipelines because we lose the
   information about precision of counts across the dynamic
   range. However, because we keep the set of inferential replicate
   matrices, and these contain both sampling and additional
   inferential variance, we are able to track how scaling the counts
   affects the variance, and this informs the test statistic.

For the second approach, one can directly scale the inferential counts
using *limma* and `removeBatchEffect`.

Here a demonstration using simulated data:

```{r}
set.seed(1)
y <- makeSimSwishData()
```

First we perform standard scaling for sequencing depth and labeling of
feautures with too few counts. We save the mean of scaled inferential
replicates using `saveMeanScaled=TRUE` as we will use this to estimate
the batch factors.

```{r}
y <- scaleInfReps(y, saveMeanScaled=TRUE)
tail(assayNames(y),4) # 'meanScaled' assay
y <- labelKeep(y)
y <- y[mcols(y)$keep,]
```

Next estimate the batch factors on the mean of scaled inferential
count matrices.

```{r}
norm_cts_for_batch <- assay(y, "meanScaled")
# use batch factor estimation method of your choosing
# ...
```

Suppose we estimate two batch factors, `w1` and `w2`, now we remove
the variation from each inferential replicate associated with these,
using limma's `removeBatchEffect`.

```{r}
w1 <- runif(ncol(y)) # here simulated w1, use real instead
w2 <- runif(ncol(y)) # here simulated w2, use real instead
W <- data.frame(w1, w2)
infRepIdx <- grep("infRep",assayNames(y),value=TRUE)
nreps <- length(infRepIdx)
```

The following loads *limma* and removes batch associated variation on
the log2 scale (a pseudocount is added to avoid `-Inf` values).

```{r}
suppressPackageStartupMessages(library(limma))
mm <- model.matrix(~condition, colData(y))
pc <- .1
for (k in seq_len(nreps)) {
  logInfRep <- log(assay(y, infRepIdx[k]) + pc)
  logInfRep <- limma::removeBatchEffect(
                        logInfRep,
                        covariates=W,
                        design=mm)
  assay(y, infRepIdx[k]) <- exp(logInfRep)
}
```

Now we run `swish` as usual:

```{r}
y <- swish(y, x="condition")
```

While we have assessed this approach with a small number of 
estimated nuisance covariates, note that with many covariates this
approach may lead to inflation of test statistics.

## Structure of tximeta output

While `tximeta` is the safest way to provide the correct input to
`swish`, all that `swish` requires for running is a
*SummarizedExperiment* object with the following assays: `counts`,
`length`, and `infRep1`, `infRep2`, ..., `infRepN`, where `N` is
simply the number of Gibbs samples or boostraps samples, e.g. 20 in
the examples above. The counts and inferential replicates are
estimated counts from a quantification method, either at the
transcript level or summed to the gene level (simple sum). These
counts sum up to the (mapped) library size for each sample. It is
assumed that the `length` matrix gives the effective lengths for each
transcript, or average transcript length for each gene as summarized
by the functions in `tximeta`/`tximport`. If the counts should not be
corrected for effective length (e.g. 3' tagged RNA-seq), then
`lengthCorrect=FALSE` should be specified when running
`scaleInfReps`. 

Note on simulation: it is difficult to simulate inferential uncertainty
in a realistic manner without construction of reads from transcripts,
using a method like *polyester*. Constructing reads from the reference
transcriptome or a sample-specific transcriptome naturally produces
the structure of read-assignment inferential uncertainty that `swish`
and other methods control for in real RNA-seq data.

## Plotting q-values over statistics

As with *SAMseq* and *SAM*, `swish` makes use of the permutation
plug-in approach for q-value calculation. `swish` calls the `empPvals`
and `qvalue` functions from the *qvalue* package to calculate the
q-values (or optionally similar functions from the *samr* package).
If we plot the q-values against the statistic, or against the log2
fold change, one can see clusters of genes with the same q-value
(because they have the same or similar statistic). One consequence of
this is that, in order to rank the genes, rather than ranking directly
by q-value, it makes more sense to pick a q-value threshold and then
within that set of genes, to rank by the log fold change, as shown
above when the code chunk has `log2FC * sig`.

```{r}
gres <- mcols(gy)[mcols(gy)$keep,]
min(gres$qvalue, na.rm=TRUE) # min nominal FDR is not 0
with(gres, plot(stat, -log10(qvalue)))
with(gres, plot(log2FC, -log10(qvalue)))
abline(v=0, col="red")
with(gres, plot(log2FC, -log10(qvalue),
                xlim=c(-1.5,1.5), ylim=c(0,1.5)))
abline(v=0, col="red")
```

## Plotting InfRV

In the Swish paper, we describe a statistic, InfRV, which is useful
for categorizing groups of features by their inferential uncertainty.
Note that InfRV is not used in the `swish` method, but only for
visualization in the paper. Here we show how to compute and plot the
InfRV:

```{r}
y3 <- se
y3 <- y3[seqnames(y3) == "chr4",]
y3 <- y3[,y3$condition %in% c("naive","IFNg")]
y3 <- labelKeep(y3)
y3 <- y3[mcols(y3)$keep,]
y3 <- computeInfRV(y3)
mcols(y3)$meanCts <- rowMeans(assays(y3)[["counts"]])
with(mcols(y3), plot(meanCts, meanInfRV, log="xy"))
hist(log10(mcols(y3)$meanInfRV),
     col="grey50", border="white", breaks=20,
     xlab="mean InfRV", main="Txp-level inferential uncertainty")
```

## Salmon in alignment mode, how to use tximeta

If you ran Salmon in alignment mode, then a Salmon index was not
used. Without a Salmon index, there is not a checksum that tximeta
uses to identify the transcriptome. You can specify `skipMeta=TRUE` to
avoid tximeta attempting to match on the transcriptome in this case. 

If you want to combine transcripts to gene level (or some other level
of aggregation), you can set `txOut=FALSE, tx2gene=tx2gene` when
calling `tximeta` and it will pass these arguments to tximport when
importing the counts and inferential replicates.

## Permutation schemes for interactions

The following diagrams describe the permutation schemes used for the
interaction designs implemented in `swish`. The case with matched
samples (pair indicated by number, primary condition indicated by
color, the vertical line separating the pairs by secondary covariate):

```{r echo=FALSE}
n <- 8
condition <- rep(1:2,length=2*n)
group <- rep(1:2,each=n)
pair <- rep(c(1:n),each=2)
cols <- c("dodgerblue","goldenrod4")
plot(1:(2*n), rep(0,2*n), ylim=c(-.5,3.5),
     type="n", xaxt="n", yaxt="n",
     xlab="samples", ylab="permutation")
abline(v=8.5, lty=2)
axis(2, 0:3, c("orig",1:3), las=2)
text(1:(2*n), rep(0,2*n), pair, col=cols[condition], cex=2)
set.seed(1)
for (i in 1:3) {
  perms <- rep(2*sample(n,n),each=2) - rep(1:0,length=2*n)
  text(1:(2*n), rep(i,2*n), pair[perms], col=cols[condition[perms]], cex=2)
}
```

The case without matched samples (sample indicated by letter, primary
condition indicated by color, the vertical line separating the samples
by secondary covariate). Here multiple random pseudo-pairs are chosen
across condition. The permutation scheme ensures that LFCs are always 
calculated between samples from the same covariate group.

```{r echo=FALSE}
n <- 8
condition <- rep(c(1:2,1:2),each=n/2)
group <- rep(1:2,each=n)
id <- LETTERS[1:(2*n)]
cols <- c("dodgerblue","goldenrod4")
plot(1:(2*n), rep(0,2*n), ylim=c(-.5,3.5),
     type="n", xaxt="n", yaxt="n",
     xlab="samples", ylab="permutation")
abline(v=8.5, lty=2)
axis(2, 0:3, c("orig",1:3), las=2)
text(1:(2*n), rep(0,2*n), id, col=cols[condition], cex=2)
set.seed(3)
for (i in 1:3) {
  id.perms <- character(2*n)
  grp1 <- id[group==1]
  grp2 <- id[group==2]
  id.perms[c(1:4,9:12)] <- sample(id[condition==1],n)
  idx1 <- id.perms[c(1:4,9:12)] %in% grp1
  id.perms[c(5:8,13:16)][idx1] <- sample(id[condition==2 & group==1],sum(idx1))
  idx2 <- id.perms[c(1:4,9:12)] %in% grp2
  id.perms[c(5:8,13:16)][idx2] <- sample(id[condition==2 & group==2],sum(idx2))
  text(1:(2*n), rep(i,2*n), id.perms, col=cols[condition], cex=2)
}
arrows(3,1.5,1.3,1.15,length=.1)
arrows(3,1.5,4.7,1.15,length=.1)
```

## Session information

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

# References