Functions for Evaluation of Dataset Alignment

• compareCCA(): This function enables the comparison of datasets by applying Canonical Correlation Analysis (CCA). It helps assess how well two datasets align with each other, providing insights into the relationship between different single-cell experiments or conditions.

• comparePCA(): This function allows you to compare datasets using Principal Component Analysis (PCA). It evaluates how similar or different the principal components are between two datasets, offering a way to understand the underlying structure and variance in your data.

• comparePCASubspace(): Extending the comparison to specific subspaces, this function focuses on subsets of principal components. It provides a detailed analysis of how subspace structures differ or align between datasets, which is valuable for fine-grained comparative studies.

• plotPairwiseDistancesDensity(): To visualize the distribution of distances between pairs of samples, this function generates density plots. It helps in understanding the variation and relationships between samples in high-dimensional spaces.

• plotWassersteinDistance(): This function visualizes the Wasserstein distance, a metric for comparing distributions, across datasets. It provides an intuitive view of how distributions differ between datasets, aiding in the evaluation of alignment and discrepancies.

Statistical Measures to Assess Dataset Alignment

• calculateCramerPValue(): This function computes the Cramér’s V statistic, which measures the strength of association between categorical variables. Cramér’s V is particularly useful for quantifying the association strength in contingency tables, providing insight into how strongly different cell types or conditions are related.

• calculateHotellingPValue(): This function performs Hotelling’s T-squared test, a multivariate statistical test used to assess the differences between two groups of observations. It is useful for comparing the mean vectors of multiple variables and is commonly employed in single-cell data to evaluate group differences in principal component space.

• calculateNearestNeighborProbabilities(): This function calculates the probability of nearest neighbor assignments in a high-dimensional space. It is used to evaluate the likelihood of a cell belonging to a specific cluster based on its proximity to neighboring cells. This function helps in assessing clustering quality and the robustness of cell type assignments.

• calculateAveragePairwiseCorrelation(): This function computes the average pairwise correlation between variables or features across cells. It provides an overview of the relationships between different markers or genes, helping to identify correlated features and understand their interactions within the dataset.

• regressPC(): This function performs linear regression of a covariate of interest onto one or more principal components derived from single-cell data. This analysis helps quantify the variance explained by the covariate and is useful for tasks such as assessing batch effects, clustering homogeneity, and alignment of query and reference datasets.

Marker Gene Alignment

• calculateHVGOverlap(): To assess the similarity between datasets based on highly variable genes, this function computes the overlap coefficient. It measures how well the sets of highly variable genes from different datasets correspond to each other.

• calculateVarImpOverlap(): Using Random Forest, this function identifies and compares the importance of genes for differentiating cell types between datasets. It highlights which genes are most critical in each dataset and compares their importance, providing insights into shared and unique markers.

Purpose and Applications

These functions collectively offer a comprehensive toolkit for comparing and analyzing single-cell data. Whether you are assessing alignment between datasets, visualizing distance distributions, or identifying key genes, these tools are designed to enhance your ability to derive meaningful insights from complex, high-dimensional data.

In this vignette, we will guide you through the practical use of each function, demonstrate how to interpret their outputs, and show how they can be integrated into your single-cell genomics workflow.

compareCCA()

In single-cell genomics, datasets from different sources or experimental conditions often need to be compared to understand how well they align. The compareCCA function facilitates this comparison by performing Canonical Correlation Analysis (CCA) between a query dataset and a reference dataset. This function is particularly useful when datasets have different sample sizes or distributions, and it helps in assessing the similarity between them after projecting them onto a common principal component space.

compareCCA() performs the following steps:

• Project Data: Project the query dataset onto the principal component space of the reference dataset. This ensures that both datasets are represented in the same reduced-dimensional space for a meaningful comparison.
• Perform CCA: Conduct Canonical Correlation Analysis to find the relationships between the projected datasets. CCA helps to identify how the canonical variables derived from the query dataset correlate with those from the reference dataset.
• Compare Alignment: Evaluate the alignment between the datasets by computing cosine similarities between the canonical correlation vectors derived from the CCA. This provides a measure of how well the two datasets correspond to each other in the canonical space.

# Perform CCA 
cca_comparison <- compareCCA(query_data = query_data_cd4, 
                             reference_data = ref_data_cd4, 
                             query_cell_type_col = "expert_annotation", 
                             ref_cell_type_col = "expert_annotation", 
                             pc_subset = 1:5)
plot(cca_comparison)

Canonical Correlation Analysis (CCA) produces several key outputs:

• Canonical Variables: For each dataset, CCA computes a set of canonical variables that maximize the correlation between the datasets. These variables are linear combinations of the original variables.
• Canonical Correlation Coefficients: These coefficients indicate the strength of the correlation between the canonical variables of the two datasets. Higher values suggest stronger correlations.
• Cosine Similarities: After obtaining the canonical variables, the function computes the cosine similarity between these variables from the query and reference datasets. This similarity measure helps assess how closely the canonical variables align between the two datasets.

