Changes regarding python dependencies

Previous versions of APL used pytorch SVD to speed up the computation of the full SVD. This has been deprecated in favor of fast truncated SVD implementations starting with Version 1.10.1. Calling runAPL or cacomp with python = TRUE will not lead to an error, but only issue a warning. If you still want to perform a full SVD, set the dimensions to rank of the matrix. Until a faster replacement is identified, this computation will be performed by the rather slow base R svd and should therefore not be done on very large matrices. The default number of dimensions now defaults to half of the rank of the matrix.

Setup

In this vignette we will use a small data set published by Darmanis et al. (2015) consisting of 466 human adult cortical single cells sequenced on the Fluidigm platform as an example. To obtain the data necessary to follow the vignette we use the Bioconductor package scRNAseq.

Besides the package APL we will use Bioconductor packages to preprocess the data. Namely we will use SingleCellExperiment, scater and scran. However, the preprocessing could equally be performed with the single-cell RNA-seq analysis suite Seurat.

The preprocessing steps are performed according to the recommendations published in Orchestrating Single-Cell Analysis with Bioconductor by Amezquita et al. (2022). For more information about the rational behind them please refer to the book.

library(APL)
library(scRNAseq)
library(SingleCellExperiment)
library(scran)
library(scater)
set.seed(1234)

Loading the data

We start with the loading and preprocessing of the Darmanis data.

darmanis <- DarmanisBrainData()
darmanis
#> class: SingleCellExperiment 
#> dim: 22085 466 
#> metadata(0):
#> assays(1): counts
#> rownames(22085): 1/2-SBSRNA4 A1BG ... ZZZ3 tAKR
#> rowData names(0):
#> colnames(466): GSM1657871 GSM1657872 ... GSM1658365 GSM1658366
#> colData names(6): metrics age ... experiment_sample_name tissue
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):

Normalization, PCA & Clustering

Association Plots from APL should be computed based on the normalized expression data. Therefore, we first normalize the counts from the Darmanis data and calculate both PCA and UMAP for visualizations later.

For now, APL requires the data to be clustered beforehand. The darmanis data comes already annotated, so we will use the cell types stored in the cell.type metadata column instead of performing a clustering.

set.seed(100)
clust <- quickCluster(darmanis)
darmanis <- computeSumFactors(darmanis, cluster = clust, min.mean = 0.1)
darmanis <- logNormCounts(darmanis)

dec <- modelGeneVar(darmanis)
top_darmanis <- getTopHVGs(dec, n = 5000)
darmanis <- fixedPCA(darmanis, subset.row = top_darmanis)
darmanis <- runUMAP(darmanis, dimred = "PCA")

plotReducedDim(darmanis, dimred = "UMAP", colour_by = "cell.type")

Correspondence Analysis

The first step of Association Plot computations is correspondence analysis (CA). CA is a data dimensionality reduction method similar to PCA, however it allows for a simultaneous embedding of both cells and genes from the input data in the same space. In this example we perform CA on the log-normalized count matrix of the darmanis brain data.

# Computing CA on logcounts
logcounts <- logcounts(darmanis)
ca <- cacomp(
  obj = logcounts,
  top = 5000
)
#> Warning in rm_zeros(obj): Matrix contains rows with only 0s. These rows were
#> removed. If undesired set rm_zeros = FALSE.
#> No dimensions specified. Setting dimensions to: 93

# The above is equivalent to:
# ca <- cacomp(obj = darmanis,
#              assay = "logcounts",
#              top = 5000)

The function cacomp accepts as an input any matrix with non-negative entries, be it a single-cell RNA-seq, bulk RNA-seq or other data. For ease of use, cacomp accepts also SingleCellExperiment and Seurat objects, however for these we additionally have to specify via the assay and/or slot (for Seurat) parameter from where to extract the data. Importantly, in order to ensure the interpretability of the results cacomp (and related functions such as runAPL) requires that the input matrix contains both row and column names.

When performing a feature selection before CA, we can set the argument top to the desired number of genes with the highest variance across cells from the input data to retain for further analysis. By default, only the top 5,000 most variable genes are kept as a good compromise between computational time and keeping the most relevant genes. If we want to ensure however that even marker genes of smaller clusters are kept, we can increase the number of genes.

The output of cacomp is an object of class cacomp:

ca
#> cacomp object with 466 columns, 5000 rows and 93 dimensions.
#> Calc. standard coord.:  std_coords_rows, std_coords_cols
#> Calc. principal coord.: prin_coords_rows, prin_coords_cols
#> Calc. APL coord.:       
#> Explained inertia:      3.1% Dim1, 2.2% Dim2

As can be seen in the summarized output, by default both types of coordinates in the CA space (principal and standardized) are calculated. Once the coordinates for the Association Plot are calculated, they will also be shown in the output of cacomp. Slots are accessed through an accessor function:

