README.Rmd

---
title: "Visualization Quality Control"
author: "Robert M Flight"
date: "`r Sys.time()`"
output: md_document
editor_options: 
  chunk_output_type: console
---

# Visualization Quality Control

A set of useful functions for calculating various measures from high-feature datasets and visualizing them.

In addition to internal documentation, the package is also documented heavily  [here](https://moseleybioinformaticslab.github.io/visualizationQualityControl/).

This package combines my needs for visualizing sample-sample correlations using heatmaps, and novel quality control measures that apply to different types of -omics or high-feature datasets proposed by [Gierlinski et al., 2015](https://dx.doi.org/10.1093/bioinformatics/btv425), namely the `median_correlation` and `outlier_fraction` functions.

## Installation

### Dependencies

* `ComplexHeatmap`, for generating the heatmaps.
* `dendsort`, for reordering samples in heatmaps.
* `ICIKendallTau`, for calculating Kendall-tau with missing values.

These should get installed automatically.

It is recommended to install `BiocManager` first so Bioconductor dependencies are installed automatically.

```
install.packages("BiocManager")
```

### This Package

This package can be installed from the MoseleyBioinformaticsLab r-universe as it is not yet on CRAN.

```r
options(repos = c(
    moseleybioinformaticslab = 'https://moseleybioinformaticslab.r-universe.dev',
    BiocManager::repositories()))
install.packages(c("ICIKendallTau", "visualizationQualityControl"))
```

## Examples

These examples show the primary functionality. We will apply the visualizations
to a two group dataset. However, all of the functions are still applicable to datasets with more than two groups. The examples below are for a dataset where there has been a sample swapped
between the two groups (i.e. there is a problem!). If you want to see how the 
visualizations compare between a **good** dataset and a **bad** dataset, see the **quality_control** vignette.

```{r setup, echo = FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.width = 8, fig.height = 5,
                      fig.path = "man/figures/")
```


```{r example_data}
library(visualizationQualityControl)
library(ggplot2)
library(ggforce)
data("grp_cor_data")
exp_data = grp_cor_data$data
rownames(exp_data) = paste0("f", seq(1, nrow(exp_data)))
colnames(exp_data) = paste0("s", seq(1, ncol(exp_data)))

sample_info = data.frame(id = colnames(exp_data), class = grp_cor_data$class)

exp_data[, 5] = grp_cor_data$data[, 19]
exp_data[, 19] = grp_cor_data$data[, 5]
sample_classes = sample_info$class
```


### Visualize PCA Component Scores

```{r do_pca}
pca_data = prcomp(t(exp_data), center = TRUE)
pca_scores = as.data.frame(pca_data$x)
pca_scores = cbind(pca_scores, sample_info)
ggplot(pca_scores, aes(x = PC1, y = PC2, color = class)) + geom_point()
```

To see how much explained variance each PC has, you can calculate them:

```{r show_variance}
knitr::kable(visqc_score_contributions(pca_data$x))
```

### visqc_heatmap

Calculate sample-sample correlations and reorder based on within class correlations.
We recommend a transform agnostic correlation like Kendall-tau that can also handle missing data when necessary.
Here we use the {ici_kendalltau} function from our [ICIKendallTau](https://moseleybioinformaticslab.github.io/ICIKendallTau/) package.

```{r correlation}
rownames(sample_info) = sample_info$id
data_cor = ICIKendallTau::ici_kendalltau(exp_data)
data_order = similarity_reorderbyclass(data_cor$cor, sample_info[, "class", drop = FALSE], transform = "sub_1")
```

And then generate a colormapping for the sample classes and plot the correlation heatmap.

```{r colormapandheatmap}
data_legend = generate_group_colors(2)
names(data_legend) = c("grp1", "grp2")
row_data = sample_info[, "class", drop = FALSE]
row_annotation = list(class = data_legend)

library(viridis)
suppressPackageStartupMessages(library(circlize))
colormap = colorRamp2(seq(0.3, 1, length.out = 50), viridis::viridis(50))

visqc_heatmap(data_cor$cor, colormap, "Correlation", row_color_data = row_data,
              row_color_list = row_annotation, col_color_data = row_data,
              col_color_list = row_annotation, row_order = data_order$indices,
              column_order = data_order$indices)
```

### median_correlations

```{r median_correlations}
data_medcor = median_correlations(data_cor$cor, sample_info$class)
ggplot(data_medcor, aes(x = sample_id, y = med_cor)) + geom_point() + 
  facet_grid(. ~ sample_class, scales = "free_x") + ggtitle("Median Correlation")

ggplot(data_medcor, aes(x = sample_class, y = med_cor)) +
  geom_sina() +
  ggtitle("Median Correlation")
```

### outlier_fraction

```{r outlier_fraction}
data_outlier = outlier_fraction(exp_data, sample_info$class)
ggplot(data_outlier, aes(x = sample_id, y = frac)) + geom_point() + 
  facet_grid(. ~ sample_class, scales = "free_x") + ggtitle("Outlier Fraction")

ggplot(data_outlier, aes(x = sample_class, y = frac)) +
  geom_sina() +
  ggtitle("Outlier Fraction")
```

### determine_outliers

We can combine the median correlations and outlier fractions into a single score and then examine the distribution of scores to look for outliers.

```{r determine_outliers}
out_samples = determine_outliers(data_medcor, data_outlier)

ggplot(out_samples, aes(x = sample_id, y = score, color = outlier)) +
  geom_point() +
  facet_wrap(~ sample_class, scales = "free_x") +
  ggtitle("Outliers Score")

ggplot(out_samples, aes(x = sample_class, y = score, color = outlier, group = sample_class)) +
  geom_sina() +
  ggtitle("Outliers Score")
```

Here we can see the outliers by their combined score.
**However**, in this case we don't actually want to remove the samples.
In this example, what actually happened was that two samples got their `sample_class` wrong.
And we can see that by going back to the **correlation heatmap**, that this is the case by the high correlation values observed with the other class of samples.

### Correlation that Includes Missing Values

When there are missing values (either NA, or 0 depending on the case), we can use the information-content-informed Kendall-tau.
This works under the assumption that **most** missing data in -omics is because samples have values that fall below the detection limit. 
Because of this, missingness actually contributes **some** information that can be incorporated in the correlation.
The package [ICIKendallTau](https://moseleybioinformaticslab.github.io/ICIKendallTau/) provides this correlation measure.

Lets add some missingness to our data.

```{r}
exp_data = grp_cor_data$data
rownames(exp_data) = paste0("f", seq(1, nrow(exp_data)))
colnames(exp_data) = paste0("s", seq(1, ncol(exp_data)))
```

```{r add_random_missingness}
make_na = rep(FALSE, nrow(exp_data))
s1_missing = make_na
s1_missing[sample(length(make_na), 20)] = TRUE
s2_missing = make_na
s2_missing[sample(which(!s1_missing), 20)] = TRUE

exp_data2 = exp_data
exp_data2[s1_missing, 1] = NA
exp_data2[s2_missing, 1] = NA
```

```{r cor_random_missingness}
cor_random_missing = ICIKendallTau::ici_kendalltau(exp_data2)$cor
cor_random_missing[1:4, 1:4]
```

```{r cor_random_missingness_noweight}
cor_random_missing_nw = ICIKendallTau::ici_kendalltau(exp_data)$cor
cor_random_missing_nw[1:4, 1:4]
```

What happens if we make the missingness match between them? That counts as information? 
If the feature is missing in the same samples, that is
worth something?

```{r cor_same_missingness}
exp_data = grp_cor_data$data
rownames(exp_data) = paste0("f", seq(1, nrow(exp_data)))
colnames(exp_data) = paste0("s", seq(1, ncol(exp_data)))
exp_data[s1_missing, 1:2] = NA

cor_same_missing = ICIKendallTau::ici_kendalltau(exp_data)$cor
cor_same_missing[1:4, 1:4]
```

Here we can see that the correlation between sapmles S1 and S2 has actually increased over the random missing case.

## Fake Data Generation

Some fake data is stored in `grp_cor_data` that is useful for testing the `median_correlation`
function. It was generated by:

```{r fakedata, eval=FALSE}
library(fakeDataWithError)
set.seed(1234)

s1 = runif(100, 0, 1)
grp1 = add_uniform_noise(10, s1, 0.1)

model_data = data.frame(s1 = s1, s2 = grp1[, 1])

lm_1 = lm(s1 ~ s2, data = model_data)

lm_1$coefficients[2] = 0.5

s3 = predict(lm_1)
s4 = add_uniform_noise(1, s3, 0.2)

grp2 = add_uniform_noise(10, s4, 0.1)

grp_class = rep(c("grp1", "grp2"), each = 10)

grp_cor_data = list(data = cbind(grp1, grp2), class = grp_class)
```

```{r fakedata2, eval=FALSE}
library(fakeDataWithError)
set.seed(1234)

n_point = 1000
n_rep = 10

# a nice log-normal distribution of points with points along the entire range
simulated_data = c(rlnorm(n_point / 2, meanlog = 1, sdlog = 1),
                    runif(n_point / 2, 5, 100))

# go to log to have decent correlations on the "transformed" data
lsim1 = log(simulated_data)

# add some uniform noise to get lower than 1 correlations
lgrp1 = add_uniform_noise(n_rep, lsim1, .5)

# add some uniform noise to everything in normal space
sim1_error = add_uniform_noise(n_rep, simulated_data, 1, use_zero = TRUE)
# and generate the grp1 data in normal space
ngrp1 = exp(lgrp1) + sim1_error


# do regression to generate some other data
model_data = data.frame(lsim1 = lsim1, lsim2 = lgrp1[, 1])
lm_1 = lm(lsim1 ~ lsim2, data = model_data)

# reduce the correlation between them
lm_1$coefficients[2] = 0.5
lsim3 = predict(lm_1)

# and a bunch of error
lsim4 = add_uniform_noise(1, lsim3, 1.5)

# create group with added error to reduce correlation from 1
lgrp2 = add_uniform_noise(10, lsim4, .5)

# add error in original space
nsim4 = exp(lsim4)
sim4_error = add_uniform_noise(10, nsim4, 1, use_zero = TRUE)
ngrp2 = exp(lgrp2) + sim4_error

# put all data together, and make negatives zero
all_data = cbind(ngrp1, ngrp2)
all_data[(all_data < 0)] = 0

grp_class = rep(c("grp1", "grp2"), each = 10)

grp_exp_data = list(data = all_data, class = grp_class)
```