In summary, the compareCCA function allows you to compare how well two datasets are aligned by projecting them into a shared PCA space, performing CCA, and then evaluating the similarity of the canonical variables. This approach is valuable for integrative analyses and understanding the relationships between different datasets in single-cell studies.

comparePCA()

The comparePCA() function compares the PCA subspaces between the query and reference datasets. It calculates the principal angles between the PCA subspaces to assess the alignment and similarity between them. This is useful for understanding how well the dimensionality reduction spaces of your datasets match.

comparePCA() operates as follows:

• PCA Computation: It computes PCA for both query and reference datasets, reducing them into lower-dimensional spaces.
• Subspace Comparison: The function calculates the principal angles between the PCA subspaces of the query and reference datasets. These angles help determine how closely the subspaces align.
• Distance Metrics: It uses distance metrics based on principal angles to quantify the similarity between the datasets.

# Perform PCA 
pca_comparison <- comparePCA(query_data = query_data_cd4, 
                             reference_data = ref_data_cd4, 
                             query_cell_type_col = "expert_annotation", 
                             ref_cell_type_col = "expert_annotation", 
                             pc_subset = 1:5, 
                             metric = "cosine")
plot(pca_comparison)

comparePCASubspace()

In single-cell RNA-seq analysis, it is essential to assess the similarity between the subspaces spanned by the top principal components (PCs) of query and reference datasets. This is particularly important when comparing the structure and variation captured by each dataset. The comparePCASubspace() function is designed to provide insights into how well the subspaces align by computing the cosine similarity between the loadings of the top variables for each PC. This analysis helps in determining the degree of correspondence between the datasets, which is critical for accurate cell type annotation and data integration.

comparePCASubspace() performs the following operations:

• Cosine Similarity Calculation: The function computes the cosine similarity between the loadings of the top variables for each PC in both the query and reference datasets. This similarity measures how closely the two datasets align in the space spanned by their principal components.
• Selection of Top Similarities: The function selects the top cosine similarity scores and stores the corresponding PC indices from both datasets. This step identifies the most aligned principal components between the two datasets.
• Variance Explained Calculation: It then calculates the average percentage of variance explained by the selected top PCs. This helps in understanding how much of the data’s variance is captured by these components.
• Weighted Cosine Similarity Score: Finally, the function computes a weighted cosine similarity score based on the top cosine similarities and the average percentage of variance explained. This score provides an overall measure of subspace alignment between the datasets.

# Compare PCA subspaces between query and reference data
subspace_comparison <- comparePCASubspace(
    query_data = query_data_cd4,
    reference_data = ref_data_cd4, 
    query_cell_type_col = "expert_annotation", 
    ref_cell_type_col = "expert_annotation", 
    pc_subset = 1:5
)

# View weighted cosine similarity score
subspace_comparison$weighted_cosine_similarity
#> [1] 0.2578375

# Plot output for PCA subspace comparison (if a plot method is available)
plot(subspace_comparison)

In the results:

• Cosine Similarity: The values in principal_angles_cosines indicate the degree of alignment between the principal components of the query and reference datasets. Higher values suggest stronger alignment.
• Variance Explained: The average_variance_explained vector provides the average percentage of variance captured by the selected top PCs in both datasets. This helps assess the importance of these PCs in explaining data variation.
• Weighted Cosine Similarity: The weighted_cosine_similarity score combines the cosine similarity with variance explained to give a comprehensive measure of how well the subspaces align. A higher score indicates that the datasets are well-aligned in the PCA space.

By using comparePCASubspace(), you can quantify the alignment of PCA subspaces between query and reference datasets, aiding in the evaluation of dataset integration and the reliability of cell type annotations.

plotPairwiseDistancesDensity()

# Example usage of the function
plotPairwiseDistancesDensity(query_data = query_data, 
                             reference_data = reference_data, 
                             query_cell_type_col = "expert_annotation", 
                             ref_cell_type_col = "expert_annotation", 
                             cell_type_query = "CD4", 
                             cell_type_ref = "CD4", 
                             pc_subset = 1:5,
                             distance_metric = "euclidean")

calculateWassersteinDistance()

the calculateWassersteinDistance() function uses the concept of Wasserstein distances, a measure rooted in optimal transport theory, to quantitatively compare these datasets. By projecting both datasets into a common PCA space, the function enables a meaningful comparison even when they originate from different experimental conditions.

Before diving into a code example, it’s important to understand the output of this function. The result is a list containing three key components: the null distribution of Wasserstein distances computed from the reference dataset alone, the mean Wasserstein distance between the query and reference datasets, and a vector of the unique cell types present in the reference dataset. This comparison provides insights into how similar or different the query dataset is relative to the reference dataset.

