13-exprs-norm-reads.Rmd

---
knit: bookdown::preview_chapter
---

# Normalization for library size (Reads)

```{r, echo=FALSE, message=FALSE, warning=FALSE}
library(scRNA.seq.funcs)
library(scater, quietly = TRUE)
options(stringsAsFactors = FALSE)
library(knitr)
opts_chunk$set(out.width='90%', fig.align = 'center', echo=FALSE)
reads <- readRDS("blischak/reads.rds")
reads.qc <- reads[fData(reads)$use, pData(reads)$use]
endog_genes <- !fData(reads.qc)$is_feature_control
```

```{r norm-pca-raw-reads, fig.cap = "PCA plot of the blischak data"}
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "counts")
```
```{r norm-ours-rle-raw-reads, fig.cap = "Cell-wise RLE of the blischak data"}
boxplot(calc_cell_RLE(counts(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim=c(-1,1))
```

```{r norm-pca-cpm-reads, fig.cap = "PCA plot of the blischak data after CPM normalisation"}
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "exprs")
```
```{r norm-ours-rle-cpm-reads, fig.cap = "Cell-wise RLE of the blischak data", warning=FALSE}
boxplot(calc_cell_RLE(exprs(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim = c(-1,1))
```

```{r norm-pca-tmm-reads, fig.cap = "PCA plot of the blischak data after TMM normalisation"}
reads.qc <- 
    scater::normaliseExprs(reads.qc,
                           method = "TMM",
                           feature_set = endog_genes)
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "norm_counts")
```
```{r norm-ours-rle-tmm-reads, fig.cap = "Cell-wise RLE of the blischak data"}
boxplot(calc_cell_RLE(norm_counts(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim=c(-1,1))
```

```{r norm-pca-lsf-umi, fig.cap = "PCA plot of the blischak data after LSF normalisation"}
qclust <- scran::quickCluster(reads.qc, min.size = 30)
reads.qc <- scran::computeSumFactors(reads.qc, sizes = 15, clusters = qclust)
reads.qc <- scater::normalize(reads.qc)
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "exprs")
```

```{r norm-ours-rle-scran-reads, fig.cap = "Cell-wise RLE of the blischak data"}
boxplot(calc_cell_RLE(exprs(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim=c(-1,1))
```

```{r norm-pca-rle-reads, fig.cap = "PCA plot of the blischak data after RLE normalisation"}
reads.qc <- 
    scater::normaliseExprs(reads.qc,
                           method = "RLE",
                           feature_set = endog_genes)
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "norm_counts")
```

```{r norm-ours-rle-rle-reads, fig.cap = "Cell-wise RLE of the blischak data"}
boxplot(calc_cell_RLE(norm_counts(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim=c(-1,1))
```

```{r norm-pca-uq-reads, fig.cap = "PCA plot of the blischak data after UQ normalisation"}
reads.qc <- 
    scater::normaliseExprs(reads.qc,
                           method = "upperquartile", 
                           feature_set = endog_genes,
                           p = 0.99)
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "norm_counts")
```
```{r norm-ours-rle-uq-reads, fig.cap = "Cell-wise RLE of the blischak data"}
boxplot(calc_cell_RLE(norm_counts(reads.qc[endog_genes, ])),
        col = "grey50",
        ylab = "RLE",
        main = "", ylim=c(-1,1))
```

```{r norm-pca-downsample-reads, fig.cap = "PCA plot of the blischak data after downsampling"}
norm_counts(reads.qc) <- 
    scRNA.seq.funcs::Down_Sample_Matrix(counts(reads.qc))
scater::plotPCA(reads.qc[endog_genes, ],
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "norm_counts")
```
```{r norm-ours-rle-downsample-reads, fig.cap = "Cell-wise RLE of the blischak data"}
tmp <- norm_counts(reads.qc[endog_genes, ])
# ignore genes which are not detected in any cells following downsampling
boxplot(calc_cell_RLE(tmp[rowMeans(tmp) > 0, ]), 
        col = "grey50",
        ylab = "RLE",
        main = "", ylim = c(-1, 1))
```

```{r}
reads.qc <-
    scater::getBMFeatureAnnos(reads.qc,
                              filters = "ensembl_gene_id", 
                              attributes = c("ensembl_gene_id",
                                             "hgnc_symbol",
                                             "chromosome_name",
                                             "start_position",
                                             "end_position"), 
                              feature_symbol = "hgnc_symbol",
                              feature_id = "ensembl_gene_id",
                              biomart = "ENSEMBL_MART_ENSEMBL", 
                              dataset = "hsapiens_gene_ensembl",
                              host = "www.ensembl.org")
```

```{r}
reads.qc.ann <-
    reads.qc[!is.na(fData(reads.qc)$ensembl_gene_id), ]
```

```{r}
eff_length <- abs(fData(reads.qc.ann)$end_position -
                      fData(reads.qc.ann)$start_position)/1000
```

```{r}
tpm(reads.qc.ann) <-
    calculateTPM(
        reads.qc.ann,
        eff_length
    )
fpkm(reads.qc.ann) <-
    calculateFPKM(
        reads.qc.ann,
        eff_length
    )
```

```{r norm-pca-fpkm-reads, fig.cap = "PCA plot of the blischak data after FPKM normalisation"}
scater::plotPCA(reads.qc.ann,
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "fpkm")
```

```{r norm-pca-tpm-reads, fig.cap = "PCA plot of the blischak data after TPM normalisation"}
scater::plotPCA(reads.qc.ann,
                colour_by = "batch",
                size_by = "total_features",
                shape_by = "individual",
                exprs_values = "tpm")
```