cacomp_slot(ca, "std_coords_cols")[1:5, 1:5]
#>                  Dim1      Dim2       Dim3        Dim4       Dim5
#> GSM1657871  0.1067086 0.1253297 -0.5046485 -0.32461780 -2.4262520
#> GSM1657872 -1.5101357 0.4299418  0.1219273  0.07942823 -0.2362747
#> GSM1657873  0.2237619 0.2610148 -0.1955599  0.16963578 -2.1477040
#> GSM1657874  0.5680872 0.1138251 -0.4725071  0.04409559  0.6708562
#> GSM1657875  0.4739344 0.2505648 -0.4384626  0.10316896 -2.8910899

In the case of SingleCellExperiment and Seurat objects, we can alternatively set return_input = TRUE to get the input object back, with the CA results computed by “APL” and stored in the appropriate slot for dimension reduction. This also allows for using the plotting functions that come with these packages:

darmanis <- cacomp(
  obj = darmanis,
  assay = "logcounts",
  top = 5000,
  return_input = TRUE
)
#> Warning in rm_zeros(obj): Matrix contains rows with only 0s. These rows were
#> removed. If undesired set rm_zeros = FALSE.
#> No dimensions specified. Setting dimensions to: 93

plotReducedDim(darmanis,
  dimred = "CA",
  ncomponents = c(1, 2),
  colour_by = "cell.type"
)

plotReducedDim(darmanis,
  dimred = "CA",
  ncomponents = c(3, 4),
  colour_by = "cell.type"
)

However, some functions such as apl_coords() require information that cannot be stored in the single-cell container objects. It is therefore often easier to work with a cacomp object instead. We can convert Seurat or SingleCellExperiment objects which have CA results stored to a cacomp object using the function as.cacomp():

# Converting the object darmanis to cacomp
ca <- as.cacomp(darmanis)

Reducing the number of CA dimensions

When working with high-dimensional data, after singular value decomposition there will often be many dimensions which are representing the noise in the data. In order to minimize the noise, it is generally recommended to reduce the dimensionality of the data before generating Association Plots.

The number of dimensions to retain can be computed using the function pick_dims. This function offers three standard methods which we implemented:

• elbow rule (method = "elbow_rule") - the number of dimensions to retain is calculated based on scree plots generated for randomized data, and corresponds to a point in the plot where the band of randomized singular values enters the band of the original singular values,

• 80% rule (method = "maj_inertia") - only those first dimensions are retained which in total account for >= 80% of total inertia,

• average rule (method = "avg_inertia") - only those dimensions are retained which account for more inertia than a single dimension on average.

Additionally, the user can compute a scree plot to choose the number of dimensions by themselves:

pick_dims(ca, method = "scree_plot") +
  xlim(c(0, 20))
#> Warning: Removed 74 rows containing missing values or values outside the scale range
#> (`geom_col()`).
#> Warning: Removed 73 rows containing missing values or values outside the scale range
#> (`geom_line()`).

In the scree plot above we can see that the first dimension explains only ~1% of the total inertia and we observe the “jump” in the scree plot at roughly 5 dimensions. The first dimensions however explain only a small amount of the total inertia.

Here we compute the number of dimensions using the elbow rule. For demonstration, only three data permutations are computed:

pd <- pick_dims(
  ca,
  mat = logcounts(darmanis),
  method = "elbow_rule",
  reps = 3
)
#> Warning in rm_zeros(mat): Matrix contains rows with only 0s. These rows were
#> removed. If undesired set rm_zeros = FALSE.

pd
#> [1] 13

In this case the elbow rule leads to a higher number of dimensions.

# Compute the amount of inertia explained by each of the dimensions
D <- cacomp_slot(ca, "D")
expl_inertia <- (D^2 / sum(D^2)) * 100

# Compute the amount of intertia explained
# by the number of dimensions defined by elbow rule
sum(expl_inertia[seq_len(pd)])
#> [1] 21.14017

In this example the elbow rule suggests to keep 13 dimensions that explain 21.14% of the total inertia from the data.

Finally, we can reduce the dimensionality of the data to the desired number of dimensions:

ca <- subset_dims(ca, dims = pd)

Association Plots

When working with single-cell transcriptomics data we are often interested in which genes are associated to a cluster of cells. To reveal such genes we can compute an Association Plot for a selected cluster of cells. In the following example we want to generate an Association Plot for the cluster of endothelial cells:

# Specifying a cell cluster of interest
endo <- which(darmanis$cell.type == "endothelial")

# Calculate Association Plot coordinates for endothelial cells
ca <- apl_coords(ca, group = endo)

After computing the coordinates of genes and cells in the Association Plot we are able to plot the results using the apl function.