What the function reveals is particularly valuable: if the mean Wasserstein distance between the query and reference is significantly different from the null distribution, it indicates a notable variation between the datasets, reflecting differences in cell-type compositions, technical artifacts, or biological questions of interest.

# Generate the Wasserstein distance density plot
wasserstein_data <- calculateWassersteinDistance(
    query_data = query_data_cd4,
    reference_data = ref_data_cd4, 
    query_cell_type_col = "expert_annotation", 
    ref_cell_type_col = "expert_annotation", 
    pc_subset = 1:10,
)
plot(wasserstein_data)

calculateCramerPValue()

The calculateCramerPValue function is designed to perform the Cramer test for comparing multivariate empirical cumulative distribution functions (ECDFs) between two samples in single-cell data. This test is particularly useful for assessing whether the distributions of principal components (PCs) differ significantly between the reference and query datasets for specific cell types.

To use this function, you first need to provide two key inputs: reference_data and query_data, both of which should be SingleCellExperiment objects containing numeric expression matrices. If query_data is not supplied, the function will only use the reference_data. You should also specify the column names for cell type annotations in both datasets via ref_cell_type_col and query_cell_type_col. If cell_types is not provided, the function will automatically include all unique cell types found in the datasets. The pc_subset parameter allows you to define which principal components to include in the analysis, with the default being the first five PCs.

The function performs the following steps: it first projects the data into PCA space, subsets the data according to the specified cell types and principal components, and then applies the Cramer test to compare the ECDFs between the reference and query datasets. The result is a named vector of p-values from the Cramer test for each cell type, which indicates whether there is a significant difference in the distributions of PCs between the two datasets.

Here’s an example of how to use the calculateCramerPValue() function:

# Perform Cramer test for the specified cell types and principal components
cramer_test <- calculateCramerPValue(
    reference_data = reference_data,
    query_data = query_data,
    ref_cell_type_col = "expert_annotation", 
    query_cell_type_col = "SingleR_annotation",
    cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"),
    pc_subset = 1:5
)

# Display the Cramer test p-values
print(cramer_test)
#>          CD4          CD8 B_and_plasma      Myeloid 
#>    0.0000000    0.0000000    0.1378621    0.4575425

In this example, the function compares the distributions of the first five principal components between the reference and query datasets for cell types such as “CD4”, “CD8”, “B_and_plasma”, and “Myeloid”. The output is a vector of p-values indicating whether the differences in PC distributions are statistically significant.

calculateHotellingPValue()

The calculateHotellingPValue() function is designed to compute Hotelling’s T-squared test statistic and corresponding p-values for comparing multivariate means between reference and query datasets in the context of single-cell RNA-seq data. This statistical test is particularly useful for assessing whether the mean vectors of principal components (PCs) differ significantly between the two datasets, which can be indicative of differences in the cell type distributions.

To use this function, you need to provide two SingleCellExperiment objects: query_data and reference_data, each containing numeric expression matrices. You also need to specify the column names for cell type annotations in both datasets with query_cell_type_col and ref_cell_type_col. The cell_types parameter allows you to choose which cell types to include in the analysis, and if not specified, the function will automatically include all cell types present in the datasets. The pc_subset parameter determines which principal components to consider, with the default being the first five PCs. Additionally, n_permutation specifies the number of permutations for calculating p-values, with a default of 500.

The function works by first projecting the data into PCA space and then performing Hotelling’s T-squared test for each specified cell type to compare the means between the reference and query datasets. It uses permutation testing to determine the p-values, indicating whether the observed differences in means are statistically significant. The result is a named numeric vector of p-values for each cell type.

Here is an example of how to use the calculateHotellingPValue() function:

# Calculate Hotelling's T-squared test p-values
p_values <- calculateHotellingPValue(
    query_data = query_data, 
    reference_data = reference_data, 
    query_cell_type_col = "SingleR_annotation", 
    ref_cell_type_col = "expert_annotation",
    pc_subset = 1:10
)

# Display the p-values rounded to five decimal places
print(round(p_values, 5))
#>          CD4          CD8 B_and_plasma      Myeloid 
#>         0.00         0.00         0.09         0.00

In this example, the function calculates p-values for Hotelling’s T-squared test, focusing on the first ten principal components to assess whether the multivariate means differ significantly between the reference and query datasets for each cell type. The p-values indicate the likelihood of observing the differences by chance and help identify significant disparities in cell type distributions.

calculateAveragePairwiseCorrelation()

The calculateAveragePairwiseCorrelation() function is designed to compute the average pairwise correlations between specified cell types in single-cell gene expression data. This function operates on SingleCellExperiment objects and is ideal for single-cell analysis workflows. It calculates pairwise correlations between query and reference cells using a specified correlation method, and then averages these correlations for each cell type pair. This helps in assessing the similarity between cells in the reference and query datasets and provides insights into the reliability of cell type annotations.

