---
title: "Spatial quantile normalization for co-expression analysis"
shorttitle: "spqn guide"
author:
  - Yi Wang
  - Kasper D. Hansen
package: spqn
abstract: >
  A guide to using spqn for co-expression analysis
bibliography: spqn.bib
vignette: >
  %\VignetteIndexEntry{spqn User's Guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}  
output:
  BiocStyle::html_document:
    toc_float: true
---

# Introduction

The spqn package contains an implementation of spatial quantile normalization (SpQN) as described in [@spqnPreprint]. In this work we describe how correlation matrices estimated from gene expression data exhibits a _mean-correlation_ relationship, in which highly expressed genes appear more correlated than lowly expressed genes. This relationship is not observed in protein-protein interaction networks, which suggests it is technical rather than biological. To correct for this bias in estimation, we developed SpQN, a method for normalizing the correlation matrix. In our work, we show how correcting for this relationship increases the correlations involving transcription factors, an important class of gene regulators.

We believe spatial quantile normalization has applications outside of co-expression analysis including beyond correlation matrices. However, the _implementation_ in this package is currently focused on correlation matrices.

The package contains various plotting functions for assessing the mean-correlation relationship, mirroring plots in [@spqnPreprint].


# Preliminaries

Let's load the package and some example data.

```{r load, message=FALSE}
library(spqn)
library(spqnData)
library(SummarizedExperiment)
```

This dataset contains 4,000 genes from the GTEx tissue "Adipose Subcutaneous" which we use as exemplar tissue in our preprint [@spqnPreprint]. The data is from GTEx v6p. The genes have been restricted to a random sample of 4,000 genes from genes which are

- either protein-coding or lincRNAs
- expressed (median expression greater than 0 on the $\log_2(\text{RPKM})$ scale

The data is counts (from the GTEx website) of bulk RNA-seq data which has been transformed to the expression matrix in the $\log_2(\text{RPKM})$ scale. Each gene has been scaled to have mean 0 and variance 1, and 4 principal components have been removed from the data matrix. For SpQN we need access to the gene expression level, which we take as the mean or median expression level across samples. So this quantity needs to be computed (and stored) prior to removing principal components. For this example data, we have stored the average expression value on the log-RPKM scale in `rowData(gtex.4k)$ave_logrpkm`.

To compare the distribution of correlation across different expression levels, we get the correlation matrix and average gene expression level ($\log_2(\text{RPKM})$) for each gene.
```{r cor_m, message=FALSE}
data(gtex.4k)
cor_m <- cor(t(assay(gtex.4k)))
ave_logrpkm <- rowData(gtex.4k)$ave_logrpkm
```

# Preparing your own data

**Note**: code in this section is not being evaluated.

We assume you have counts (or something similar) stored in a `SummarizedExperiment`. In our case, we have downloaded the GTEx v6p counts from the GTEx portal website (code is available in the `spqnData` package inside `scripts`). 

It makes sense to do some correction for library size. In this data we have simply divided by the number of reads in millions. There are multiple normalization strategies for estimating a library size correction factor. The literature shows that using such methods are *essential* when comparing expression across different tissues / cell types / conditions. However, here we are looking at correlation patterns **within** tissue.

We also need to remove genes which are not expressed. The threshold we have been using for GTEx is fairly liberal as it removes less than 50% of protein-coding codes. We apply this cutoff on the log2RPKM scale as it is necessary to consider gene length when applying a cutoff across genes. However, the exact way of doing this is probably not too important.

Finally, following [@Parsana2019] we strongly recommend removing principal components from the data matrix prior to analysis. A fast convenience function for this is `removePrincipalComponents()` from the WGCNA package. How many principal components to remove is an open question. In [@Freytag2015] the authors recommend looking at control genes to assess this, under the hypothesis that random genes should be uncorrelated and we have existing groups of genes which are believe to be highly correlated in all tissues, such as the PRC2 complex [@Freytag2015] or the ribosomal RNA genes [@Boukas2019]. Another approach is to use the `num.sv()` function from the SVA package. In our experience, this yields a substantially higher number of PCs than the graphical methods. Finally, as we shown [@spqnPreprint], the QC plots we present below might also help making this decision. At the end of the day, it is still an open problem.

Example code is as follows, the starting point is a `SummarizedExperiment` called `gtex` with the relevant data. It has two relevant columns in `rowData(gtex)`. One column is `gene_length` and another is `gene_type`.

```{r, examplePreproc, eval=FALSE}
# only keep coding genes and lincRNA in the counts matrix
gtex <- gtex[rowData(gtex)$gene_type %in%
                          c("lincRNA", "protein_coding"), ]

# normalize the counts matrix by log2rpkm
cSums <- colSums(assay(gtex))
logrpkm <- sweep(log2(assay(gtex) + 0.5), 2, FUN = "-",
                 STATS = log2(cSums / 10^6))
logrpkm <- logrpkm - log2(rowData(gtex)$gene_length / 1000)

# only keep those genes with median log2rpkm above 0
wh.expressed  <- which(rowMedians(logrpkm) > 0)

gtex.0pcs <- gtex[wh.expressed,]
logrpkm.0pcs <- logrpkm[wh.expressed,]
ave_logrpkm <- rowMeans(logrpkm.0pcs)
logrpkm <- logrpkm - ave_logrpkm # mean centering
logrpkm  <- logrpkm / matrixStats::rowSds(logrpkm) # variance scaling
assays(gtex.0pcs) <- SimpleList(logrpkm = logrpkm.0pcs)
rowData(gtex.0pcs)$ave_logrpkm <- ave_logrpkm

# remove PCs from the gene expression matrix after scaling each gene to have mean=0 and variance=1
gtex.4pcs <- gtex.0pcs
assay(gtex.4pcs) <- removePrincipalComponents(t(scale(t(logrpkm.0pcs))), n = 4)
```


# Examine the mean-correlation relationship

We now recreate a number of diagnostic plots which illustrate the mean-correlation relationship. These plots are similar to that is being presented in [@spqnPreprint].

To compare the distribution of correlation across different expression levels, we re-arrrange the correlation matrix by sortng the row and column by the expression level, and partitioned the correlation matrix into same-size bins, then used each bin to estimate the distribution of correlation for genes with that expression level. Here we use 10 by 10 bins. Each bin corresponds to two groups of genes that correspodnd to two expression levels. 

We first looked into the bins on the diagnal, where the two groups of genes have the same expression level. We looked at the distribution of the correlations within each of these bins. We find that the correlations of highly expressed genes are more variant than those of lowly expressed genes. We call such relationship between gene expression level and the distribution of correlation as mean-correlation relationship.
```{r, message=FALSE}
plot_signal_condition_exp(cor_m, ave_logrpkm, signal=0)
```


Next, we examine the impact of this mean-correlation relationship on the downstream analysis. 
A number of studies predict the co-expressed genes by setting a threshold. We therefore examine how the mean-correlation relationship would bias such prediction. We set 0.1% highest correlations as predicted signals in each bin. We compare the distibution of signals for bins on the diagonal. We find that the the predicted signals are stronger for high expressed genes and weaker for low expressed genes. Because of such bias, the true signals among low expressed genes could be hidden and may not be detected. 
```{r, message = FALSE}
plot_signal_condition_exp(cor_m, ave_logrpkm, signal=0.001)
```

Next we looked further into the bins across the entire correlation matrix. We quantify the variance inside a bin using inter-quantile range (IQR), which is the distance from the 1st and 3rd quantile. We use the width of box in the plot to denote IQR. We observe that the IQR would increase with the increase of expression level for either gene.
```{r, message=FALSE}
IQR_list <- get_IQR_condition_exp(cor_m, rowData(gtex.4k)$ave_logrpkm)
plot_IQR_condition_exp(IQR_list)
```

We next plot the marginal relationship between IQR of correlation and the expression level of the lower expression group for each bin, and we also see a positive relationship.
```{r, message = FALSE}
IQR_unlist <- unlist(lapply(1:10, function(ii) IQR_list$IQR_cor_mat[ii, ii:10]))
plot(rep(IQR_list$grp_mean, times = 1:10),
     IQR_unlist,
     xlab="min(average(log2RPKM))", ylab="IQR", cex.lab=1.5, cex.axis=1.2, col="blue")
```

We further examined the difference of the distributions using Q-Q plot. We used the bins with the second highest expression level on the diagonal as reference, and compared the distribution with other bins in the diagonal. We observe both scale and shape difference in the distributions.
```{r, message = FALSE}
par(mfrow = c(3,3))
for(j in c(1:8,10)){
    qqplot_condition_exp(cor_m, ave_logrpkm, j, j)
}
```

# Using SpQN to remove mean-correlation relationship
To remove the mean-correlation relationship and avoid the bias in the co-expression analysis, we develop a normalization method for correlation matrix, named as spatial quantile normalization (SpQN). The main idea is to use bins to partition the correlation matrix with genes sorted by exp level, and then apply quantile normalization for each bin. To smooth the normalization, we used a larger bin to approximate the empirical distribution of the inner bin. And we applied 1-dimension quantile normalization on the inner bins using its approximated distribution. The target distribution we used here is the bin with the second highest expression level on the diagonal. 

Using SpQN is relatively easy. In the function *normalize_correlation*, imput matrix *cor_m* is the correlation matrix. Imput vector *ave_exp* is the average expression level for the genes for the normalized expression matrix, corresponding to the row/column of *cor_m*. For other types of data, *ave_exp* can be the vector corresponding to the row/column of the correlation matrix, whose dependency with the distribution of correlations need to be removed. Parameter *ngrp* is an integer, denoting the number of bins in each row/column to be used to partition the correlation matrix in the process of SpQN. For example, setting *ngrp*=20 means to partition the correlation matrix into 20 by 20 bins. Parameter *size_grp* is an integer, denoting the size of the outter bins to be used to appriximate the distribution of the inner bins, in order to smooth the normalization. Note that the product of *size_grp* and *ngrp* must be equal or larger than than the row/column number of *cor_m*, and there is no smoothness in the normalization when they are equal. Parameter *ref_grp* is the location of the reference bin on the diagonal, whose distribution will be used as target distribution in the normalization.
```{r spqn, message = FALSE}
cor_m_spqn <- normalize_correlation(cor_m, ave_exp=ave_logrpkm, ngrp=20, size_grp=300, ref_grp=18)
```

# Assess the impact of normalization

To evaluate the performance of SpQN, we make the same plots as before, but using the normalized correlation matrix.

Compare distributions of correlations for bins on the diagonal of the correlation matrix.
```{r, message = FALSE}
plot_signal_condition_exp(cor_m_spqn, ave_logrpkm, signal=0)
```

Distribution of predicted signals and background correlations for bins on the diagonal of the correlation matrix.
```{r, message = FALSE}
plot_signal_condition_exp(cor_m_spqn, ave_logrpkm, signal=0.001)
```

Compare IQRs of bins across the entire correlation matrix.
```{r, message=FALSE}
IQR_spqn_list <- get_IQR_condition_exp(cor_m_spqn, rowData(gtex.4k)$ave_logrpkm)
plot_IQR_condition_exp(IQR_spqn_list)
```

Q-Q plot to compare the distribution for the bin with the 2nd highest expression level on the diagonal with other bins on the diagonal.
```{r, message = FALSE}
par(mfrow = c(3,3))
for(j in c(1:8,10)){
    qqplot_condition_exp(cor_m_spqn, ave_logrpkm, j, j)
}
```

Marginal relationship between IQR and the expression level of the lower expressed group for each bin.
```{r, message = FALSE}
IQR_unlist <- unlist(lapply(1:10, function(ii) IQR_spqn_list$IQR_cor_mat[ii, ii:10]))
plot(rep(IQR_spqn_list$grp_mean, times = 1:10),
     IQR_unlist,
     xlab="min(average(log2RPKM))", ylab="IQR", cex.lab=1.5, cex.axis=1.2, col="blue")
```

# SessionInfo

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

# References