# endothelial marker genes
marker_genes <- c("APOLD1", "TM4SF1", "SULT1B1", "ESM1", "SELE")

# Plot APL
apl(ca,
  row_labs = TRUE,
  rows_idx = marker_genes,
  type = "ggplot"
) # type = "plotly" for an interactive plot

In the Association Plot all genes are represented by blue circles. The further to the right a gene is located the more associated it is with the chosen cluster of cells and the lower the y-axis value, the more specific it is for the selected cluster. Additionally, it is possible to highlight in the Association Plot any set of genes. In the example above we highlighted five genes (APOLD1, TM4SF1, SULT1B1, ESM1, SELE) which are known to be marker genes for endothelial cells. As we can see in the plot, they are located in the right part of the plot, which confirms their specificity for endothelial cells.

By default we plot only the genes in the Association Plot. To also display the cells in the Association Plot, use the argument show_cols = TRUE. This way we can identify other cells which show similar expression profiles to the cells of interest. Cells that belong to the cluster of interest will be colored in red, and all remaining cells will be colored in violet. Furthermore, an interactive plot in which you can hover over genes to see their name can be created by setting type = "plotly".

Association Plots with the S_α-scores

The S_α-score allows us to rank genes by their specificity for a selected cell cluster, and is computed for each gene from the Association Plot separately. The higher the S_α-score of a gene, the more characteristic its expression for the investigated cell cluster. The S_α-scores can be computed using the apl_score function. To display the S_α-scores in the Association Plot, we can use the argument show_score = TRUE in the apl function:

# Compute S-alpha score
# For the calculation the input matrix is also required.
ca <- apl_score(ca,
  mat = logcounts(darmanis),
  reps = 5
)

apl(ca,
  show_score = TRUE,
  type = "ggplot"
)

By default, only genes that have a S_α-score larger than 0 are colored as these tend to be genes of interest and we consider them as cluster-specific genes. This cutoff can be easily changed through the score_cutoff argument to apl().

The S_α-scores are stored in the "APL_score" slot and can be accessed as follows:

head(cacomp_slot(ca, "APL_score"))
#>             Rowname    Score Row_num Rank
#> IFITM1       IFITM1 2.750320    2081    1
#> LOC646268 LOC646268 2.743558      81    2
#> OR2A20P     OR2A20P 2.743558    4888    3
#> OR4C46       OR4C46 2.743558    4889    4
#> PRKCDBP     PRKCDBP 2.722349    2390    5
#> TWIST1       TWIST1 2.637943    1688    6

To see the expression of genes with the highest S_α-scores (or any selected genes) across all cell types from the data we can use plotting functions provided by scater:

scores <- cacomp_slot(ca, "APL_score")

plotExpression(darmanis,
  features = head(scores$Rowname, 3),
  x = "cell.type",
  colour_by = "cell.type"
)

plotReducedDim(darmanis,
  dimred = "UMAP",
  colour_by = scores$Rowname[1]
)

As expected, the 3 most highly scored genes are over-expressed in the endothelial cells. Due to the small size of the data set and number of cells in the cluster (only 20 out of 466 cells are endothelial cells) some cluster specific genes are only expressed in a few cells. Most data sets nowadays are significantly larger so this should not be a major concern and it can further be mitigated by performing a more stringent feature selection before CA.

Visualization of CA

In addition to Association Plots “APL” produces also other forms of the output. For instance, we can use “APL” to generate a two- and three-dimensional correspondence analysis projection of the data. The so-called biplot visualizes both cells and genes from the input data and can be created using the function ca_biplot. Alternatively, a three-dimensional data projection plot can be generated using the function ca_3Dplot. To generate such biplots a cacomp object is required.

# Specifying a cell cluster of interest
endo <- which(darmanis$cell.type == "endothelial")

# Creating a static plot
ca_biplot(ca, col_labels = endo, type = "ggplot")
#> Warning in ggplot2::geom_point(data = rows, ggplot2::aes_(x = as.name(rnmx), :
#> Ignoring unknown aesthetics: text
#> Warning in ggplot2::geom_point(data = cols, ggplot2::aes_(x = as.name(cnmx), :
#> Ignoring unknown aesthetics: text
#> Warning in ggplot2::geom_point(data = cols[col_labels, ], ggplot2::aes_(x =
#> as.name(cnmx), : Ignoring unknown aesthetics: text

# Creating an interactive plot
# ca_biplot(ca, type = "plotly", col_labels = platelets)

# 3D plot
# ca_3Dplot(ca, col_labels = platelets)

The above described plots give us a quick overview of the first 2 dimensions of the data (more dimensions can be plotted). As shown in the commented-out code, to interactively explore the projection of the data type = "plotly" can be set.