To use the calculateAveragePairwiseCorrelation() function, you need to supply it with two SingleCellExperiment objects: one for the query cells and one for the reference cells. The function also requires column names specifying cell type annotations in both datasets, and optionally a vector of cell types to include in the analysis. Additionally, you can specify a subset of principal components to use, or use the raw data directly if pc_subset is set to NULL. The correlation method can be either “spearman” or “pearson”.

Here’s an example of how to use this function:

# Compute pairwise correlations between specified cell types
cor_matrix_avg <- calculateAveragePairwiseCorrelation(
  query_data = query_data, 
  reference_data = reference_data, 
  query_cell_type_col = "SingleR_annotation", 
  ref_cell_type_col = "expert_annotation", 
  cell_types = c("CD4", "CD8", "B_and_plasma"), 
  pc_subset = 1:5,
  correlation_method = "spearman"
)

# Visualize the average pairwise correlation matrix
plot(cor_matrix_avg)

In this example, calculateAveragePairwiseCorrelation() computes the average pairwise correlations for the cell types “CD4”, “CD8”, and “B_and_plasma”. This function is particularly useful for understanding the relationships between different cell types in your single-cell datasets and evaluating how well cell types in the query data correspond to those in the reference data.

regressPC()

The regressPC() function performs linear regression of a covariate of interest onto one or more principal components using data from a SingleCellExperiment object. This method helps quantify the variance explained by a covariate, which can be useful in applications such as quantifying batch effects, assessing clustering homogeneity, and evaluating alignment between query and reference datasets in cell type annotation settings.

The function calculates the R-squared value from the linear regression of the covariate onto each principal component. The variance contribution of the covariate effect per principal component is computed as the product of the variance explained by the principal component and the R-squared value. The total variance explained by the covariate is obtained by summing these contributions across all principal components.

Here is how you can use the regressPC() function:

# Perform regression analysis using only reference data
regress_res <- regressPC(
    reference_data = reference_data,
    ref_cell_type_col = "expert_annotation",
    cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"),
    pc_subset = 1:15
)

# Plot results showing R-squared values
plot(regress_res, plot_type = "r_squared")

# Plot results showing p-values
plot(regress_res, plot_type = "p-value")

In this example, regressPC() is used to perform regression analysis on principal components 1 to 15 from the reference dataset. The results are then visualized using plots showing the R-squared values and p-values.

If you also have a query dataset and want to compare it with the reference dataset, you can include it in the analysis:

# Perform regression analysis using both reference and query data
regress_res <- regressPC(
    reference_data = reference_data,
    query_data = query_data,
    ref_cell_type_col = "expert_annotation",
    query_cell_type_col = "SingleR_annotation",
    cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"),
    pc_subset = 1:15
)

# Plot results showing R-squared values
plot(regress_res, plot_type = "r_squared")

# Plot results showing p-values
plot(regress_res, plot_type = "p-value")

In this case, the function projects the query data onto the principal components of the reference data and performs regression analysis. The results are visualized in the same way, providing insights into how well the principal components from the query data align with those from the reference data.

This function is useful for understanding the impact of different covariates on principal components and for comparing how cell types in different datasets are represented in the principal component space.

calculateHVGOverlap()

The calculateHVGOverlap() function computes the overlap coefficient between two sets of highly variable genes (HVGs) from a reference dataset and a query dataset. The overlap coefficient is a measure of similarity between the two sets, reflecting how much the HVGs in one dataset overlap with those in the other.

# Load library to get top HVGs
library(scran)

# Select the top 500 highly variable genes from both datasets
ref_var_genes <- getTopHVGs(reference_data, n = 500)
query_var_genes <- getTopHVGs(query_data, n = 500)

# Calculate the overlap coefficient between the reference and query HVGs
overlap_coefficient <- calculateHVGOverlap(reference_genes = ref_var_genes, 
                                           query_genes = query_var_genes)

# Display the overlap coefficient
overlap_coefficient
#> [1] 0.93

calculateVarImpOverlap()

# RF function to compare which genes are best at differentiating cell types
rf_output <- calculateVarImpOverlap(reference_data = reference_data, 
                                    query_data = query_data, 
                                    query_cell_type_col = "expert_annotation", 
                                    ref_cell_type_col = "expert_annotation", 
                                    n_tree = 500,
                                    n_top = 50)

# Comparison table
rf_output$var_imp_comparison
#>              CD4-CD8     CD4-B_and_plasma          CD4-Myeloid 
#>                 0.84                 0.84                 0.88 
#>     CD8-B_and_plasma          CD8-Myeloid B_and_plasma-Myeloid 
#>                 0.84                 0.80                 0.76