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.
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, namely the median_correlation
and outlier_fraction
functions.
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 can be installed from the MoseleyBioinformaticsLab r-universe as it is not yet on CRAN.
options(repos = c( moseleybioinformaticslab = 'https://moseleybioinformaticslab.r-universe.dev', BiocManager::repositories())) install.packages(c("ICIKendallTau", "visualizationQualityControl"))
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.
knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.width = 8, fig.height = 5, fig.path = "man/figures/")
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
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:
knitr::kable(visqc_score_contributions(pca_data$x))
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 package.
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.
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)
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")
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")
We can combine the median correlations and outlier fractions into a single score and then examine the distribution of scores to look for 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.
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 provides this correlation measure.
Lets add some missingness to our 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)))
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
cor_random_missing = ICIKendallTau::ici_kendalltau(exp_data2)$cor cor_random_missing[1:4, 1:4]
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?
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.
Some fake data is stored in grp_cor_data
that is useful for testing the median_correlation
function. It was generated by:
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)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.