In this tutorial, we will construct a Symphony reference from two PBMC datasets from 2 technologies (10x 3’v1 and 3’v2), then map a third dataset from a new technology (10x 5’) with Symphony. The analysis follows from Fig. 2 of the paper (but with downsampled datasets to fit within CRAN limits on subdirectory size).
Install Symphony with standard commands.
install.packages('symphony')
Once Symphony is installed, load it up!
library(symphony)
# Other packages for this tutorial
suppressPackageStartupMessages({
# Analysis
library(harmony)
library(irlba)
library(data.table)
library(dplyr)
# Plotting
library(ggplot2)
library(ggthemes)
library(ggrastr)
library(RColorBrewer)
})
Get the expression and metadata.
load('../data/pbmcs_exprs_small.rda')
load('../data/pbmcs_meta_small.rda')
dim(pbmcs_exprs_small)
## [1] 1767 1200
dim(pbmcs_meta_small)
## [1] 1200 7
%>% head(4) pbmcs_meta_small
## cell_id donor nUMI nGene percent_mito cell_type
## 17642 fivePrime_GATCGCGAGACCGGAT 5p 6580 2020 0.05714286 B
## 4820 threepfresh_GATCGTACAAAGGTGC 3pv2 4426 1207 0.03366471 T_CD8
## 13347 fivePrime_AACTCAGGTTACGCGC 5p 4750 1670 0.05642105 B
## 10634 threepv1_CGTGATGAGTCCTC 3pv1 4121 1228 0.01504489 T_CD4
## cell_type_broad
## 17642 B
## 4820 T
## 13347 B
## 10634 T
Subset the dataset into reference and query.
= which(pbmcs_meta_small$donor == "5p") # use 5' dataset as the query
idx_query = pbmcs_exprs_small[, -idx_query]
ref_exp_full = pbmcs_meta_small[-idx_query, ]
ref_metadata = pbmcs_exprs_small[, idx_query]
query_exp = pbmcs_meta_small[idx_query, ] query_metadata
There are two options for how to build a Symphony reference. Option 1 (buildReferenceFromHarmonyObj
) is the more modular option, meaning that the user has more control over the preprocessing steps prior to reference compression. Option 2 (buildReference
) builds a reference starting from expression, automating the procedure more but offering less flexibility.
We’ll demonstrate both options below.
This option consists of more steps than Option 2 but allows your code to be more modular and flexible if you want to do your own preprocessing steps before the Harmony integration step. We recommend this option for most users.
It is important to generate vargenes_means_sds
(containing variable gene means and standard deviations used to scale the genes) as well as save the loadings for the PCA step.
Starting with the reference expression,
1:5, 1:2] # Sparse matrix with the normalized genes x cells matrix ref_exp_full[
## 5 x 2 sparse Matrix of class "dgCMatrix"
## threepfresh_GATCGTACAAAGGTGC threepv1_CGTGATGAGTCCTC
## LYZ 1.181536 1.766987
## FTL 2.678021 2.370840
## HLA-DRA 1.708152 .
## CD74 1.181536 1.766987
## S100A4 1.708152 2.113817
Select variable genes and subset reference expression by variable genes (the command below will select the top 1,000 genes per batch, then pool them)
= vargenes_vst(ref_exp_full, groups = as.character(ref_metadata[['donor']]), topn = 1000)
var_genes = ref_exp_full[var_genes, ]
ref_exp dim(ref_exp)
## [1] 1600 770
Calculate and save the mean and standard deviations for each gene
= tibble(symbol = var_genes, mean = Matrix::rowMeans(ref_exp))
vargenes_means_sds $stddev = rowSDs(ref_exp, vargenes_means_sds$mean)
vargenes_means_sdshead(vargenes_means_sds)
## # A tibble: 6 x 3
## symbol mean stddev
## <chr> <dbl> <dbl>
## 1 LYZ 1.78 1.91
## 2 HLA-DRA 1.87 1.69
## 3 FTL 3.55 1.22
## 4 CD74 2.54 1.53
## 5 S100A9 1.62 1.83
## 6 S100A4 2.53 1.46
Scale data using calculated gene means and standard deviations
= scaleDataWithStats(ref_exp, vargenes_means_sds$mean, vargenes_means_sds$stddev, 1) ref_exp_scaled
Run PCA (using SVD), save gene loadings (s$u
)
set.seed(0)
= irlba(ref_exp_scaled, nv = 20)
s = diag(s$d) %*% t(s$v) # [pcs by cells]
Z_pca_ref = s$u loadings
Run Harmony integration. It is important to set return_object = TRUE
.
set.seed(0)
= harmony::HarmonyMatrix(
ref_harmObj data_mat = t(Z_pca_ref), ## PCA embedding matrix of cells
meta_data = ref_metadata, ## dataframe with cell labels
theta = c(2), ## cluster diversity enforcement
vars_use = c('donor'), ## variable to integrate out
nclust = 100, ## number of clusters in Harmony model
max.iter.harmony = 20, ## max number of iterations
return_object = TRUE, ## return the full Harmony model object
do_pca = FALSE ## don't recompute PCs
)
## Harmony 1/20
## Harmony 2/20
## Harmony 3/20
## Harmony 4/20
## Harmony 5/20
## Harmony 6/20
## Harmony 7/20
## Harmony 8/20
## Harmony 9/20
## Harmony converged after 9 iterations
To run the next function buildReferenceFromHarmonyObj()
, you need to input the saved gene loadings (loadings
) and vargenes_means_sds
.
# Compress a Harmony object into a Symphony reference
= buildReferenceFromHarmonyObj(
reference # output object from HarmonyMatrix()
ref_harmObj, # reference cell metadata
ref_metadata, # gene names, means, and std devs for scaling
vargenes_means_sds, # genes x PCs matrix
loadings, verbose = TRUE, # verbose output
do_umap = TRUE, # set to TRUE to run UMAP
save_uwot_path = './testing_uwot_model_1') # file path to save uwot model
## Save metadata, vargenes (S), and loadings (U)
## Save R, Z_orig, Z_corr, and betas from Harmony object
## Calculate final L2 normalized reference centroids (Y_cos)
## Calculate reference compression terms (Nr and C)
## UMAP
## File already exists at that path... overwriting...
## Warning: invalid uid value replaced by that for user 'nobody'
## Saved uwot model
## Finished nicely.
Save Symphony reference for later mapping (modify with your desired output path)
saveRDS(reference, './testing_reference1.rds')
Let’s take a look at what the reference object contains: * meta_data: metadata * vargenes: variable genes, means, and standard deviations used for scaling * loadings: gene loadings for projection into pre-Harmony PC space * R: Soft cluster assignments * Z_orig: Pre-Harmony PC embedding * Z_corr: Harmonized PC embedding * centroids: locations of final Harmony soft cluster centroids * cache: pre-calculated reference-dependent portions of the mixture model * umap: UMAP coordinates * save_uwot_path: path to saved uwot model (for query UMAP projection into reference UMAP coordinates)
str(reference)
## List of 12
## $ meta_data :'data.frame': 770 obs. of 7 variables:
## ..$ cell_id : chr [1:770] "threepfresh_GATCGTACAAAGGTGC" "threepv1_CGTGATGAGTCCTC" "threepv1_AAGAATCTTGGAGG" "threepfresh_CTGAAACCACACTGCG" ...
## ..$ donor : chr [1:770] "3pv2" "3pv1" "3pv1" "3pv2" ...
## ..$ nUMI : int [1:770] 4426 4121 1535 11956 6273 1840 2702 4260 3436 3966 ...
## ..$ nGene : int [1:770] 1207 1228 560 2685 1615 584 945 1188 1360 1486 ...
## ..$ percent_mito : num [1:770] 0.0337 0.015 0.0261 0.0278 0.0289 ...
## ..$ cell_type : chr [1:770] "T_CD8" "T_CD4" "B" "DC" ...
## ..$ cell_type_broad: chr [1:770] "T" "T" "B" "DC" ...
## $ vargenes : tibble [1,600 × 3] (S3: tbl_df/tbl/data.frame)
## ..$ symbol: chr [1:1600] "LYZ" "HLA-DRA" "FTL" "CD74" ...
## ..$ mean : Named num [1:1600] 1.78 1.87 3.55 2.54 1.62 ...
## .. ..- attr(*, "names")= chr [1:1600] "LYZ" "HLA-DRA" "FTL" "CD74" ...
## ..$ stddev: Named num [1:1600] 1.91 1.69 1.22 1.53 1.83 ...
## .. ..- attr(*, "names")= chr [1:1600] "LYZ" "HLA-DRA" "FTL" "CD74" ...
## $ loadings : num [1:1600, 1:20] 0.1247 0.0811 0.1181 0.0621 0.118 ...
## $ R : num [1:100, 1:770] 7.10e-11 3.34e-02 2.47e-06 3.32e-02 5.24e-10 ...
## $ Z_orig : num [1:20, 1:770] -4.332 0.976 -3.186 0.138 -0.154 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:20] "PC_1" "PC_2" "PC_3" "PC_4" ...
## .. ..$ : chr [1:770] "4820" "10634" "8538" "4091" ...
## $ Z_corr : num [1:20, 1:770] -4.334 0.977 -3.521 -0.833 0.222 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:20] "harmony_1" "harmony_2" "harmony_3" "harmony_4" ...
## .. ..$ : chr [1:770] "4820" "10634" "8538" "4091" ...
## $ betas : num [1:3, 1:20, 1:100] -4.681788 0.011596 -0.011596 6.017041 -0.000601 ...
## $ centroids : num [1:20, 1:100] -0.21171 0.2733 0.36088 0.10265 0.00121 ...
## $ cache :List of 2
## ..$ : num [1:100, 1] 1.91 7.81 8.21 7.81 5.64 ...
## ..$ : num [1:100, 1:20] -8.9 -37.9 -34.6 -38 -21.2 ...
## $ centroids_pc :'data.frame': 100 obs. of 20 variables:
## ..$ harmony_1 : num [1:100] -4.66 -4.86 -4.21 -4.86 -3.75 ...
## ..$ harmony_2 : num [1:100] 6.019 0.792 3.614 0.793 -11.907 ...
## ..$ harmony_3 : num [1:100] 7.949 -3.73 -0.606 -3.729 4.487 ...
## ..$ harmony_4 : num [1:100] 2.261 -0.331 0.449 -0.33 -0.267 ...
## ..$ harmony_5 : num [1:100] 0.0266 0.1799 -0.3958 0.1795 1.2473 ...
## ..$ harmony_6 : num [1:100] 1.995 -0.295 0.933 -0.295 -0.557 ...
## ..$ harmony_7 : num [1:100] -0.878 -0.887 1.859 -0.889 -0.662 ...
## ..$ harmony_8 : num [1:100] -1.692 -0.858 2.703 -0.861 0.441 ...
## ..$ harmony_9 : num [1:100] -3.4883 -0.0192 1.9045 -0.0262 -1.9898 ...
## ..$ harmony_10: num [1:100] -4.18151 0.07435 -0.14573 0.07504 -0.00894 ...
## ..$ harmony_11: num [1:100] -2.1811 0.0718 -1.7295 0.0733 -0.5961 ...
## ..$ harmony_12: num [1:100] -3.3362 0.0511 0.4145 0.0513 1.5912 ...
## ..$ harmony_13: num [1:100] -2.348 -0.338 0.489 -0.337 -0.132 ...
## ..$ harmony_14: num [1:100] -6.7619 -0.0479 -0.4441 -0.0482 4.0588 ...
## ..$ harmony_15: num [1:100] 7.69 -0.241 0.328 -0.245 0.437 ...
## ..$ harmony_16: num [1:100] 0.9247 0.0268 0.7702 0.025 -1.2835 ...
## ..$ harmony_17: num [1:100] 3.9477 0.0965 0.6243 0.097 -1.0399 ...
## ..$ harmony_18: num [1:100] 9.9037 -0.0592 0.8881 -0.0584 -3.5468 ...
## ..$ harmony_19: num [1:100] 6.7289 -0.0746 -0.6049 -0.0739 2.7701 ...
## ..$ harmony_20: num [1:100] -5.959 -0.319 0.62 -0.323 -4.712 ...
## $ umap :List of 1
## ..$ embedding: num [1:770, 1:2] 4.05 3.98 -11.82 -4.81 -5.02 ...
## .. ..- attr(*, "scaled:center")= num [1:2] -0.486 -7.26
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "UMAP1" "UMAP2"
## $ save_uwot_path: chr "./testing_uwot_model_1"
The harmonized embedding is located in the Z_corr
slot of the reference object.
dim(reference$Z_corr)
## [1] 20 770
$Z_corr[1:5, 1:5] reference
## 4820 10634 8538 4091 11682
## harmony_1 -4.3340226 -4.2434240 -4.4843546 8.9367094 12.47204697
## harmony_2 0.9768383 1.2437083 -13.2156410 -2.4408398 -0.09292606
## harmony_3 -3.5205608 -3.2172170 4.6302237 0.2076788 0.29597939
## harmony_4 -0.8331292 -0.5120756 0.2940588 -0.4150630 -4.93994465
## harmony_5 0.2217745 -0.8703953 1.4084306 -3.8374885 -0.03576764
Visualize reference UMAP
= readRDS('./testing_reference1.rds')
reference = cbind(ref_metadata, reference$umap$embedding)
umap_labels plotBasic(umap_labels, title = 'Reference', color.mapping = pbmc_colors)
This option computes a reference object starting from expression in a unified pipeline, automating the preprocessing steps.
# Build reference
set.seed(0)
= symphony::buildReference(
reference
ref_exp_full,
ref_metadata,vars = c('donor'), # variables to integrate over
K = 100, # number of Harmony clusters
verbose = TRUE, # verbose output
do_umap = TRUE, # can set to FALSE if want to run umap separately later
do_normalize = FALSE, # set to TRUE if input counts are not normalized yet
vargenes_method = 'vst', # method for variable gene selection ('vst' or 'mvp')
vargenes_groups = 'donor', # metadata column specifying groups for variable gene selection
topn = 1000, # number of variable genes to choose per group
d = 20, # number of PCs
save_uwot_path = './testing_uwot_model_2' # file path to save uwot model
)
## Finding variable genes using vst method
## Total 1600 genes for downstream steps
## Scaling and PCA
## Running Harmony integration
## Harmony 1/20
## Harmony 2/20
## Harmony 3/20
## Harmony 4/20
## Harmony 5/20
## Harmony 6/20
## Harmony 7/20
## Harmony 8/20
## Harmony 9/20
## Harmony 10/20
## Harmony converged after 10 iterations
## Computing reference compression terms
## Running UMAP
## File already exists at that path... overwriting...
## Warning: invalid uid value replaced by that for user 'nobody'
## Saved uwot model
# Save reference (modify with your desired output path)
saveRDS(reference, './testing_reference2.rds')
Visualize reference UMAP
= readRDS('./testing_reference2.rds')
reference = cbind(ref_metadata, reference$umap$embedding)
umap_labels plotBasic(umap_labels, title = 'Reference', color.mapping = pbmc_colors)
In order to map a new query dataset onto the reference, you will need a reference object saved from the steps above, as well as query cell expression and metadata.
The query dataset is assumed to have been normalized in the same manner as the reference cells (here, default is log(CP10k+1) normalization).
# Read in Symphony reference to map to
= readRDS('./testing_reference1.rds')
reference # Map query
= mapQuery(query_exp, # query gene expression (genes x cells)
query # query metadata (cells x attributes)
query_metadata, # Symphony reference object
reference, do_normalize = FALSE, # perform log(CP10k) normalization on query
do_umap = TRUE) # project query cells into reference UMAP
## Scaling and synchronizing query gene expression
## Found 1600 out of 1600 reference variable genes in query dataset
## Project query cells using reference gene loadings
## Clustering query cells to reference centroids
## Correcting query batch effects
## UMAP
## All done!
Note: Symphony assumes that the query is normalized in the same manner as the reference. Our default implementation currently uses log(CP10k+1) normalization.
Let’s take a look at what the query object contains: * Z: query cells in reference Harmonized embedding * Zq_pca: query cells in pre-Harmony reference PC embedding (prior to correction) * R: query cell soft cluster assignments * Xq: query cell design matrix for correction step * umap: query cells projected into reference UMAP coordinates (using uwot) * meta_data: metadata
str(query)
## List of 7
## $ exp :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. ..@ i : int [1:99120] 1 2 3 5 10 13 14 15 26 27 ...
## .. ..@ p : int [1:431] 0 267 497 727 942 1197 1384 1608 1877 2047 ...
## .. ..@ Dim : int [1:2] 1767 430
## .. ..@ Dimnames:List of 2
## .. .. ..$ : chr [1:1767] "LYZ" "FTL" "HLA-DRA" "CD74" ...
## .. .. ..$ : chr [1:430] "fivePrime_GATCGCGAGACCGGAT" "fivePrime_AACTCAGGTTACGCGC" "fivePrime_ACATCAGGTCTTGTCC" "fivePrime_GCCTCTATCATGTCCC" ...
## .. ..@ x : num [1:99120] 2.15 3.9 5.41 3.4 3.17 ...
## .. ..@ factors : list()
## $ meta_data:'data.frame': 430 obs. of 7 variables:
## ..$ cell_id : chr [1:430] "fivePrime_GATCGCGAGACCGGAT" "fivePrime_AACTCAGGTTACGCGC" "fivePrime_ACATCAGGTCTTGTCC" "fivePrime_GCCTCTATCATGTCCC" ...
## ..$ donor : chr [1:430] "5p" "5p" "5p" "5p" ...
## ..$ nUMI : int [1:430] 6580 4750 4432 5831 5558 5360 4749 5152 4929 4194 ...
## ..$ nGene : int [1:430] 2020 1670 1668 2027 1890 1388 1661 1805 1417 1575 ...
## ..$ percent_mito : num [1:430] 0.0571 0.0564 0.0399 0.0326 0.0587 ...
## ..$ cell_type : chr [1:430] "B" "B" "B" "T_CD4" ...
## ..$ cell_type_broad: chr [1:430] "B" "B" "B" "T" ...
## $ Z : num [1:20, 1:430] -3.375 -13.148 5.845 -0.137 0.288 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:20] "harmony_1" "harmony_2" "harmony_3" "harmony_4" ...
## .. ..$ : chr [1:430] "17642" "13347" "13638" "17934" ...
## $ Zq_pca : num [1:20, 1:430] -3.089 -14.01 6.805 1.207 -0.721 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:20] "PC_1" "PC_2" "PC_3" "PC_4" ...
## .. ..$ : chr [1:430] "17642" "13347" "13638" "17934" ...
## $ R : num [1:100, 1:430] 2.90e-09 2.79e-11 3.66e-11 2.79e-11 1.67e-03 ...
## $ Xq :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## .. ..@ i : int [1:860] 0 1 0 1 0 1 0 1 0 1 ...
## .. ..@ p : int [1:431] 0 2 4 6 8 10 12 14 16 18 ...
## .. ..@ Dim : int [1:2] 2 430
## .. ..@ Dimnames:List of 2
## .. .. ..$ : NULL
## .. .. ..$ : NULL
## .. ..@ x : num [1:860] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..@ factors : list()
## $ umap : num [1:430, 1:2] -11.9 -12.26 -12.34 3.64 -11.58 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "UMAP1" "UMAP2"
Predict query cell types using k-NN. Setting confidence = TRUE also returns the prediction confidence scores (proportion of neighbors with winning vote).
= knnPredict(query, # query object
query # reference object
reference, $meta_data$cell_type, # reference cell labels for training
referencek = 5, # number of reference neighbors to use for prediction
confidence = TRUE)
Query cell type predictions are now in the cell_type_pred_knn column.
head(query$meta_data)
## cell_id donor nUMI nGene percent_mito cell_type
## 17642 fivePrime_GATCGCGAGACCGGAT 5p 6580 2020 0.05714286 B
## 13347 fivePrime_AACTCAGGTTACGCGC 5p 4750 1670 0.05642105 B
## 13638 fivePrime_ACATCAGGTCTTGTCC 5p 4432 1668 0.03993682 B
## 17934 fivePrime_GCCTCTATCATGTCCC 5p 5831 2027 0.03258446 T_CD4
## 14432 fivePrime_AGGGATGGTGGGTATG 5p 5558 1890 0.05865419 B
## 18442 fivePrime_GGGACCTCATGCCACG 5p 5360 1388 0.03750000 T_CD4
## cell_type_broad cell_type_pred_knn cell_type_pred_knn_prob
## 17642 B B 1.0
## 13347 B B 1.0
## 13638 B B 1.0
## 17934 T T_CD4 1.0
## 14432 B B 1.0
## 18442 T T_CD4 0.8
# Sync the column names for both data frames
$meta_data$cell_type_pred_knn = NA
reference$meta_data$cell_type_pred_knn_prob = NA
reference$meta_data$ref_query = 'reference'
reference$meta_data$ref_query = 'query'
query
# Add the UMAP coordinates to the metadata
= rbind(query$meta_data, reference$meta_data)
meta_data_combined = rbind(query$umap, reference$umap$embedding)
umap_combined = cbind(meta_data_combined, umap_combined)
umap_combined_labels
# Plot UMAP visualization of all cells
plotBasic(umap_combined_labels, title = 'Reference and query cells', color.by = 'ref_query')
Plot the reference and query side by side.
plotBasic(umap_combined_labels, title = 'Reference and query cells',
color.mapping = pbmc_colors, facet.by = 'ref_query')
And that’s a wrap! If you run into issues or have questions about Symphony or this tutorial, please open an issue on GitHub.
sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.5 (Santiago)
##
## Matrix products: default
## BLAS/LAPACK: /PHShome/jbk37/anaconda3/envs/seurat4/lib/libopenblasp-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 ggrastr_0.2.3 ggthemes_4.2.4 ggplot2_3.3.5
## [5] dplyr_1.0.7 data.table_1.14.0 irlba_2.3.3 Matrix_1.3-3
## [9] harmony_0.1.0 Rcpp_1.0.7 symphony_0.1.0
##
## loaded via a namespace (and not attached):
## [1] ggrepel_0.9.1 lattice_0.20-44 tidyr_1.1.3 class_7.3-19
## [5] assertthat_0.2.1 digest_0.6.27 utf8_1.2.1 RSpectra_0.16-0
## [9] R6_2.5.0 evaluate_0.14 highr_0.9 pillar_1.6.1
## [13] rlang_0.4.11 jquerylib_0.1.4 rmarkdown_2.8 labeling_0.4.2
## [17] stringr_1.4.0 munsell_0.5.0 uwot_0.1.10 compiler_4.0.5
## [21] vipor_0.4.5 xfun_0.23 pkgconfig_2.0.3 ggbeeswarm_0.6.0
## [25] htmltools_0.5.1.1 tidyselect_1.1.1 tibble_3.1.2 RANN_2.6.1
## [29] codetools_0.2-18 fansi_0.5.0 crayon_1.4.1 withr_2.4.2
## [33] grid_4.0.5 jsonlite_1.7.2 gtable_0.3.0 lifecycle_1.0.0
## [37] DBI_1.1.1 magrittr_2.0.1 scales_1.1.1 cli_3.0.1
## [41] stringi_1.6.2 farver_2.1.0 bslib_0.2.5.1 ellipsis_0.3.2
## [45] generics_0.1.0 vctrs_0.3.8 cowplot_1.1.1 RcppAnnoy_0.0.18
## [49] tools_4.0.5 Cairo_1.5-12.2 glue_1.4.2 beeswarm_0.4.0
## [53] purrr_0.3.4 yaml_2.2.1 colorspace_2.0-2 knitr_1.33
## [57] sass_0.4.0