• About
• Part 1 — Data Loading & Metadata
• 1. Loading the Data
• 2. Create the sample sheet
• 3. Add meta-data
  • 3.1. How many cells are in each sample?
  • 3.2. How many genes are in the expression matrices?
  • 3.3. What information is now part of the meta-data of the objects?
• Part 2 — Preprocessing, Batch Correction & Clustering
• 4. Preprocessing
  • 4.1 Preprocessing
    • 4.1.1. Correct order of preprocessing steps:
  • 4.1.2 Steps Before and After Merging (Task 4.2):
  • 4.1.3. Performing preprocessing on the data
    • 4.1.3.1. Filtering
    • Mitochondrial Content
    • 4.1.4. Choosing Cut-off parameters
    • 4.1.5. Apply Filtering
  • 4.1.6. Doublet removal with DoubletFinder
    • 4.1.6.1. Explanation of why we are performing doublet removal
    • 4.1.6.2. Normalization
    • 4.1.6.3. Feature Selection
    • 4.1.6.4. Scale
    • 4.1.6.5. Run PCA
    • 4.1.6.6. Elbow Point
    • 4.1.6.7. Choosing the Number of PCs
  • 4.2 Merging and Batch Correction
    • 4.2.1 1. Merge All Four Samples Without Batch Correction
    • 4.2.2 Merge All Four Samples Using Seurat’s Data Integration Method (Batch Correction)
    • 4.2.4. Batch effects in single-cell data can arise from various sources, including:
  • 4.1.8. Re-Feature Selection
    • 4.1.8.1. Purpose
    • 4.1.8.2. How Are Features Selected?
• 5. Dimensionality Reduction
  • 5.1. Dimensionality Reduction
    • 5.1.1. Scaling the Data
    • 5.1.2 Running PCA
    • 5.1.3 Determining the Number of Principal Components
    • 5.1.3.1 Choosing elbow point (number of dimensions)
  • 5.1.4. Why Use PCA with UMAP
  • 5.2. Clustering
    • 5.2.2. Run UMAP
    • 5.2.1. Visualization of Clusters
    • Was batch correction necessary:
    • 4.2.5. Batch correction is crucial
  • 5.2.2. Keeping Results for Subsequent Tasks
• Part 3 — Cell-Type Annotation, Differential Expression & Pathways
• 6. Cell Type Annotation
  • 6.1. Automatic Cell Type Annotation with SingleR
  • 6.2. Manual Annotation
    • 6.2.1. Perform Differential Expression Analysis for Cell-Type Annotation
    • 6.2.2. Use Markers from Table 2 to Identify Cell Types and Assign Names to Clusters
  • Strategy for Manual cell annotation
    • Calculate Module Scores Using AddModuleScore
    • Assign Cell Types Based on Highest Module Score
    • 3. Visualize the UMAP with Assigned Cell Types
    • 6.2.3 Automatic vs Manual Annotation Comparison
  • 6.3. Cell-Type Proportions (Bonus)
    • 6.3.1. Compute Cell-Type Proportions for Each Sample
    • 6.3.2. Plot Cell-Type Proportions
    • 6.3.3. Explain Sample Variability in Cell-Type Proportions
• 7. Differential Analysis
  • 7.1. Differential Expression Analysis on cell-types
    • 7.1.1. Subset the Data for Relevant Cell Types
    • 7.1.3. Perform Differential Expression Analysis
    • 7.1.3.4 Prepare Data for Volcano Plots
    • 7.1.3.5. Create Volcano Plots
    • 7.1.1. Memory Formation on Cells
  • 7.2. Plot Differentially Expressed Genes
• 8. Pathway Analysis
  • 8.1. Differential Expression Analysis on groups
    • 8.1.1. Differential Expression Analysis for BMMC vs CD34 (All Cell Types)
    • 8.1.2. Extract Top 5 DEGs
    • 8.1.3. Differential Expression Analysis for Monocytes in BMMC vs CD34
  • 8.2. Pathway analysis on groups
  • 8.3. Biological interpretation
• Part 4 — Trajectory & Cell-Cell Communication
• 9. Trajectory Analysis
  • 9.1. Select Subset
    • 9.1.1. Selecting a Group of Cells for Trajectory Analysis
    • 9.1.2. Performing Trajectory Analysis with Monocle 3
    • 9.1.2.2. Preprocess
  • 9.1.3. UMAP, Clustering and Visualisation
    • 9.1.5. Plot Trajectory Colored by Cluster
  • 9.2. Select Root Nodes Manually
  • 9.2.1 Selecting Root Nodes and Performing Trajectory Analysis
    • 9.2.1.1. Selecting Root Nodes
    • Observation
    • Cluster 2
    • Observation:
  • 9.3 Select Root Nodes Automatically
    • Observation
    • Tracing the Biological Path
    • Manual vs Automatic root node selection
• 10. Cell-Cell Communication
  • 10.0. Subsetting Common Cell Types
  • 10.1. Create CellChat objects for each group using the common cell types
  • 10.3. Set the CellChat database, preprocesssing with visualization.
  • Common pathway
  • Interaction and Interaction Strength
• 11. Summary
• Session Info

About

Author: Md Mobashir Rahman

This notebook analyzes four human bone marrow scRNA-seq samples (Granja et al. 2019) end-to-end: QC and filtering, doublet removal, normalization, batch correction, dimensionality reduction and clustering, automatic and manual cell-type annotation, differential expression, pathway enrichment, trajectory analysis, and cell-cell communication. It originated as a project for the Single Cell Bioinformatics course at Saarland University and has since been cleaned up for portfolio presentation. See the repository README for a narrative summary and key figures.

Acknowledgement: Modern GPT tools (ChatGPT, Perplexity AI and GitHub Copilot) were frequently used to explore strategies, troubleshoot code, and as a replacement for search engines.

knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  fig.dpi = 130
)

# Set a global seed for reproducibility
set.seed(1234)

# Check all required libraries are installed
suppressPackageStartupMessages({
library(dplyr)
library(spatstat.core)
library(Seurat)
library(patchwork)
library(DoubletFinder)
library(SingleR)
library(enrichR)
library(CellChat)
library(SingleCellExperiment)
library(SeuratWrappers)
library(monocle3)
library(celldex)
library(knitr)
})

Part 1 — Data Loading & Metadata

1. Loading the Data

We will load the expression matrices from the dataset and construct a Seurat object. We will need to load two files: one containing data on Bone Marrow Mononuclear Cells (BMMC) and the other on CD34+ Enriched Bone Marrow Cells (CD34).

# Load the .rds files

# Bone Marrow Mononuclear Cells (BMMC)
bmmc_d1t1 <- readRDS("../data/GSM4138872_scRNA_BMMC_D1T1.rds")
bmmc_d1t2 <- readRDS("../data/GSM4138873_scRNA_BMMC_D1T2.rds")

# CD34+ Enriched Bone Marrow Cells (CD34)
cd34_d2t1 <- readRDS("../data/GSM4138874_scRNA_CD34_D2T1.rds")
cd34_d3t1 <- readRDS("../data/GSM4138875_scRNA_CD34_D3T1.rds")

2. Create the sample sheet

We will label each sample with the corresponding metadata from the table 1 provided in the assignment sheet.

# Define the sample metadata based on Table 1
metadata <- data.frame(
  Sample = c("BMMC_D1T1", "BMMC_D1T2", "CD34_D2T1", "CD34_D3T1"),
  Donor = c("D1", "D1", "D2", "D3"),
  Replicate = c("T1", "T2", "T1", "T1"),
  Sex = c("F", "F", "M", "F"),
  row.names = c("BMMC_D1T1", "BMMC_D1T2", "CD34_D2T1", "CD34_D3T1")
)

3. Add meta-data

# Create Seurat objects with metadata added
bmmc_d1t1_seurat <- CreateSeuratObject(counts = bmmc_d1t1, project = "BMMC_D1T1")
bmmc_d1t1_seurat@meta.data$Sample <- "BMMC_D1T1"
bmmc_d1t1_seurat@meta.data$Donor <- "D1"
bmmc_d1t1_seurat@meta.data$Replicate <- "T1"
bmmc_d1t1_seurat@meta.data$Sex <- "F"

bmmc_d1t2_seurat <- CreateSeuratObject(counts = bmmc_d1t2, project = "BMMC_D1T2")
bmmc_d1t2_seurat@meta.data$Sample <- "BMMC_D1T2"
bmmc_d1t2_seurat@meta.data$Donor <- "D1"
bmmc_d1t2_seurat@meta.data$Replicate <- "T2"
bmmc_d1t2_seurat@meta.data$Sex <- "F"

cd34_d2t1_seurat <- CreateSeuratObject(counts = cd34_d2t1, project = "CD34_D2T1")
cd34_d2t1_seurat@meta.data$Sample <- "CD34_D2T1"
cd34_d2t1_seurat@meta.data$Donor <- "D2"
cd34_d2t1_seurat@meta.data$Replicate <- "T1"
cd34_d2t1_seurat@meta.data$Sex <- "M"

cd34_d3t1_seurat <- CreateSeuratObject(counts = cd34_d3t1, project = "CD34_D3T1")
cd34_d3t1_seurat@meta.data$Sample <- "CD34_D3T1"
cd34_d3t1_seurat@meta.data$Donor <- "D3"
cd34_d3t1_seurat@meta.data$Replicate <- "T1"
cd34_d3t1_seurat@meta.data$Sex <- "F"

For each sample, we report the following information:

3.1. How many cells are in each sample?

3.2. How many genes are in the expression matrices?

# Function to extract and print required information 
report_sample_info <- function(seurat_object, sample_name) {
  cat(sample_name, ":\n")
  
  # 1. Report the number of cells
  num_cells <- ncol(seurat_object)
  cat("Number of cells in the sample:", num_cells, "\n")
  
  # 2. Report the number of genes
  num_genes <- nrow(seurat_object)
  cat("Number of genes in the sample:", num_genes, "\n")
  
}

# Report information for each sample
report_sample_info(bmmc_d1t1_seurat, "BMMC_D1T1")

## BMMC_D1T1 :
## Number of cells in the sample: 6270 
## Number of genes in the sample: 20287

report_sample_info(bmmc_d1t2_seurat, "BMMC_D1T2")

## BMMC_D1T2 :
## Number of cells in the sample: 6332 
## Number of genes in the sample: 20287

report_sample_info(cd34_d2t1_seurat, "CD34_D2T1")

## CD34_D2T1 :
## Number of cells in the sample: 2424 
## Number of genes in the sample: 20287

report_sample_info(cd34_d3t1_seurat, "CD34_D3T1")

## CD34_D3T1 :
## Number of cells in the sample: 5752 
## Number of genes in the sample: 20287

# Creating a named list of Seurat objects
seurat_list <- list(
  bmmc_d1t1_seurat = bmmc_d1t1_seurat,
  bmmc_d1t2_seurat = bmmc_d1t2_seurat,
  cd34_d2t1_seurat = cd34_d2t1_seurat,
  cd34_d3t1_seurat = cd34_d3t1_seurat
)

3.3. What information is now part of the meta-data of the objects?

Column name present in the metadata of all samples:

orig.ident: this often contains the sample identity if known, in our case, we know the identity of the samples.

nCount_RNA: Number of UMIs per cell.

nFeature_RNA: number of genes detected per cell.

Sample: Sample names comprising of cell type, donor ID and replicate ID.

Donor: The Donor column identifies the individual source of the sample (e.g., D1, D2, D3).

Replicate: The Replicate column (e.g., T1, T2) indicates separate samples collected from the same donor. For example: BMMC_D1T1 and BMMC_D1T2 both come from the same donor (D1), but are labeled as T1 and T2, indicating they are separate replicates.

Sex: Biological sexual identity of the origin human subject.

Part 2 — Preprocessing, Batch Correction & Clustering

4. Preprocessing

4.1 Preprocessing

4.1.1. Correct order of preprocessing steps:

1. Filtering is the initial step to remove low-quality cells and genes from each dataset. This ensures that only high-quality data is used for downstream analyses.

2. Normalization adjusts for technical variations such as differences in sequencing depth between cells, making the gene expression levels comparable across all cells.

3. Feature Selection identifies the most variable genes across the dataset, which are most informative for downstream analyses like dimensionality reduction and clustering.

4. Doublet Removal is performed after normalization and feature selection because tools like DoubletFinder require normalized data and principal component analysis (PCA) results to accurately detect doublets. (as documented in its Github docuementation).

5. Batch Correction/Merging

6. Re-feature Selection

7. Re-scaling

Then proceed to dimensionality reduction.

4.1.2 Steps Before and After Merging (Task 4.2):

After removing doublet (running doubletfinder on individual samples separatedly), we merged the individual samples.

4.1.3. Performing preprocessing on the data

4.1.3.1. Filtering

Before filtering, lets get some QC Metrics.

4.1.3.1.1. QC Metrics

# Merging the dataset to get the QC metrics.

merged_seurat <- merge(
  x = bmmc_d1t1_seurat,
  y = list(bmmc_d1t2_seurat, cd34_d2t1_seurat, cd34_d3t1_seurat),
  add.cell.ids = c("bmmc_d1t1", "bmmc_d1t2", "cd34_d2t1", "cd34_d3t1"),
  project = "SCB-DS1"
)

# Loop over each Seurat object
for (name in names(seurat_list)) {
  seurat_object <- seurat_list[[name]]
  
  # Calculate log10GenesPerUMI
  seurat_object$log10GenesPerUMI <- log10(seurat_object$nFeature_RNA) / log10(seurat_object$nCount_RNA)
  
  # Calculate percent.mt using the PercentageFeatureSet function
  seurat_object[["percent.mt"]] <- PercentageFeatureSet(object = seurat_object, pattern = "^MT-|^mt-|^mt:|^Mito-")
  
  # Calculate mitoRatio
  seurat_object$mitoRatio <- seurat_object$percent.mt / 100
  
  # Save the modified object back to the list
  seurat_list[[name]] <- seurat_object
}

# The seurat_list now contains Seurat objects with the new attributes added.

Mitochondrial Content

grep("^MT-|^mt-|^mt:|^Mito-", rownames(merged_seurat), value = TRUE)

## character(0)

Troubleshooting Missing Mitochondrial Sequences

After testing multiple patterns (e.g., ^MT-, ^mt-, ^mt:, ^Mito-), we found no trace of mitochondrial sequences in our dataset. Here are possible reasons and steps to troubleshoot:

1. Dataset-Specific Naming Convention

2. Lack of Mitochondrial Genes in the Dataset

   • In some cases, mitochondrial genes might not be included in the dataset, possibly due to upstream filtering.

3. Custom or Non-Standard Gene Annotation

   • The dataset might use a custom annotation, such as gene IDs rather than standard symbols.

4. Low or No Mitochondrial Gene Expression

   • Mitochondrial gene expression might be low or absent, especially if sequencing depth is low or if certain cell types rely on glycolysis rather than mitochondrial respiration.

5. Pattern Matching Issue Unexpected formatting in gene names (e.g., trailing spaces) may interfere with pattern matching.

6. Dataset-Specific Filtering of Low-Expression Genes

Some datasets exclude genes with low expression. If mitochondrial genes are filtered out, they will not be present.

We will assume there are no mitochondrial sequences are not available, as there is not other documentations are available for the dataset.

Otherwise, the above code would have been used to find mitochondrial gene percentage and mitoratio. But for our project, mitochondrial filtering will not be necessary.

4.1.3.1.2. Cell Counts

# View the distribution of cell counts
table(merged_seurat$Sample) %>%
    barplot(main="Cell Counts per Sample",
            xlab="Sample", ylab="Cell Count", col="steelblue", las=2)

4.1.3.1.2. UMI Counts per Cell

# Visualize the distribution of UMI counts per cell using orig.ident_mod
VlnPlot(merged_seurat, features = "nCount_RNA", group.by = "Sample", pt.size = 0.1) +
    ggtitle("UMI Counts per Cell") +
    ylab("UMI Count") +
    theme_minimal()

Observation: Based on this violin plot of UMI counts per cell, it looks like the majority of cells have UMI counts clustered below 5,000 (except for CD34_D2T1). However, there are also cells with UMI counts exceeding 7,500 or even 10,000, which might represent high-expressing or possibly doublet cells.

4.1.3.1.3. Genes Detected per Cell

# Visualize the distribution of genes detected per cell using orig.ident_mod
VlnPlot(merged_seurat, features = "nFeature_RNA", group.by = "Sample", pt.size = 0.1) +
    ggtitle("Genes Detected per Cell") +
    ylab("Gene Count") +
    theme_minimal()

Most cells have gene counts around 1,000 to 2,500, with a few outliers reaching above 3,500 or even 4,000. However, we have to be careful of CD34_D2T1 sample which shows more distributed pattern.

4.1.3.1.4. UMIs vs. Genes Detected

# Scatter plot of UMIs vs. genes detected per cell using orig.ident_mod
FeatureScatter(merged_seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "Sample") +
    ggtitle("UMIs vs. Genes Detected") +
    xlab("UMI Count") +
    ylab("Gene Count") +
    theme_minimal()

This scatter plot of UMIs vs. genes detected per cell shows a strong positive correlation between UMI counts and the number of genes detected, which is typical in high-quality single-cell RNA-seq data.

4.1.4. Choosing Cut-off parameters

We have used the following MAD based method for filtering our dataset previously. The stringent filtration resulted in loss of marker data and subsequent failure of manually annotating the cell types where we only obtained four cell types. We will keep the previous filtration procedure here for documentation.

nFeature_RNA_values <- merged_seurat@meta.data$nFeature_RNA
nCount_RNA_values <- merged_seurat@meta.data$nCount_RNA

# For nFeature_RNA

# Median and MAD
nFeature_median <- median(nFeature_RNA_values)
nFeature_mad <- mad(nFeature_RNA_values)

# Mean and Standard Deviation
nFeature_mean <- mean(nFeature_RNA_values)
nFeature_sd <- sd(nFeature_RNA_values)

# Determining min/max threshold

k <- 6

# For nFeature_RNA using median and MAD
nFeature_lower <- nFeature_median - k * nFeature_mad

# Ensure lower thresholds are not negative
nFeature_lower <- max(nFeature_lower, min(nFeature_RNA_values))

cat("lower threshold for Genes Detected:", nFeature_lower, "\n")

# Calculate Median and MAD for nCount_RNA
nCount_median <- median(nCount_RNA_values)
nCount_mad <- mad(nCount_RNA_values)

# Calculate Mean and Standard Deviation for nCount_RNA
nCount_mean <- mean(nCount_RNA_values)
nCount_sd <- sd(nCount_RNA_values)

# Choosing k = 3 for thresholding
k <- 6

# Determine min/max threshold for nCount_RNA using median and MAD
nCount_lower <- nCount_median - k * nCount_mad
nCount_upper <- nCount_median + k * nCount_mad

# Ensure lower threshold is not negative
nCount_lower <- max(nCount_lower, min(nCount_RNA_values))

cat("lower threshold for UMI:", nCount_lower, "\n")
cat("upper threshold for UMI:", nCount_upper, "\n")

4.1.5. Apply Filtering

The code defines a function, filter_seurat, that filters cells in a Seurat object (sobj) based on thresholds derived from count and feature distributions. The overview of the code is adapted from this discussion in Seurat’s Github page Here’s a summary of each step:

1. Define Count Thresholds:
   • If the minimum RNA count per cell (nCount_RNA) is greater than or equal to a specified minimum (minCov, defaulting to 1000), countLOW is set to this minimum count.
   • Otherwise, countLOW is set to the 1st percentile of nCount_RNA, capturing the lower limit for filtering.
   • countHIGH is set to the 99th percentile of nCount_RNA to exclude high outliers.
2. Define Feature Count Threshold:
   • featureLOW is set to the 1st percentile of detected gene features (nFeature_RNA), filtering out cells with low feature counts.
3. Filter Cells:
   • The function removes cells with nFeature_RNA and nCount_RNA values outside the specified range and cells with a mitochondrial gene percentage (percent.mt) of 20% or more (In this particular project, mitochondrial gene percentage based filtration is unnecessary, as the dataset lacks mitochondrial sequence data)
4. Apply the Function:
   • The function is applied to each Seurat object in seurat_list, creating a list of filtered Seurat objects (seurat_list_filtered).

This function helps clean up Seurat objects by excluding cells with extreme counts, low feature counts, or high mitochondrial content.

# Define the filtering function with added print statements
filter_seurat <- function(sobj, minCov = 1000) {
  # Determine countLOW based on coverage distribution
  if (min(sobj$nCount_RNA) >= minCov) {
    countLOW <- min(sobj$nCount_RNA)
  } else {
    countLOW <- quantile(sobj$nCount_RNA, prob = 0.01)
  }
  
  # Set countHIGH to the 99th percentile to exclude high outliers
  countHIGH <- quantile(sobj$nCount_RNA, prob = 0.95)
  
  # Set featureLOW to the 1st percentile to exclude low outliers in feature count
  featureLOW <- quantile(sobj$nFeature_RNA, prob = 0.05)
  
  # Print thresholds for this Seurat object using project name as identifier
  cat("Thresholds for", sobj@project.name, ":\n")
  cat(" - countLOW:", countLOW, "\n")
  cat(" - countHIGH:", countHIGH, "\n")
  cat(" - featureLOW:", featureLOW, "\n")
  
  # Apply filtering based on calculated thresholds and mt percent criteria
  sobj_filtered <- subset(sobj, 
                          subset = nFeature_RNA > featureLOW & 
                                   nCount_RNA > countLOW & 
                                   nCount_RNA < countHIGH & 
                                   percent.mt < 20)
  
  # Print resulting cell count after filtering
  cat(" - Resulting cell count:", ncol(sobj_filtered), "\n\n")
  
  return(sobj_filtered)
}

# Apply the function to each Seurat object in the list
seurat_list_filtered <- lapply(seurat_list, filter_seurat)

## Thresholds for BMMC_D1T1 :
##  - countLOW: 1008 
##  - countHIGH: 7337.4 
##  - featureLOW: 840 
##  - Resulting cell count: 5645 
## 
## Thresholds for BMMC_D1T2 :
##  - countLOW: 1012 
##  - countHIGH: 7326.15 
##  - featureLOW: 830 
##  - Resulting cell count: 5698 
## 
## Thresholds for CD34_D2T1 :
##  - countLOW: 1002 
##  - countHIGH: 8788.55 
##  - featureLOW: 778 
##  - Resulting cell count: 2179 
## 
## Thresholds for CD34_D3T1 :
##  - countLOW: 1002 
##  - countHIGH: 6310.45 
##  - featureLOW: 755 
##  - Resulting cell count: 5171

4.1.5.1. Compare before and after filtration

# Load required libraries
library(Seurat)
library(ggplot2)
library(patchwork)

# Define a function to create violin plots for comparison
plot_violin_comparison <- function(sobj_before, sobj_after, sobj_name) {
  # Define common theme elements for consistent styling
  common_theme <- theme(
    plot.title = element_text(hjust = 0.5, size = 14),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.position = "none",
    plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm")
  )
  
  # Create feature plot (before)
  plot_before <- VlnPlot(sobj_before, 
                        features = c("nFeature_RNA", "nCount_RNA"),
                        pt.size = 0.1, 
                        ncol = 2) &
    common_theme &
    ggtitle(paste0(sobj_name, "\nBefore Filtering"))
  
  # Create feature plot (after)
  plot_after <- VlnPlot(sobj_after, 
                        features = c("nFeature_RNA", "nCount_RNA"),
                        pt.size = 0.1, 
                        ncol = 2) &
    common_theme &
    ggtitle(paste0(sobj_name, "\nAfter Filtering"))
  
  # Stack plots vertically with consistent spacing
  combined_plot <- (plot_before / plot_after) +
    plot_layout(heights = c(1, 1))
  
  return(combined_plot)
}

# Generate comparison plots for each Seurat object
sample_names <- c("BMMC_d1t1", "BMMC_d1t2", "CD34_d2t1", "CD34_d3t1")
comparison_plots <- mapply(
  function(sobj, sobj_filtered, sobj_name) {
    plot_violin_comparison(sobj, sobj_filtered, sobj_name)
  },
  seurat_list,
  seurat_list_filtered,
  sample_names,
  SIMPLIFY = FALSE
)

# Combine all plots in a grid with consistent spacing
final_plot <- wrap_plots(comparison_plots, 
                        ncol = 2,  # Two columns for balanced layout
                        guides = "collect") +
  plot_layout(guides = "collect") &
  theme(plot.margin = margin(20, 20, 20, 20))

# Display the final plot
final_plot

4.1.6. Doublet removal with DoubletFinder

4.1.6.1. Explanation of why we are performing doublet removal

Purpose: Remove potential doublets to ensure that each observation represents a single cell.

4.1.6.2. Normalization

Which Normalization method is used by the Seurat Normalization function by default?

– Answer: LogNormalize

4.1.6.3. Feature Selection

When using DoubletFinder for detecting doublets in single-cell RNA sequencing data, it is generally recommended to select between 2000 and 3000 features. This range is considered optimal for capturing the most informative genes while minimizing noise, which is crucial for accurate doublet detection.

4.1.6.4. Scale

4.1.6.5. Run PCA

We chose to run 50 PC to determine the elbow point.

# Apply preprocessing and PCA to each Seurat object
seurat_list_filtered <- lapply(names(seurat_list_filtered), function(sample_name) {
  seurat_obj <- seurat_list_filtered[[sample_name]]
  
  # Step 1: Normalization
  seurat_obj <- NormalizeData(seurat_obj)
  
  # Step 2: Feature Selection
  seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 3000)
  
  # Step 3: Scaling
  seurat_obj <- ScaleData(seurat_obj, features = VariableFeatures(object = seurat_obj))
  
  # Step 4: Run PCA
  seurat_obj <- RunPCA(
    seurat_obj,
    npcs = 50,        # Number of PCs to compute
    verbose = FALSE    # Suppress verbose output
  )
  
  #print a message
  message(paste("Preprocessing and PCA completed for", sample_name))
  
  return(seurat_obj)
})

# Fix Seurat Object names

sample_names <- c("bmmc_d1t1", "bmmc_d1t2", "cd34_d2t1", "cd34_d3t1")

# Ensure the number of names matches the number of Seurat objects
if (length(sample_names) == length(seurat_list_filtered)) {
  names(seurat_list_filtered) <- sample_names
  message("Specific sample names assigned to seurat_list_filtered.")
} else {
  warning("Length of sample_names does not match the number of Seurat objects. Assigning default names instead.")
  names(seurat_list_filtered) <- paste0("Sample_", seq_along(seurat_list_filtered))
}

4.1.6.6. Elbow Point

The Elbow Plot is commonly used to help decide the number of principal components (PCs) to retain. By plotting the variance explained by each PC, we look for an “elbow point” where the added value of each subsequent PC diminishes significantly.

4.1.6.7. Choosing the Number of PCs

To determine the optimal number of PCs to retain, we consider both the traditional elbow method and the cumulative variance explained. This ensures that a significant portion of the data’s variability is captured while maintaining computational efficiency. Ideally we should have chosen 80% to 90% of cumulative variance explained to determine the elbow points but cost and time of computation increases exponentially. So, we chose 75% cumulative variance explained for the sake of computational efficiency.

# Initialize lists to store cumulative variance data and elbow points
cumulative_variance_list <- list()
elbow_points <- data.frame(
  Sample = character(),
  ElbowPoint = numeric(),
  stringsAsFactors = FALSE
)

# Set the threshold for cumulative variance explained.
threshold <- 75

# Loop over each Seurat object in seurat_list_filtered
for (i in seq_along(seurat_list_filtered)) {
  seurat_obj <- seurat_list_filtered[[i]]
  sample_name <- names(seurat_list_filtered)[i]
  
  # Check if PCA was successfully run
  if (!"pca" %in% names(seurat_obj@reductions) ||
      is.null(seurat_obj[["pca"]]@stdev)) {
    warning(paste("PCA has not been run on", sample_name, ". Skipping variance calculation."))
    next  # Skip to the next iteration
  }
  
  # Get the standard deviation of each PC and calculate the variance explained
  stdev <- seurat_obj[["pca"]]@stdev
  variance_explained <- (stdev^2) / sum(stdev^2) * 100
  
  # Compute cumulative variance explained
  cumulative_variance_explained <- cumsum(variance_explained)
  
  # Create a data frame to store cumulative variance explained for each PC
  pca_cumulative_variance_df <- data.frame(
    PC = 1:length(cumulative_variance_explained),
    CumulativeVarianceExplained = cumulative_variance_explained,
    Sample = sample_name
  )
  
  # Determine the elbow point where cumulative variance exceeds the threshold
  elbow_point <- which(cumulative_variance_explained >= threshold)[1]
  
  if (is.na(elbow_point)) {
    warning(paste("Threshold of", threshold, "% not reached in", sample_name))
    elbow_point <- NA
  }
  
  # Store the data frame in the list
  cumulative_variance_list[[i]] <- pca_cumulative_variance_df
  
  # Record the elbow point
  elbow_points <- rbind(elbow_points, data.frame(Sample = sample_name, ElbowPoint = elbow_point))
}

# Combine all cumulative variance data frames into one
combined_variance_df <- do.call(rbind, cumulative_variance_list)

# Remove any rows with NA elbow points
elbow_points_clean <- na.omit(elbow_points)

# Plot cumulative variance explained for all samples
p <- ggplot(combined_variance_df, aes(x = PC, y = CumulativeVarianceExplained, color = Sample)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  geom_hline(yintercept = threshold, color = "red", linetype = "dashed") +
  # Add vertical dashed lines for elbow points
  geom_vline(data = elbow_points_clean, aes(xintercept = ElbowPoint, color = Sample), linetype = "dashed") +
  # Annotate elbow points
  geom_text(
    data = elbow_points_clean, 
    aes(x = ElbowPoint, y = threshold + 2, label = paste("PC", ElbowPoint), color = Sample),
    angle = 90, vjust = -0.5, hjust = 0, size = 3, show.legend = FALSE
  ) +
  labs(
    title = "Cumulative Variance Explained by Principal Components",
    x = "Principal Component",
    y = "Cumulative Variance Explained (%)"
  ) +
  theme_minimal() +
  theme(
    legend.title = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold")
  )

# Display the plot
print(p)

# print the elbow points for each sample
print(elbow_points_clean)

##      Sample ElbowPoint
## 1 bmmc_d1t1         11
## 2 bmmc_d1t2         11
## 3 cd34_d2t1         19
## 4 cd34_d3t1         17

# Define the DoubletFinder processing function
run_doublet_finder <- function(seurat_obj, sample_name, elbow_point, 
                                pN = 0.25, doublet_rate = 0.075, sct = FALSE) {
  cat("\nProcessing Sample:", sample_name, "\n")
  
  # Step 1: Parameter Sweep to Estimate Optimal pK
  cat("Performing parameter sweep to estimate optimal pK...\n")
  tryCatch({
    sweep.res.list <- paramSweep_v3(seurat_obj, PCs = 1:elbow_point, sct = sct)
    sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
    bcmvn <- find.pK(sweep.stats)
  }, error = function(e) {
    stop(paste("Error during parameter sweep for", sample_name, ":", e$message))
  })
  
  # Validate bcmvn data frame
  if (!all(c("pK", "BCmetric") %in% colnames(bcmvn))) {
    stop(paste("Expected columns 'pK' and 'BCmetric' not found in sweep results for", sample_name))
  }
  
  # Step 2: Select Optimal pK (pK with maximum BCmetric)
  optimal_pK <- as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))
  if (is.na(optimal_pK)) {
    stop(paste("Optimal pK could not be determined for", sample_name))
  }
  
  cat("Optimal pK for", sample_name, ":", optimal_pK, "\n")
  
  # Step 3: Plot BCmetric vs pK
  pK_plot <- ggplot(bcmvn, aes(x = as.numeric(as.character(pK)), y = BCmetric)) +
    geom_point(color = "blue") +
    geom_line(color = "blue") +
    geom_vline(xintercept = optimal_pK, color = "red", linetype = "dashed") +
    annotate("text", x = optimal_pK, y = max(bcmvn$BCmetric, na.rm = TRUE), 
             label = paste("Optimal pK:", optimal_pK), hjust = -0.1, vjust = -0.5, color = "red") +
    labs(title = paste("Optimal pK Value Identification for", sample_name),
         x = "pK",
         y = "BCmetric") +
    theme_minimal()
  
  print(pK_plot)
  
  # Step 4: Estimate Expected Number of Doublets
  nExp <- round(doublet_rate * ncol(seurat_obj))
  cat("Estimated number of doublets (nExp) for", sample_name, ":", nExp, "\n")
  
  # Step 5: Run DoubletFinder
  cat("Running DoubletFinder...\n")
  tryCatch({
    seurat_obj <- doubletFinder_v3(seurat_obj, 
                                   PCs = 1:elbow_point, 
                                   pN = pN, 
                                   pK = optimal_pK, 
                                   nExp = nExp, 
                                   reuse.pANN = FALSE, 
                                   sct = sct)
  }, error = function(e) {
    stop(paste("Error during DoubletFinder execution for", sample_name, ":", e$message))
  })
  
  # Step 6: Identify the Classification Column
  classification_col <- paste0("DF.classifications_", pN, "_", optimal_pK, "_", nExp)
  
  # Verify the classification column exists
  if (!(classification_col %in% colnames(seurat_obj@meta.data))) {
    stop(paste("Classification column", classification_col, "not found in metadata for", sample_name))
  }
  
  # Display the classification counts
  cat("Doublet Classification Counts for", sample_name, ":\n")
  print(table(seurat_obj@meta.data[[classification_col]]))
  
  # Step 7: Filter Out Doublets (Keep Only Singlets)
  cat("Filtering out doublets, retaining only singlets...\n")
  seurat_obj_filtered <- subset(seurat_obj, subset = seurat_obj@meta.data[[classification_col]] == "Singlet")
  
  # Confirm the number of cells after filtering
  cat("Number of cells after doublet removal for", sample_name, ":", ncol(seurat_obj_filtered), "\n")
  
  # Return the filtered Seurat object, optimal_pK, and the pK plot
  return(list(
    seurat_obj = seurat_obj_filtered,
    optimal_pK = optimal_pK,
    pK_plot = pK_plot
  ))
}

# Initialize lists to store results
filtered_seurat_list <- list()
optimal_pK_values <- data.frame(
  Sample = character(),
  Optimal_pK = numeric(),
  stringsAsFactors = FALSE
)

# Ensure that 'seurat_list_filtered' is named correctly
# This step is crucial to match samples with their elbow points
names(seurat_list_filtered) <- elbow_points_clean$Sample

# Iterate over each Seurat object and apply the DoubletFinder function
for (sample_name in names(seurat_list_filtered)) {
  seurat_obj <- seurat_list_filtered[[sample_name]]
  
  # Retrieve the elbow_point for this sample from elbow_points_clean
  elbow_point <- elbow_points_clean$ElbowPoint[elbow_points_clean$Sample == sample_name]
  
  # Check if elbow_point is available
  if (length(elbow_point) == 0 || is.na(elbow_point)) {
    warning(paste("Elbow point not found for", sample_name, ". Skipping DoubletFinder."))
    next  # Skip to the next iteration
  }
  
  # Run the DoubletFinder function
  result <- tryCatch({
    run_doublet_finder(seurat_obj, sample_name, elbow_point, 
                       pN = 0.25, doublet_rate = 0.075, sct = FALSE)
  }, error = function(e) {
    warning(paste("DoubletFinder failed for", sample_name, ":", e$message))
    return(NULL)
  })
  
  # If the result is NULL due to an error, skip storing
  if (is.null(result)) {
    next
  }
  
  # Store the filtered Seurat object
  filtered_seurat_list[[sample_name]] <- result$seurat_obj
  
  # Record the optimal pK value
  optimal_pK_values <- rbind(optimal_pK_values, 
                             data.frame(Sample = sample_name, Optimal_pK = result$optimal_pK))
}

## 
## Processing Sample: bmmc_d1t1 
## Performing parameter sweep to estimate optimal pK...
## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## Optimal pK for bmmc_d1t1 : 0.005

## Estimated number of doublets (nExp) for bmmc_d1t1 : 423 
## Running DoubletFinder...
## [1] "Creating 1882 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
## Doublet Classification Counts for bmmc_d1t1 :
## 
## Doublet Singlet 
##     423    5222 
## Filtering out doublets, retaining only singlets...
## 
## Processing Sample: bmmc_d1t2 
## Performing parameter sweep to estimate optimal pK...
## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## Optimal pK for bmmc_d1t2 : 0.005

## Estimated number of doublets (nExp) for bmmc_d1t2 : 427 
## Running DoubletFinder...
## [1] "Creating 1899 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
## Doublet Classification Counts for bmmc_d1t2 :
## 
## Doublet Singlet 
##     427    5271 
## Filtering out doublets, retaining only singlets...
## 
## Processing Sample: cd34_d2t1 
## Performing parameter sweep to estimate optimal pK...
## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## Optimal pK for cd34_d2t1 : 0.12

## Estimated number of doublets (nExp) for cd34_d2t1 : 163 
## Running DoubletFinder...
## [1] "Creating 726 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
## Doublet Classification Counts for cd34_d2t1 :
## 
## Doublet Singlet 
##     163    2016 
## Filtering out doublets, retaining only singlets...
## 
## Processing Sample: cd34_d3t1 
## Performing parameter sweep to estimate optimal pK...
## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## Optimal pK for cd34_d3t1 : 0.005

## Estimated number of doublets (nExp) for cd34_d3t1 : 388 
## Running DoubletFinder...
## [1] "Creating 1724 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
## Doublet Classification Counts for cd34_d3t1 :
## 
## Doublet Singlet 
##     388    4783 
## Filtering out doublets, retaining only singlets...

# Confirm the filtered list
print(filtered_seurat_list)

## list()

Here, we have determined pK values that corresponds to each of the sample’s uniquely identified elbow point.

4.2 Merging and Batch Correction

The batch correction process in Seurat allows the integration of multiple datasets, each potentially generated from different experimental batches or conditions. Batch correction is essential when we want to combine data from separate sources while minimizing batch effects, which can obscure the biological signal of interest. The batch correction process in Seurat aligns multiple datasets from different experimental batches to minimize technical differences. It begins with normalization and identification of highly variable features in each dataset. Shared variable features across datasets are selected as common points for alignment. Using these features, Seurat identifies “anchors,” which are pairs of biologically similar cells across batches. Anchors serve as reference points, allowing Seurat to map datasets onto a shared space. The integration step then adjusts gene expression values across datasets, reducing batch-specific effects while retaining biological signals. Seurat generates an “integrated” assay with batch-corrected data, which replaces the batch-biased individual datasets.

4.2.1 1. Merge All Four Samples Without Batch Correction

# Merge all four samples without batch correction
merged_seurat_no_batch <- merge(
  x = seurat_list_filtered[[1]],
  y = seurat_list_filtered[-1],
  add.cell.ids = names(filtered_seurat_list),
  project = "Merged_NoBatchCorrection"
)

# Verify the merged object
print(merged_seurat_no_batch)

## An object of class Seurat 
## 20287 features across 18693 samples within 1 assay 
## Active assay: RNA (20287 features, 0 variable features)

4.2.2 Merge All Four Samples Using Seurat’s Data Integration Method (Batch Correction)

# Select features for integration
integration_features <- SelectIntegrationFeatures(object.list = seurat_list_filtered, nfeatures = 3000)

# Find integration anchors
anchors <- FindIntegrationAnchors(object.list = seurat_list_filtered, 
                                  anchor.features = integration_features, 
                                  reduction = "rpca", 
                                  dims = 1:50)

# Integrate data
merged_seurat_batch_corrected <- IntegrateData(anchorset = anchors, dims = 1:50)

# Switch to the integrated assay
DefaultAssay(merged_seurat_batch_corrected) <- "integrated"

# Verify the integrated object
# print(merged_seurat_batch_corrected)

4.2.4. Batch effects in single-cell data can arise from various sources, including:

Sequencing Depth: Different batches may have variations in the number of reads per cell.

Library Preparation: Variations in reagent quality or handling can cause batch effects.

Cell Preparation: Differences in cell handling (e.g., sorting, staining) across batches can influence the data.

Instrument Settings: Variability in instrument calibration across runs can introduce batch effects.

# merged_seurat_no_batch <- NormalizeData(merged_seurat_no_batch, normalization.method = "LogNormalize", scale.factor = 10000)
# merged_seurat_batch_corrected <- NormalizeData(merged_seurat_batch_corrected, normalization.method = "LogNormalize", scale.factor = 10000)

4.1.8. Re-Feature Selection

4.1.8.1. Purpose

The purpose of Feature Selection is to identify a subset of genes (features) that exhibit the most significant variability across cells. This variability is often indicative of biological differences, such as distinct cell types, cell states, or responses to different conditions. Selecting these features improves the efficiency and effectiveness of downstream analyses, such as clustering, dimensionality reduction, and identifying cell types.

4.1.8.2. How Are Features Selected?

In Seurat, features are typically selected based on high variability across cells using the FindVariableFeatures() function, which implements several methods for feature selection, such as vst, mean.var.plot, and dispersion.

1. Calculate Mean and Variance of Each Gene

Seurat calculates the average expression (mean) and variability (variance) of each gene across all cells.

2. Standardize Variability

To ensure comparability across genes with different mean expression levels, Seurat uses a method such as Variance Stabilizing Transformation (VST) to adjust the variance based on the mean expression. This transformation helps in selecting genes with true biological variability rather than variability caused by technical noise.

3. Rank Genes by Variability

After transformation, Seurat ranks genes by their standardized variability, identifying those that deviate the most from expected variance at a given mean expression level.

4. Select Top Variable Genes

Based on the user-defined nfeatures parameter (commonly set to 2000-3000), Seurat selects the top genes with the highest variance as “highly variable genes.” This subset is then used for scaling, dimensionality reduction (PCA and UMAP), and clustering.

# Identify highly variable features (genes)
merged_seurat_no_batch <- FindVariableFeatures(
  object = merged_seurat_no_batch, 
  selection.method = "vst",   
  nfeatures = 3000,              
  verbose = FALSE,
  assay = "RNA"
)

merged_seurat_batch_corrected <- FindVariableFeatures(
  object = merged_seurat_batch_corrected, 
  selection.method = "vst",   
  nfeatures = 3000,              
  verbose = FALSE,
  assay = "RNA"
)

5. Dimensionality Reduction

5.1. Dimensionality Reduction

We will perform dimensionality reduction using Principal Component Analysis (PCA) followed by Uniform Manifold Approximation and Projection (UMAP) on both the merged (non-batch corrected) and integrated (batch corrected) datasets.

5.1.1. Scaling the Data

# Scale data
merged_seurat_no_batch <- ScaleData(merged_seurat_no_batch, verbose = FALSE)
merged_seurat_batch_corrected <- ScaleData(merged_seurat_batch_corrected, verbose = FALSE)

5.1.2 Running PCA

# Run PCA
merged_seurat_no_batch <- RunPCA(merged_seurat_no_batch, features = VariableFeatures(object = merged_seurat_no_batch, npcs = 50))
merged_seurat_batch_corrected <- RunPCA(merged_seurat_batch_corrected, verbose = FALSE, npcs = 50)

5.1.3 Determining the Number of Principal Components

# Elbow plot for non-batch corrected data
ElbowPlot(merged_seurat_no_batch, ndims = 50) + ggtitle("Elbow Plot - Non-batch Corrected Data")

# Elbow plot for batch-corrected data
ElbowPlot(merged_seurat_batch_corrected, ndims = 50) + ggtitle("Elbow Plot - Batch Corrected Data")

5.1.3.1 Choosing elbow point (number of dimensions)

The Elbow Plot displays the standard deviation of each PC. We look for a point (the “elbow”) where the rate of decrease sharply changes, indicating that additional PCs contribute less to explaining variance.

Explanation: In the Elbow Plot, we observed that after PC30, the variance explained by each additional PC diminishes significantly. Therefore, we chose to use the first 30 PCs for UMAP and clustering. This balances capturing sufficient variance while reducing noise from less informative PCs.

elbow_point <- 30

5.1.4. Why Use PCA with UMAP

PCA serves as a linear dimensionality reduction technique that helps in denoising the data and capturing the primary axes of variation. However, it may not capture complex, nonlinear relationships in the data. UMAP, on the other hand, is a nonlinear dimensionality reduction method that preserves both local and global data structure.

By combining PCA with UMAP:

Noise Reduction: PCA filters out noise by focusing on components that explain the most variance.

Computational Efficiency: Reducing dimensions with PCA before UMAP speeds up the computation.

Capturing Nonlinear Structure: UMAP effectively captures complex patterns and relationships in the data that PCA might miss.

Using both methods results in a more meaningful clustering and visualization.

5.2. Clustering

We performed clustering on the reduced data to identify groups of similar cells.

library(Matrix)
library(igraph)

# For batch-corrected data
merged_seurat_batch_corrected <- FindNeighbors(merged_seurat_batch_corrected, dims = 1:elbow_point)
merged_seurat_batch_corrected <- FindClusters(merged_seurat_batch_corrected, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 18693
## Number of edges: 722266
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9393
## Number of communities: 23
## Elapsed time: 2 seconds

# For non-batch corrected data
merged_seurat_no_batch <- FindNeighbors(merged_seurat_no_batch, dims = 1:elbow_point)
merged_seurat_no_batch <- FindClusters(merged_seurat_no_batch, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 18693
## Number of edges: 689256
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9411
## Number of communities: 24
## Elapsed time: 2 seconds

The resolution parameter controls the granularity of the clustering:

• Lower values lead to fewer clusters.
• Higher values result in more clusters.

We could have adjusted the resolution < 0.21 to achieve between 7-15 clusters but decided to use 0.5 value instead for better resolution (and thus more clusters).

# For batch-corrected data
num_clusters_batch_corrected <- length(unique(merged_seurat_batch_corrected$seurat_clusters))
cat("Number of clusters in batch-corrected data:", num_clusters_batch_corrected, "\n")

## Number of clusters in batch-corrected data: 23

# For non-batch corrected data
num_clusters_no_batch <- length(unique(merged_seurat_no_batch$seurat_clusters))
cat("Number of clusters in non-batch corrected data:", num_clusters_no_batch, "\n")

## Number of clusters in non-batch corrected data: 24

5.2.2. Run UMAP

# Run UMAP based on the PCA results
merged_seurat_no_batch <- RunUMAP(merged_seurat_no_batch, dims = 1:elbow_point)

merged_seurat_batch_corrected <- RunUMAP(merged_seurat_batch_corrected, dims = 1:elbow_point)

5.2.1. Visualization of Clusters

We visualized the clustering results using UMAP plots.

# Generate UMAP plots for batch-corrected and non-batch corrected data
p1 <- DimPlot(merged_seurat_batch_corrected, reduction = "umap", group.by = "seurat_clusters") + ggtitle("Batch-Corrected Data")
p2 <- DimPlot(merged_seurat_no_batch, reduction = "umap", group.by = "seurat_clusters") + ggtitle("Non-Batch Corrected Data")

# Arrange the plots
p1

p2

Explanation: The UMAP plots display cells in two-dimensional space, colored by their assigned clusters. The clustering results show distinct groups of cells, indicating successful identification of cell populations.

Was batch correction necessary:

Batch correction improves clustering by:

• Aligning Similar Cell Types across samples, reducing batch-related variability.
• Enhancing Cluster Compactness for clearer separation, making clusters more distinct.
• Reducing Artifacts like stretched or fragmented clusters, improving biological interpretation.

4.2.5. Batch correction is crucial

The results could be biased by technical artifacts, misrepresenting true biological variability.

Downstream analyses, such as clustering and differential expression, may be inaccurate if technical differences are mistaken for biological differences.

To demonstrate the effect of each parameter, we could create additional visualizations for each suspected source of batch effect (e.g., comparing UMAPs colored by sequencing depth or library preparation batch). However, the primary goal is to achieve a plot where cells cluster by type rather than batch, as seen in the “With Batch Correction” UMAP plot on the right.

Parameters that need to be looked at for batch correction:

Number of Unique Molecular Identifiers (nUMI)
Number of Detected Genes (nGene)
Percentage of Mitochondrial Genes (percent.mt)
Batch Labels or Experimental Conditions
Highly Variable Features (HVFs)
Dimensions for Principal Component Analysis (PCA)
Clustering Resolution
Integration Anchors for Batch Correction

5.2.2. Keeping Results for Subsequent Tasks

The results of the dimensionality reduction and clustering are stored in the Seurat objects (merged_seurat_batch_corrected and merged_seurat_no_batch). These objects retain the PCA, UMAP embeddings, and cluster assignments, which will be used for downstream analyses such as marker gene identification and differential expression.

Part 3 — Cell-Type Annotation, Differential Expression & Pathways

6. Cell Type Annotation

6.1. Automatic Cell Type Annotation with SingleR

In this section, we use SingleR, a tool for automatic cell type annotation, to assign cell type labels to our integrated Seurat object. We utilize the built-in reference dataset HumanPrimaryCellAtlasData from the celldex package. Finally, we visualize the annotation results on a UMAP plot.

# Load the Human Primary Cell Atlas reference data
reference <- readRDS("../data/celldex_annot.RDS")

# Extract the normalized expression data from the integrated assay
# SingleR requires a normalized expression matrix
# We use the "integrated" assay's data slot
expression_data <- GetAssayData(merged_seurat_batch_corrected, assay = "integrated", slot = "data")

# Run SingleR to annotate cell types
# 'labels' parameter uses the main cell type annotations from the reference
singleR_results <- SingleR(test = expression_data, 
                           ref = reference, 
                           labels = reference$label.main, 
                           assay.type.test = "logcounts", 
                           assay.type.ref = "logcounts")

# Add SingleR labels to the Seurat object's metadata
merged_seurat_batch_corrected$SingleR.labels <- singleR_results$labels

# Verify the annotation results
table(merged_seurat_batch_corrected$SingleR.labels)

## 
##            B_cell                BM        BM & Prog.               CMP 
##              1315                 7                91              1069 
##                DC      Erythroblast               GMP        HSC_-G-CSF 
##                 4               341              1688               104 
##               MEP          Monocyte         Myelocyte       Neutrophils 
##               374              3346               144                29 
##           NK_cell  Pre-B_cell_CD34-  Pro-B_cell_CD34+     Pro-Myelocyte 
##               897              1523              1174               648 
##           T_cells Tissue_stem_cells 
##              5937                 2

# Plot UMAP with SingleR cell type labels
DimPlot(merged_seurat_batch_corrected, 
        reduction = "umap", 
        group.by = "SingleR.labels", 
        label = TRUE, 
        label.size = 2,   
        repel = TRUE) +
  ggtitle("UMAP Plot with SingleR Cell Type Annotations") +
  theme_minimal()

6.2. Manual Annotation

6.2.1. Perform Differential Expression Analysis for Cell-Type Annotation

# Set the default assay
DefaultAssay(merged_seurat_batch_corrected) <- "integrated"

# Perform differential expression analysis
all_markers <- FindAllMarkers(
  object = merged_seurat_batch_corrected,
  assay = "integrated",
  only.pos = FALSE,
  min.pct = 0.1,
  logfc.threshold = 0.25, # Default 0.25, could used lesser value to detect weaker signals,
  random.seed = 789
)

# View the top markers for each cluster
print(head(all_markers))

##         p_val avg_log2FC pct.1 pct.2 p_val_adj cluster    gene
## S100A9      0   4.685399 1.000 0.660         0       0  S100A9
## S100A12     0   4.563143 0.960 0.519         0       0 S100A12
## S100A8      0   4.262029 0.999 0.665         0       0  S100A8
## VCAN        0   3.837171 0.978 0.618         0       0    VCAN
## FCN1        0   3.832863 0.990 0.571         0       0    FCN1
## CD14        0   3.401156 0.913 0.473         0       0    CD14

6.2.2. Use Markers from Table 2 to Identify Cell Types and Assign Names to Clusters

6.2.2.1. Define Cell Type Markers

We used acronyms for the cell types here. Please refer to the assignment sheet for full names.

# Define a list of marker genes for each cell type
cell_type_markers <- list(
  HSC = c("CD34", "CD38", "Sca1", "Kit"),
  LMPP = c("CD38", "CD52", "CSF3R", "ca1", "Kit", "CD34", "Flk2"),
  CLP = c("IL7R"),
  GMP_Neutrophils = c("ELANE"),
  CMP = c("IL3", "GM-CSF", "M-CSF"),
  B = c("CD19"),
  Pre_B = c("CD19", "CD34"),
  Plasma = c("SDC1", "IGHA1", "IGLC1", "MZB1", "JCHAIN"),
  CD8_T = c("CD3D", "CD3E", "CD8A", "CD8B"),
  CD4_T = c("CD3D", "CD3E", "CD4"),
  NK = c("FCGR3A", "NCAM1", "NKG7", "KLRB1"),
  Erythrocytes = c("GATA1", "HBB", "HBA1", "HBA2"),
  pDC = c("IRF8", "IRF4", "IRF7"),
  cDC = c("CD1C", "CD207", "ITGAM", "NOTCH2", "SIRPA"),
  CD14_Monocytes = c("CD14", "CCL3", "CCL4", "IL1B"),
  CD16_Monocytes = c("FCGR3A", "CD68", "S100A12"),
  Basophils = c("GATA2")
)

6.2.2.4. Missing Markers in the datasets

missing_genes_merged_seurat <- setdiff(unique(unlist(cell_type_markers)), rownames(merged_seurat))
cat("Missing gene markers in unfiltered merged seurat object:", missing_genes_merged_seurat)

## Missing gene markers in unfiltered merged seurat object: Sca1 Kit ca1 Flk2 GM-CSF M-CSF IGHA1 IGLC1 JCHAIN

missing_genes_merged_seurat_no_batch <- setdiff(unique(unlist(cell_type_markers)), rownames(merged_seurat_no_batch))
cat("Missing gene markers in filtered merged seurat object (non-batch corrected):", missing_genes_merged_seurat)

## Missing gene markers in filtered merged seurat object (non-batch corrected): Sca1 Kit ca1 Flk2 GM-CSF M-CSF IGHA1 IGLC1 JCHAIN

missing_genes_merged_seurat_batch_corrected <- setdiff(unique(unlist(cell_type_markers)), rownames(merged_seurat_batch_corrected[["integrated"]]))
cat("Missing gene markers in filtered merged seurat object (batch-corrected):", missing_genes_merged_seurat_batch_corrected)

## Missing gene markers in filtered merged seurat object (batch-corrected): Sca1 Kit ca1 Flk2 IL3 GM-CSF M-CSF IGHA1 IGLC1 JCHAIN CD207

It looks like there are two marker genes that were lost in the batch correction procedure.

unique_missing_in_merged <- setdiff(missing_genes_merged_seurat_no_batch, missing_genes_merged_seurat_batch_corrected)
unique_missing_in_integrated <- setdiff(missing_genes_merged_seurat_batch_corrected, missing_genes_merged_seurat_no_batch)


cat("Number of genes uniquely missing after batch correction:", length(unique_missing_in_integrated), "\n")

## Number of genes uniquely missing after batch correction: 2

print(unique_missing_in_integrated)

## [1] "IL3"   "CD207"

Strategy for Manual cell annotation

Calculate Module Scores Using AddModuleScore

# Ensure all marker genes are present in the dataset
features_list <- lapply(cell_type_markers, function(markers) {
  markers[markers %in% rownames(merged_seurat_batch_corrected)]
})

# Remove empty gene sets
features_list <- features_list[sapply(features_list, length) > 0]

# Add module scores for each cell type
merged_seurat_batch_corrected <- AddModuleScore(
  object = merged_seurat_batch_corrected,
  features = features_list,
  name = "CellTypeScore"
)

# Get the names of the module score columns
score_column_indices <- grep("^CellTypeScore", colnames(merged_seurat_batch_corrected@meta.data))
score_column_names <- paste0( names(features_list))
colnames(merged_seurat_batch_corrected@meta.data)[score_column_indices] <- score_column_names

Assign Cell Types Based on Highest Module Score

# Extract module scores
module_scores <- merged_seurat_batch_corrected@meta.data[, score_column_names]

# Assign cell types to cells based on the highest module score
merged_seurat_batch_corrected$AssignedCellType <- apply(
  module_scores, 1, function(x) names(x)[which.max(x)]
)

3. Visualize the UMAP with Assigned Cell Types

# Plot UMAP with assigned cell types and adjusted label size
DimPlot(
  merged_seurat_batch_corrected, 
  reduction = "umap", 
  group.by = "AssignedCellType", 
  label = TRUE, 
  repel = TRUE,
  label.size = 2  
) + 
  ggtitle("UMAP Plot with Assigned Cell Types") + 
  theme_minimal()

6.2.3 Automatic vs Manual Annotation Comparison

6.2.3.1. Confusion Matrix

confusion_matrix <- table(
  Manual = merged_seurat_batch_corrected$AssignedCellType,
  SingleR = merged_seurat_batch_corrected$SingleR.labels
)

library(pheatmap)
confusion_matrix_numeric <- as.matrix(confusion_matrix)

# Plot the heatmap
pheatmap(
  confusion_matrix_numeric,
  display_numbers = TRUE,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  fontsize_number = 3,
  main = "Confusion Matrix: Manual vs SingleR Annotations"
)

6.2.3.2. Calculate Agreement Metrics

library(caret)

# Prepare the labels
manual_labels <- as.factor(merged_seurat_batch_corrected$AssignedCellType)
singleR_labels <- as.factor(merged_seurat_batch_corrected$SingleR.labels)

# Ensure both factors have the same levels
levels(singleR_labels) <- union(levels(singleR_labels), levels(manual_labels))
levels(manual_labels) <- union(levels(singleR_labels), levels(manual_labels))

# Compute the confusion matrix and statistics
confusion_stats <- confusionMatrix(singleR_labels, manual_labels)

# Extract summary statistics
summary_stats <- list(
  Accuracy = confusion_stats$overall["Accuracy"],
  Kappa = confusion_stats$overall["Kappa"],
  Precision = confusion_stats$byClass["Pos Pred Value"],  # Positive Predictive Value
  Recall = confusion_stats$byClass["Sensitivity"],
  F1_Score = confusion_stats$byClass["F1"]
)

# Print the summary
print(summary_stats)

## $Accuracy
##  Accuracy 
## 0.1030867 
## 
## $Kappa
##      Kappa 
## 0.05672409 
## 
## $Precision
## [1] NA
## 
## $Recall
## [1] NA
## 
## $F1_Score
## [1] NA

6.2.3.3. Observation

This confusion matrix and associated statistics compare two sets of cell type annotations: manual annotations (Reference) and SingleR predictions (Prediction). Here’s what the data and metrics indicate:

’Accuracy (0.1031): The proportion of correctly classified cells (e.g., both predicted and manually labeled as the same cell type). Kappa (0.0567): Measures agreement between the two annotations, adjusted for chance. A low value indicates weak agreement. No Information Rate (NIR, 0.1675): The accuracy expected by randomly guessing the most frequent class.

The low accuracy and agreement suggest discrepancies in definitions or marker usage between the manual annotations and the SingleR reference database.

6.2.6 Visualize Gene Expression of Marker Genes

Let’s choose CD3D (T cell marker), ELANE (GMP_Neutrophils), and NKG7 (NK cell marker).

6.2.2.7. Violin Plots

# Violin plots for selected genes
VlnPlot(merged_seurat_batch_corrected, 
        features = c("CD3D", "ELANE", "NKG7"), 
        group.by = "AssignedCellType", 
        pt.size = 0) +
  theme_minimal()

6.2.2.8. Feature Plots (UMAP)

# Feature plots for selected genes
FeaturePlot(merged_seurat_batch_corrected, 
            features = c("CD3D", "ELANE", "NKG7"), 
            reduction = "umap", 
            cols = c("lightgrey", "blue")) +
  theme_minimal()

6.3. Cell-Type Proportions (Bonus)

6.3.1. Compute Cell-Type Proportions for Each Sample

# Compute cell-type proportions
cell_type_counts <- table(merged_seurat_batch_corrected$AssignedCellType, merged_seurat_batch_corrected$Sample)
cell_type_proportions <- prop.table(cell_type_counts, margin = 2) * 100

# Convert to data frame for plotting
cell_type_proportions_df <- as.data.frame(cell_type_proportions)
colnames(cell_type_proportions_df) <- c("CellType", "Sample", "Proportion")

6.3.2. Plot Cell-Type Proportions

library(RColorBrewer)
# Combine "Set1" and "Set2" color palettes
set1_colors <- brewer.pal(9, "Set1")    # Set1 has up to 9 colors
set2_colors <- brewer.pal(8, "Set2")    # Set2 has up to 8 colors
combined_colors <- c(set1_colors, set2_colors)  # Combine both palettes

# Plot with combined colors
ggplot(cell_type_proportions_df, aes(x = Sample, y = Proportion, fill = CellType)) +
  geom_bar(stat = "identity", position = "fill") +
  ylab("Proportion (%)") +
  ggtitle("Cell-Type Proportions per Sample") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(values = combined_colors)

6.3.3. Explain Sample Variability in Cell-Type Proportions

# Summarize the number of cells per assigned cell type
cell_type_counts <- table(merged_seurat_batch_corrected$AssignedCellType)

# Convert the table to a data frame for better display
cell_type_df <- as.data.frame(cell_type_counts)
colnames(cell_type_df) <- c("Cell_Type", "Cell_Count")

# Calculate median UMI counts per cell type
umi_counts <- merged_seurat_batch_corrected@meta.data %>%
  group_by(AssignedCellType) %>%
  summarize(Median_UMI_Count = median(nCount_RNA))

# Calculate average number of genes expressed per cell type
genes_expressed <- merged_seurat_batch_corrected@meta.data %>%
  group_by(AssignedCellType) %>%
  summarize(Average_Genes_Expressed = mean(nFeature_RNA))

# Calculate average mitochondrial percentage per cell type
mt_percentage <- merged_seurat_batch_corrected@meta.data %>%
  group_by(AssignedCellType) %>%
  summarize(Average_Mitochondrial_Percentage = mean(percent.mt))

# Calculate number of clusters per cell type
clusters_per_cell_type <- merged_seurat_batch_corrected@meta.data %>%
  group_by(AssignedCellType) %>%
  summarize(Clusters = n_distinct(seurat_clusters))

# Merge all metrics into the cell_type_df
cell_type_df <- cell_type_df %>%
  left_join(umi_counts, by = c("Cell_Type" = "AssignedCellType")) %>%
  left_join(genes_expressed, by = c("Cell_Type" = "AssignedCellType")) %>%
  left_join(mt_percentage, by = c("Cell_Type" = "AssignedCellType")) %>%
  left_join(clusters_per_cell_type, by = c("Cell_Type" = "AssignedCellType"))

# Display the updated table
knitr::kable(
  cell_type_df,
  caption = "Identified Cell Types with Additional Statistics",
  format = "markdown",
  digits = 2
)

Observations:

Dominant Cell Populations:

• Common Lymphoid Progenitors (CLP) are the most abundant cell type, comprising 17.38% of the dataset with 3,249 cells. This high proportion indicates a significant lymphoid lineage presence, suggesting active lymphopoiesis or a focus on lymphoid cells in the sampled tissue.
• Granulocyte-Monocyte Progenitors (GMP_Neutrophils) represent the second-largest group at 12.94% with 2,419 cells. Their abundance points to robust granulocytic activity, which could be associated with an immune response or bone marrow function.

Transcriptional Activity Indicators:

• Erythrocytes exhibit the highest Median UMI Count of 4,358.0 and the highest Average Genes Expressed at 1,905.43. This suggests that erythroid cells in the dataset are highly transcriptionally active, likely due to hemoglobin production during erythropoiesis.
• Basophils also show high transcriptional activity with a Median UMI Count of 3,747.5 and an Average Genes Expressed of 1,896.42. Basophils play roles in inflammatory responses, and their activity might indicate an allergic reaction or parasitic infection.

Low Representation of Certain Cell Types:

• are the least abundant, making up only 1.31% of the dataset with 244 cells. This low number could reflect a reduced monocytic activity.
• Lymphoid-Primed Multipotent Progenitors (LMPP) and Natural Killer (NK) Cells also have lower percentages (1.98% and 2.54%, respectively), suggesting they are less prevalent in the sampled tissue.

Cluster Distribution and Cell Heterogeneity:

• The number of clusters associated with each cell type ranges from 9 to 19, indicating varying degrees of heterogeneity within cell populations.
• CLP cells are distributed across 19 clusters, the highest among all cell types, highlighting significant heterogeneity and possibly different maturation or activation states within this group.
• CD16_Monocytes are found in 9 clusters, suggesting a relatively more homogeneous population compared to other cell types.

Comparative Cell Activity:

• GMP_Neutrophils have a high Median UMI Count (3,255.0) and a high Average Genes Expressed (1,701.75), indicating active transcriptional profiles, possibly due to their role in innate immunity and rapid response to infection.

-Plasma Cells show significant transcriptional activity (Median UMI Count of 2,927.0), aligning with their role in producing large amounts of antibodies.

Biological Implications:

• The prominence of progenitor cells (CLP, GMP_Neutrophils, HSC) suggests that the sample may be derived from bone marrow or a hematopoietic organ, indicating active hematopoiesis.
• The presence of various immune cells (Basophils, pDCs, Plasma Cells, T Cells, B Cells) points to a diverse immune environment, which could be indicative of an immune response or a state of immunological readiness.

7. Differential Analysis

7.1. Differential Expression Analysis on cell-types

We will use the manual annotations stored under the AssignedCellTypes for the differential expression analysis.

7.1.1. Subset the Data for Relevant Cell Types

# Modify the @meta.data slot
merged_seurat_batch_corrected@meta.data <- merged_seurat_batch_corrected@meta.data %>%
  mutate(CellGroup = case_when(
    AssignedCellType %in% c("CD4_T", "CD8_T") ~ "T_cells",
    AssignedCellType == "B" ~ "B_cells",
    AssignedCellType %in% c("CD14_Monocytes", "CD16_Monocytes") ~ "Monocytes",
    TRUE ~ NA_character_
  ))

# Verify the new metadata
head(merged_seurat_batch_corrected@meta.data)

##                              orig.ident nCount_RNA nFeature_RNA    Sample Donor
## BMMC_D1T1:AAACCCAAGATGCAGC-1       BMMC       2433         1352 BMMC_D1T1    D1
## BMMC_D1T1:AAACCCACAAACTCGT-1       BMMC       5106         2001 BMMC_D1T1    D1
## BMMC_D1T1:AAACCCACAGTGTACT-1       BMMC       3589         1441 BMMC_D1T1    D1
## BMMC_D1T1:AAACCCATCGCTATTT-1       BMMC       3603         1809 BMMC_D1T1    D1
## BMMC_D1T1:AAACGAACACCCAATA-1       BMMC       2065         1106 BMMC_D1T1    D1
## BMMC_D1T1:AAACGAACAGCAGTCC-1       BMMC       2783         1380 BMMC_D1T1    D1
##                              Replicate Sex log10GenesPerUMI percent.mt
## BMMC_D1T1:AAACCCAAGATGCAGC-1        T1   F        0.9246442          0
## BMMC_D1T1:AAACCCACAAACTCGT-1        T1   F        0.8902846          0
## BMMC_D1T1:AAACCCACAGTGTACT-1        T1   F        0.8885197          0
## BMMC_D1T1:AAACCCATCGCTATTT-1        T1   F        0.9158690          0
## BMMC_D1T1:AAACGAACACCCAATA-1        T1   F        0.9181987          0
## BMMC_D1T1:AAACGAACAGCAGTCC-1        T1   F        0.9115596          0
##                              mitoRatio integrated_snn_res.0.5 seurat_clusters
## BMMC_D1T1:AAACCCAAGATGCAGC-1         0                      4               4
## BMMC_D1T1:AAACCCACAAACTCGT-1         0                      0               0
## BMMC_D1T1:AAACCCACAGTGTACT-1         0                      0               0
## BMMC_D1T1:AAACCCATCGCTATTT-1         0                     14              14
## BMMC_D1T1:AAACGAACACCCAATA-1         0                     17              17
## BMMC_D1T1:AAACGAACAGCAGTCC-1         0                      4               4
##                                SingleR.labels manual_annotations        HSC
## BMMC_D1T1:AAACCCAAGATGCAGC-1          T_cells               LMPP -0.2085936
## BMMC_D1T1:AAACCCACAAACTCGT-1 Pre-B_cell_CD34-               LMPP -0.2465950
## BMMC_D1T1:AAACCCACAGTGTACT-1         Monocyte               LMPP  0.4229277
## BMMC_D1T1:AAACCCATCGCTATTT-1           B_cell              CD8_T  0.5148590
## BMMC_D1T1:AAACGAACACCCAATA-1 Pro-B_cell_CD34+               LMPP  0.5262178
## BMMC_D1T1:AAACGAACAGCAGTCC-1          T_cells               LMPP -0.1922874
##                                    LMPP        CLP GMP_Neutrophils            B
## BMMC_D1T1:AAACCCAAGATGCAGC-1 -0.1167382  3.1682689      -0.5261399 -0.079430080
## BMMC_D1T1:AAACCCACAAACTCGT-1 -0.3291861 -0.8512616      -0.4199632 -0.081091526
## BMMC_D1T1:AAACCCACAGTGTACT-1  0.6323619 -0.7382036      -0.4199048 -0.061618157
## BMMC_D1T1:AAACCCATCGCTATTT-1  0.4255064  0.6745708      -0.6319220  1.362775836
## BMMC_D1T1:AAACGAACACCCAATA-1 -0.2103655 -0.3884754      -0.4742133 -0.272816692
## BMMC_D1T1:AAACGAACAGCAGTCC-1  0.1173536  0.9333633      -0.4261055 -0.004743017
##                                   Pre_B     Plasma      CD8_T      CD4_T
## BMMC_D1T1:AAACCCAAGATGCAGC-1 -0.1117151 -0.2268927  0.5233883  1.1195594
## BMMC_D1T1:AAACCCACAAACTCGT-1 -0.1500222 -0.2341316 -0.5226325 -0.6692314
## BMMC_D1T1:AAACCCACAGTGTACT-1 -0.1536805 -0.2544634 -0.5079602 -0.6507918
## BMMC_D1T1:AAACCCATCGCTATTT-1  0.5369676  0.4410234  0.4934183  0.6208000
## BMMC_D1T1:AAACGAACACCCAATA-1  0.6529881  1.0962325 -0.2753458 -0.3820186
## BMMC_D1T1:AAACGAACAGCAGTCC-1 -0.1095558 -0.1775215  0.9705960  1.8400567
##                                       NK Erythrocytes         pDC        cDC
## BMMC_D1T1:AAACCCAAGATGCAGC-1 -0.10592471   -0.2510767 -0.19364466 -0.1346062
## BMMC_D1T1:AAACCCACAAACTCGT-1 -0.19700281   -0.4282465  0.02641115  0.3435554
## BMMC_D1T1:AAACCCACAGTGTACT-1  0.05971542   -0.4118805 -0.40121880  0.4726568
## BMMC_D1T1:AAACCCATCGCTATTT-1 -0.23476821   -0.3168708  0.55381030  0.1685869
## BMMC_D1T1:AAACGAACACCCAATA-1 -0.16380162   -0.1636095  0.37842650 -0.1680493
## BMMC_D1T1:AAACGAACAGCAGTCC-1 -0.03121210    0.1350610 -0.18757173 -0.1118684
##                              CD14_Monocytes CD16_Monocytes   Basophils
## BMMC_D1T1:AAACCCAAGATGCAGC-1    -0.10938439     -0.3115009 -0.07299451
## BMMC_D1T1:AAACCCACAAACTCGT-1     0.43058666      1.4023286 -0.15366612
## BMMC_D1T1:AAACCCACAGTGTACT-1     0.31791125      1.6491804 -0.27602482
## BMMC_D1T1:AAACCCATCGCTATTT-1    -0.20229324     -0.2885670 -0.20698645
## BMMC_D1T1:AAACGAACACCCAATA-1    -0.08361119     -0.2263020 -0.16210423
## BMMC_D1T1:AAACGAACAGCAGTCC-1    -0.10049795     -0.2797639 -0.11043150
##                              AssignedCellType CellGroup
## BMMC_D1T1:AAACCCAAGATGCAGC-1              CLP      <NA>
## BMMC_D1T1:AAACCCACAAACTCGT-1   CD16_Monocytes Monocytes
## BMMC_D1T1:AAACCCACAGTGTACT-1   CD16_Monocytes Monocytes
## BMMC_D1T1:AAACCCATCGCTATTT-1                B   B_cells
## BMMC_D1T1:AAACGAACACCCAATA-1           Plasma      <NA>
## BMMC_D1T1:AAACGAACAGCAGTCC-1            CD4_T   T_cells

table(merged_seurat_batch_corrected@meta.data$CellGroup, useNA = "ifany")

## 
##   B_cells Monocytes   T_cells      <NA> 
##       842      2617      1998     13236

# Subset the data
subset_data <- subset(
  merged_seurat_batch_corrected,
  subset = CellGroup %in% c("B_cells", "T_cells", "Monocytes")
)

# Set 'CellGroup' as the active identity
Idents(subset_data) <- "CellGroup"

7.1.3. Perform Differential Expression Analysis

7.1.3.1. B Cells vs T Cells

# Find markers between B cells and T cells
markers_B_vs_T <- FindMarkers(
  object = subset_data,
  ident.1 = "B_cells",
  ident.2 = "T_cells",
  min.pct = 0.25,
  logfc.threshold = 0.25
)

7.1.3.2. T Cells vs Monocytes

# Find markers between T cells and Monocytes
markers_T_vs_Mono <- FindMarkers(
  object = subset_data,
  ident.1 = "T_cells",
  ident.2 = "Monocytes",
  min.pct = 0.25,
  logfc.threshold = 0.25
)

7.1.3.4 Prepare Data for Volcano Plots

# Add gene names to the results
markers_B_vs_T$gene <- rownames(markers_B_vs_T)
markers_T_vs_Mono$gene <- rownames(markers_T_vs_Mono)

# Ensure adjusted p-values are available
if (!"p_val_adj" %in% colnames(markers_B_vs_T)) {
  markers_B_vs_T$p_val_adj <- p.adjust(markers_B_vs_T$p_val, method = "BH")
}
if (!"p_val_adj" %in% colnames(markers_T_vs_Mono)) {
  markers_T_vs_Mono$p_val_adj <- p.adjust(markers_T_vs_Mono$p_val, method = "BH")
}

7.1.3.5. Create Volcano Plots

Volcano Plot for B Cells vs T Cells

logFC_cutoff <- 0.25
pval_cutoff <- 0.05
markers_B_vs_T$minus_log10_pval <- -log10(markers_B_vs_T$p_val_adj)

# Create a column to indicate significant genes
markers_B_vs_T$Significance <- "Not Significant"
markers_B_vs_T$Significance[
  (markers_B_vs_T$p_val_adj < pval_cutoff) & 
  (abs(markers_B_vs_T$avg_log2FC) > logFC_cutoff)] <- "Significant"

# Volcano plot
ggplot(markers_B_vs_T, aes(x = avg_log2FC, y = minus_log10_pval)) +
  geom_point(aes(color = Significance), alpha = 0.8, size = 1.5) +
  scale_color_manual(values = c("#F25C54", "#B3DEC1")) +
  xlab("Average Log2 Fold Change") +
  ylab("-Log10 Adjusted P-value") +
  ggtitle("Volcano Plot: B Cells vs T Cells") +
  theme_minimal()

Volcano Plot for T Cells vs Monocytes

markers_T_vs_Mono$minus_log10_pval <- -log10(markers_T_vs_Mono$p_val_adj)

# Create a column to indicate significant genes
markers_T_vs_Mono$Significance <- "Not Significant"
markers_T_vs_Mono$Significance[
  (markers_T_vs_Mono$p_val_adj < pval_cutoff) & 
  (abs(markers_T_vs_Mono$avg_log2FC) > logFC_cutoff)] <- "Significant"

# Volcano plot
ggplot(markers_T_vs_Mono, aes(x = avg_log2FC, y = minus_log10_pval)) +
  geom_point(aes(color = Significance), alpha = 0.8, size = 1.5) +
  scale_color_manual(values = c("#F25C54", "#B3DEC1")) +
  xlab("Average Log2 Fold Change") +
  ylab("-Log10 Adjusted P-value") +
  ggtitle("Volcano Plot: T Cells vs Monocytes") +
  theme_minimal()

7.1.1. Memory Formation on Cells

Naive T cells are mature T cells that have not yet encountered their specific antigen, whereas memory T cells are antigen-experienced cells that persist long-term and respond more rapidly upon re-exposure to the same antigen. Memory T cells arise from naive T cells following antigen exposure, undergoing differentiation processes that enable them to provide quicker and more specific immune responses during subsequent encounters with the same pathogen.

7.2. Plot Differentially Expressed Genes

markers_B_vs_T[] <- lapply(markers_B_vs_T, function(col) {
  if (is(col, "Rle")) as.vector(col) else col
})

markers_B_vs_T <- markers_B_vs_T %>%
  mutate(
    p_val_adj = as.numeric(p_val_adj),
    p_val_adj = ifelse(p_val_adj == 0 | is.na(p_val_adj), 1e-300, p_val_adj)
  )

markers_T_vs_Mono[] <- lapply(markers_T_vs_Mono, function(col) {
  if (is(col, "Rle")) as.vector(col) else col
})

markers_T_vs_Mono <- markers_T_vs_Mono %>%
  mutate(
    p_val_adj = as.numeric(p_val_adj),
    p_val_adj = ifelse(p_val_adj == 0 | is.na(p_val_adj), 1e-300, p_val_adj)
  )

# Top 5 genes for B Cells vs T Cells
top5_B_vs_T <- markers_B_vs_T %>%
  arrange(p_val_adj) %>%
  slice(1:5) %>%
  mutate(
    Comparison = "B_vs_T",
    CellType = ifelse(avg_log2FC > 0, "B_cells", "T_cells"),
    Gene = gene
  ) %>%
  select(Gene, CellType, avg_log2FC, p_val_adj, Comparison)

# Top 5 genes for T Cells vs Monocytes
top5_T_vs_Mono <- markers_T_vs_Mono %>%
  arrange(p_val_adj) %>%
  slice(1:5) %>%
  mutate(
    Comparison = "T_vs_Monocytes",
    CellType = ifelse(avg_log2FC > 0, "T_cells", "Monocytes"),
    Gene = gene
  ) %>%
  select(Gene, CellType, avg_log2FC, p_val_adj, Comparison)

# Combine the data
top_genes <- bind_rows(top5_B_vs_T, top5_T_vs_Mono)

# Calculate -log10 of adjusted p-values
top_genes <- top_genes %>%
  mutate(minus_log10_pval = -log10(p_val_adj))

# Print final output
print(top_genes)

# Create the plot
ggplot(top_genes, aes(x = CellType, y = Gene)) +
  geom_point(aes(size = minus_log10_pval, color = avg_log2FC)) +
  scale_color_gradient2(
    low = "blue",
    mid = "white",
    high = "red",
    midpoint = 0,
    limits = c(min(top_genes$avg_log2FC), max(top_genes$avg_log2FC))
  ) +
  facet_wrap(~ Comparison, scales = "free_x") +
  theme_minimal() +
  labs(
    x = "Cell Type",
    y = "Gene",
    color = "Log2 Fold Change",
    size = "-Log10 Adjusted P-value",
    title = "Top 5 Differentially Expressed Genes for Each Comparison"
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

8. Pathway Analysis

8.1. Differential Expression Analysis on groups

8.1.1. Differential Expression Analysis for BMMC vs CD34 (All Cell Types)

# Differential Expression Analysis for BMMC vs CD34 (All Cell Types)
bmmc_vs_cd34 <- FindMarkers(merged_seurat_batch_corrected, ident.1 = "BMMC", ident.2 = "CD34", group.by = "orig.ident")

8.1.2. Extract Top 5 DEGs

bmmc_vs_cd34$gene <- rownames(bmmc_vs_cd34)

top5_bmmc_vs_cd34 <- bmmc_vs_cd34 %>%
  arrange(p_val_adj) %>%
  head(5) %>%
  select(gene, avg_log2FC, p_val_adj)

knitr::kable(top5_bmmc_vs_cd34, caption = "Top 5 Differentially Expressed Genes (BMMC vs CD34)")

8.1.3. Differential Expression Analysis for Monocytes in BMMC vs CD34

# Subset for Monocyte cells only
monocytes <- subset(merged_seurat_batch_corrected, subset = SingleR.labels == "Monocyte")

# Differential Expression Analysis for Monocytes in BMMC vs CD34
monocyte_bmmc_vs_cd34 <- FindMarkers(monocytes, ident.1 = "BMMC", ident.2 = "CD34", group.by = "orig.ident")

monocyte_bmmc_vs_cd34$gene <- rownames(monocyte_bmmc_vs_cd34)
                                       
# Extract Top 5 DEGs by p-value for Monocytes
top5_monocyte_bmmc_vs_cd34 <- monocyte_bmmc_vs_cd34 %>%
  arrange(p_val_adj) %>%
  head(5) %>%
  select(gene, avg_log2FC, p_val_adj)

knitr::kable(top5_monocyte_bmmc_vs_cd34, caption = "Top 5 Differentially Expressed Genes in Monocytes (BMMC vs CD34)")

8.2. Pathway analysis on groups

Idents(merged_seurat_batch_corrected) <- "orig.ident"

top_pathways <- DEenrichRPlot(
  object = merged_seurat_batch_corrected,
  ident.1 = "BMMC",
  ident.2 = "CD34",
  balanced = TRUE,            # Show both upregulated and downregulated terms
  logfc.threshold = 0.25,     # Only consider genes with a minimum log fold change
  test.use = "wilcox",        # Use Wilcoxon Rank Sum test (default)
  p.val.cutoff = 0.05,        # Only include significant DE genes (adjusted p < 0.05)
  enrich.database = "GO_Biological_Process_2021",  # Choose GO Biological Process terms
  num.pathway = 10,            # Display top 10 enriched pathways
  max.genes = 500
)

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
## Parsing results... Done.

8.3. Biological interpretation

# Filter the results to include only significant genes
significant_genes <- bmmc_vs_cd34 %>%
  filter(p_val_adj < 0.05) %>%
  arrange(p_val_adj) %>%
  pull(gene)

# Use EnrichR to perform enrichment analysis on significant genes
enrichr_results <- enrichr(
  significant_genes,
  databases = "GO_Biological_Process_2021"
)

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
## Parsing results... Done.

# Extract the GO Biological Process results and find the top pathway
go_results <- enrichr_results[["GO_Biological_Process_2021"]]

# Arrange the results by adjusted p-value and select the top pathway
top_pathway <- go_results %>%
  arrange(Adjusted.P.value) %>%
  head(1)

kable(top_pathway, 
      format = "html", 
      col.names = c("Term", "Overlap", "P.value", "Adjusted P-value", "Old P-value", "Old Adjusted P-value", "Odds Ratio", "Combined Score", "Genes"), 
      caption = "Pathway with the lowest P-value")

Neutrophil degranulation (GO:0043312) is the controlled release of secretory granules from neutrophils, which contain pre-stored substances like proteases, lipases, and inflammatory mediators. QuickGo

Neutrophil degranulation, BMMCs, and CD34+ cells are interconnected components of the immune system that work together in response to pathogens or inflammatory stimuli, influencing each other’s activity through cytokine signaling and cellular interactions.

The interactions is crucial in conditions like allergies or asthma, where both neutrophil activity and mast cell degranulation contribute to the pathophysiology. Dysregulation of these processes can lead to excessive inflammation or tissue damage. Proteome studies has revealed that some proteins (MPO, MMP9, DEFA1) associated with ‘neutrophil degranulation’ showed the presence of ‘signal sequence’ suggesting their potential as circulatory markers for early detection of Gall Bladder Cancer. [(Gautam et al, 2023)][https://pmc.ncbi.nlm.nih.gov/articles/PMC9853450/]. SAARS-CoV-2 infection has been associated with the result of dzsregulated neutrophil degranulation as well. (Raid et al (2022)

Part 4 — Trajectory & Cell-Cell Communication

9. Trajectory Analysis

9.1. Select Subset

9.1.1. Selecting a Group of Cells for Trajectory Analysis

Selected Group: Common Lymphoid Progenitors (CLP)

Rationale for Selection:

• CLPs are progenitor cells that differentiate into various lymphoid lineages, including B cells, NK cells and T cells. Studying CLPs allows us to explore the early stages of lymphoid development.
• Heterogeneity: The dataset contains a substantial number of CLP cells spread across multiple clusters, indicating heterogeneity and the presence of different developmental stages.
• Biological Significance: Understanding the trajectory of CLP cells can provide insights into immune system development and potential dysregulation in diseases like leukemia.

Alternative Group: Hematopoietic Stem Cells (HSC)

• Reason: HSCs are the foundational multipotent stem cells in hematopoiesis, giving rise to all blood cell types. Analyzing HSCs can reveal early differentiation events and lineage commitment.

9.1.2. Performing Trajectory Analysis with Monocle 3

9.1.2.1. Prepare Data

# Subset the Seurat object
CLP_cells <- subset(merged_seurat_batch_corrected, subset = AssignedCellType == "CLP")

expression_matrix <- GetAssayData(CLP_cells, slot = "counts", assay = "RNA")
cell_metadata <- CLP_cells@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(expression_matrix))
rownames(gene_annotation) <- rownames(expression_matrix)

# Create Monocle 3 CellDataSet
CLP_cds <- new_cell_data_set(
  expression_data = expression_matrix,
  cell_metadata = cell_metadata,
  gene_metadata = gene_annotation
)

9.1.2.2. Preprocess

CLP_cds <- preprocess_cds(CLP_cds, num_dim = 50)

9.1.3. UMAP, Clustering and Visualisation

# Reduce dimensions using UMAP
CLP_cds <- reduce_dimension(CLP_cds, reduction_method = "UMAP")

# Cluster the cells using Monocle 3's clustering algorithm
CLP_cds <- cluster_cells(CLP_cds, resolution = 1e-3)

# Plot clusters
p1 <- plot_cells(CLP_cds, color_cells_by = "cluster", show_trajectory_graph = FALSE)

# Plot partitions
p2 <- plot_cells(CLP_cds, color_cells_by = "partition", show_trajectory_graph = FALSE)

# Combine plots side by side
library(patchwork)
p1 + p2

### 9.1.4. Learn the Trajectory Graph

# Learn the trajectory graph
CLP_cds <- learn_graph(CLP_cds, use_partition = TRUE, verbose = FALSE)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

9.1.5. Plot Trajectory Colored by Cluster

set.seed(12346465)
# Plot trajectory colored by clusters
plot_cells(CLP_cds,
           color_cells_by = "cluster",
           label_cell_groups = FALSE,
           label_leaves = FALSE,
           label_branch_points = FALSE) +
  ggtitle("Trajectory Clustering of CLP Cells") +
  theme_minimal()

9.2. Select Root Nodes Manually

table(clusters(CLP_cds))

## 
##    1    2    3    4    5    6    7    8    9 
## 1120 1002  490  206  108   66   60   54   26

set.seed(12308)
CLP_markers <- "IL7R"
plot_cells(CLP_cds,
           genes = "IL7R",
           show_trajectory_graph = TRUE,
           label_cell_groups = TRUE,
           label_leaves = TRUE,
           label_branch_points = FALSE) +
  ggtitle("Expression of IL7R in UMAP") +
  theme_minimal()

9.2.1 Selecting Root Nodes and Performing Trajectory Analysis

9.2.1.1. Selecting Root Nodes

Cluster 1

# Select specific cells in the chosen starting cluster
root_cells <- colnames(CLP_cds)[clusters(CLP_cds) == "1"]

# Order cells based on specific root cells
CLP_cds <- order_cells(CLP_cds, root_cells = root_cells)

# Order cells along the trajectory with the specified root cells
CLP_cds <- order_cells(CLP_cds, root_cells = root_cells)

p1 <- plot_cells(CLP_cds,
                 color_cells_by = "pseudotime",
                 label_cell_groups = FALSE,
                 label_leaves = TRUE,
                 label_branch_points = TRUE) +
  ggtitle("Trajectory of CLP Cells with Pseudotime") +
  theme_minimal()

p2 <- plot_cells(CLP_cds,
                 color_cells_by = "partition",
                 label_cell_groups = TRUE,
                 label_leaves = TRUE,
                 label_branch_points = TRUE) +
  ggtitle("Trajectory of CLP Cells by Cluster") +
  theme_minimal()

# Plot them side by side
gridExtra::grid.arrange(p1, p2, ncol = 2)

Observation

Cells are colored according to pseudotime values, progressing from dark purple (low pseudotime) to yellow (high pseudotime). This gradient visually represents the order of cells along a hypothetical developmental or differentiation path.

Cells are colored based on assigned partitions or clusters, each represented by a unique color. For example:

Cluster 1 is in red, which is the most prominent. Other clusters, such as 2 (yellow), 3 (green), 4 (blue), and 5 (purple), represent smaller or more isolated cell populations.

Cluster 2

# Select specific cells in the chosen starting cluster
root_cells <- colnames(CLP_cds)[clusters(CLP_cds) == "2"]

# Order cells based on specific root cells
CLP_cds <- order_cells(CLP_cds, root_cells = root_cells)

# Order cells along the trajectory with the specified root cells
CLP_cds <- order_cells(CLP_cds, root_cells = root_cells)

p1 <- plot_cells(CLP_cds,
                 color_cells_by = "pseudotime",
                 label_cell_groups = FALSE,
                 label_leaves = TRUE,
                 label_branch_points = TRUE) +
  ggtitle("Trajectory of CLP Cells with Pseudotime") +
  theme_minimal()

p2 <- plot_cells(CLP_cds,
                 color_cells_by = "partition",
                 label_cell_groups = TRUE,
                 label_leaves = TRUE,
                 label_branch_points = TRUE) +
  ggtitle("Trajectory of CLP Cells by Cluster") +
  theme_minimal()

# Plot them side by side
gridExtra::grid.arrange(p1, p2, ncol = 2)

Observation:

Cell Differentiation: This trajectory analysis suggests a potential differentiation pathway for CLP cells. Cells begin in one state (low pseudotime) and transition through several stages (higher pseudotime), representing different phases of development or activation.

Why is the selection of the root nodes important for the algorithm?

The root node determines the starting point of pseudotime, defining the directionality of the trajectory and how cells are ordered along it. An appropriate root ensures the trajectory reflects true biological progression, such as differentiation or development, while an incorrect root may lead to reversed or biologically inconsistent trajectories.

Which points are a good choice for root nodes of the analysis and why?

Good root nodes represent early progenitor or stem-like cells, characterized by high expression of early markers (e.g., IL7R for CLP cells). These cells are biologically at the start of differentiation pathways and typically cluster near the trajectory’s origin, ensuring pseudotime aligns with known developmental stages.

9.3 Select Root Nodes Automatically

CLP_cds <- cluster_cells(CLP_cds)
# Specify the marker gene to identify the root cluster
marker_gene <- "IL7R" 

# Extract expression data for the marker gene across all cells
expression_data <- exprs(CLP_cds)[marker_gene, ]

# Identify the cluster with the highest average expression of the marker gene
root_cluster <- names(which.max(tapply(expression_data, CLP_cds@clusters@listData[["UMAP"]][["clusters"]], mean)))

# Select the cells that belong to the identified root cluster
root_cells <- colnames(CLP_cds[, clusters(CLP_cds) == root_cluster])
# Order cells for trajectory analysis, specifying the root cells
CLP_cds_auto <- order_cells(CLP_cds, root_cells = root_cells)

# Highlight the path to B cells
# Identify cells expressing B cell markers
b_cell_markers <- c("CD19", "MS4A1")
b_cell_expression <- rowData(CLP_cds)$gene_short_name %in% b_cell_markers

# Plot the trajectory emphasizing the B cell path
plot_cells(CLP_cds,
           genes = b_cell_markers,
           show_trajectory_graph = TRUE,
           label_cell_groups = FALSE,
           label_leaves = TRUE,
           label_branch_points = FALSE) +
  ggtitle("Trajectory Path from CLP to B Cells") +
  theme_minimal()

Observation

The plot tells the story of cell differentiation, showing how progenitor cells (low CD19/MS4A1) move through intermediate states and eventually become mature B cells (high CD19/MS4A1), evidenced by by two key B cell markers: CD19 (source: table 2) and MS4A1proteinatlas.

Tracing the Biological Path

Beginning of Trajectory:

These clusters are located near the start of the black trajectory line and have low expression (purple shading) of CD19 and MS4A1, representing early progenitor or undifferentiated CLP cells.

Middle of the Path:

As we move along the black line, clusters start showing moderate expression (shading transition toward green) of CD19, indicating these cells are beginning to commit to the B cell lineage.**

End of the Path: These clusters show strong green shading for both CD19 and MS4A1, marking the final stages of B cell differentiation. Cells here are mature B cells, having fully acquired the B cell marker expression.

Manual vs Automatic root node selection

We think selecting cluster 1 (high IL7R expression) as the root node provided a meaningful trajectory from CLP to mature B cells. Both manual and automatic methods led to consistent pseudotime ordering.

10. Cell-Cell Communication

10.0. Subsetting Common Cell Types

bmmc_seurat <- subset(merged_seurat_batch_corrected, subset = orig.ident == "BMMC")
cd34_seurat <- subset(merged_seurat_batch_corrected, subset = orig.ident == "CD34")

cat("BMMC cells:", ncol(bmmc_seurat), "| CD34 cells:", ncol(cd34_seurat), "\n")

## BMMC cells: 11343 | CD34 cells: 7350

# Subset the data by origin to get cell types for each group
bmmc_cell_types <- unique(merged_seurat_batch_corrected@meta.data[merged_seurat_batch_corrected@meta.data$orig.ident == "BMMC", "AssignedCellType"])
cd34_cell_types <- unique(merged_seurat_batch_corrected@meta.data[merged_seurat_batch_corrected@meta.data$orig.ident == "CD34", "AssignedCellType"])

# Find the common cell types
common_cell_types <- intersect(bmmc_cell_types, cd34_cell_types)
common_cell_types

##  [1] "CLP"             "CD16_Monocytes"  "B"               "Plasma"         
##  [5] "CD4_T"           "CD8_T"           "GMP_Neutrophils" "cDC"            
##  [9] "pDC"             "CD14_Monocytes"  "Erythrocytes"    "NK"             
## [13] "LMPP"            "Basophils"       "HSC"             "Pre_B"

# Task 10 restricts the analysis to cell types present in BOTH groups.
# Verify that explicitly rather than assuming it.
cat("Cell types in BMMC:", length(bmmc_cell_types),
    "| in CD34:", length(cd34_cell_types),
    "| common:", length(common_cell_types), "\n")

## Cell types in BMMC: 16 | in CD34: 16 | common: 16

cat("BMMC-only:", paste(setdiff(bmmc_cell_types, cd34_cell_types), collapse = ", "), "\n")

## BMMC-only:

cat("CD34-only:", paste(setdiff(cd34_cell_types, bmmc_cell_types), collapse = ", "), "\n")

## CD34-only:

10.1. Create CellChat objects for each group using the common cell types

# Create CellChat object for BMMC samples
cellchat_bmmc <- createCellChat(object = bmmc_seurat, group.by = "AssignedCellType", assay = "RNA")

## [1] "Create a CellChat object from a Seurat object"
## The `meta.data` slot in the Seurat object is used as cell meta information 
## Set cell identities for the new CellChat object 
## The cell groups used for CellChat analysis are  B Basophils CD14_Monocytes CD16_Monocytes CD4_T CD8_T cDC CLP Erythrocytes GMP_Neutrophils HSC LMPP NK pDC Plasma Pre_B

# Create CellChat object for CD34 samples
cellchat_cd34 <- createCellChat(object = cd34_seurat, group.by = "AssignedCellType", assay = "RNA")

## [1] "Create a CellChat object from a Seurat object"
## The `meta.data` slot in the Seurat object is used as cell meta information 
## Set cell identities for the new CellChat object 
## The cell groups used for CellChat analysis are  B Basophils CD14_Monocytes CD16_Monocytes CD4_T CD8_T cDC CLP Erythrocytes GMP_Neutrophils HSC LMPP NK pDC Plasma Pre_B

10.3. Set the CellChat database, preprocesssing with visualization.

We chose the following pathway for visualisation.

MIF (Macrophage migration Inhibitory Factor)

Why this pathway: the task asks for a circle plot of one pathway in each group, so the pathway has to be significant in BMMC and CD34. MIF is the strongest such pathway in this dataset — it carries the highest summed communication probability of the ten shared pathways (BMMC 20.8, CD34 6.1; see the comparison table below).

Why it’s interesting: MIF is a pleiotropic pro-inflammatory cytokine that signals through the CD74 receptor together with CD44 and the chemokine co-receptors CXCR2/CXCR4. Those receptors are broadly expressed across hematopoietic lineages, which is consistent with MIF being one of the few pathways active in both the mature (BMMC) and progenitor-enriched (CD34+) compartments.

Note: an earlier version of this analysis visualised WNT. WNT is prominent in the CellChatDB reference catalogue but is not inferred as significant in either group here, so those plots produced no output. Selecting from the inferred pathways avoids that.

# Load the CellChat database for human
CellChatDB <- CellChatDB.human
cellchat_bmmc@DB <- CellChatDB
cellchat_cd34@DB <- CellChatDB


# Preprocess the expression data and compute communication probabilities.
# NOTE: `for (cellchat in list(...))` binds a COPY each iteration, so every
# `cellchat <- ...` below would be discarded at the end of the iteration and the
# named objects would keep only their @DB slot. We iterate over a named list and
# write the processed object back so the inferred results actually persist.
cellchat_list <- list(BMMC = cellchat_bmmc, CD34 = cellchat_cd34)

for (grp in names(cellchat_list)) {
  cellchat <- cellchat_list[[grp]]
  
  tryCatch({
    print("Subsetting data...")
    cellchat <- subsetData(cellchat)
  }, error = function(e) {
    cat("Error in subsetting data:", e$message, "\n")
  })
  
  tryCatch({
    print("Identifying over-expressed genes...")
    cellchat <- identifyOverExpressedGenes(cellchat)
  }, error = function(e) {
    cat("Error in identifying over-expressed genes:", e$message, "\n")
  })
  
  tryCatch({
    print("Identifying over-expressed interactions...")
    cellchat <- identifyOverExpressedInteractions(cellchat)
  }, error = function(e) {
    cat("Error in identifying over-expressed interactions:", e$message, "\n")
  })
  
  tryCatch({
    print("Projecting data onto PPI network...")
    cellchat <- projectData(cellchat, PPI.human)
  }, error = function(e) {
    cat("Error in projecting data onto PPI network:", e$message, "\n")
  })
  
  tryCatch({
    print("Computing communication probabilities...")
    cellchat <- computeCommunProb(cellchat)
  }, error = function(e) {
    cat("Error in computing communication probabilities:", e$message, "\n")
  })
  
  tryCatch({
    print("Filtering communications with minimum cell threshold...")
    cellchat <- filterCommunication(cellchat, min.cells = 10)
  }, error = function(e) {
    cat("Error in filtering communications:", e$message, "\n")
  })
  
  tryCatch({
    print("Computing communication probabilities for pathways...")
    cellchat <- computeCommunProbPathway(cellchat)
  }, error = function(e) {
    cat("Error in computing communication probabilities for pathways:", e$message, "\n")
  })
  
  tryCatch({
    print("Aggregating network data...")
    cellchat <- aggregateNet(cellchat)
  }, error = function(e) {
    cat("Error in aggregating network data:", e$message, "\n")
  })
  
  tryCatch({
    groupSize <- as.numeric(table(cellchat@idents))
    print("Visualizing interaction counts in a circle plot...")
    par(mfrow = c(1, 2), mar = c(5, 4, 4, 2) + 0.1, xpd = TRUE)
    netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge = F, title.name = "Number of interactions")
    
    print("Visualizing interaction weights in a circle plot...")
    netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge = F, title.name = "Interaction weights/strength")
    
    mat <- cellchat@net$weight
  }, error = function(e) {
    cat("Error in visualizing interaction counts or weights:", e$message, "\n")
  })
  
  tryCatch({
    print("Visualizing individual rows of the interaction matrix in circle plots...")
    par(mfrow = c(2, 2), mar = c(5, 4, 4, 2) + 0.1, xpd = TRUE)
    for (i in 1:nrow(mat)) {
      print(paste("Visualizing interaction for:", rownames(mat)[i]))
      mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
      mat2[i, ] <- mat[i, ]
      netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
    }
  }, error = function(e) {
    cat("Error in visualizing individual rows of the interaction matrix:", e$message, "\n")
  })
  
  # Visualise a pathway that is actually significant for THIS group.
  # MIF is the strongest pathway shared by BMMC and CD34; fall back to the
  # group's own top pathway if it is somehow absent, so the plots always render.
  pathways.show <- if ("MIF" %in% cellchat@netP$pathways) "MIF" else cellchat@netP$pathways[1]
  cat("Group", grp, "- visualising pathway:", pathways.show, "\n")
  
  tryCatch({
    print("Creating hierarchy plot for specified pathways...")
    svg("hierarchy_plot.svg", width = 10, height = 8)
    vertex.receiver = seq(1, 4) # a numeric vector
    netVisual_aggregate(cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
    dev.off()
  }, error = function(e) {
    cat("Error in creating hierarchy plot:", e$message, "\n")
  })
  
  tryCatch({
    print("Creating circle plot for specified pathways...")
    svg("circle_plot.svg", width = 10, height = 8)
    netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")
    dev.off()
  }, error = function(e) {
    cat("Error in creating circle plot:", e$message, "\n")
  })
  
  tryCatch({
    print("Creating chord diagram for specified pathways...")
    svg("chord_diagram.svg", width = 10, height = 8)
    netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord")
    dev.off()
  }, error = function(e) {
    cat("Error in creating chord diagram:", e$message, "\n")
  })
  
  tryCatch({
    print("Creating heatmap for specified pathways...")
    svg("heatmap.svg", width = 10, height = 8)
    netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
    dev.off()
  }, error = function(e) {
    cat("Error in creating heatmap:", e$message, "\n")
  })

  # persist the fully-processed object for this group
  cellchat_list[[grp]] <- cellchat
}

## [1] "Subsetting data..."
## [1] "Identifying over-expressed genes..."
## [1] "Identifying over-expressed interactions..."
## [1] "Projecting data onto PPI network..."
## [1] "Computing communication probabilities..."
## triMean is used for calculating the average gene expression per cell group. 
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-08-19 22:03:42]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-08-19 22:07:49]"
## [1] "Filtering communications with minimum cell threshold..."
## [1] "Computing communication probabilities for pathways..."
## [1] "Aggregating network data..."
## [1] "Visualizing interaction counts in a circle plot..."

## [1] "Visualizing interaction weights in a circle plot..."

## [1] "Visualizing individual rows of the interaction matrix in circle plots..."
## [1] "Visualizing interaction for: B"

## [1] "Visualizing interaction for: Basophils"

## [1] "Visualizing interaction for: CD14_Monocytes"

## [1] "Visualizing interaction for: CD16_Monocytes"

## [1] "Visualizing interaction for: CD4_T"

## [1] "Visualizing interaction for: CD8_T"

## [1] "Visualizing interaction for: cDC"

## [1] "Visualizing interaction for: CLP"

## [1] "Visualizing interaction for: Erythrocytes"

## [1] "Visualizing interaction for: GMP_Neutrophils"

## [1] "Visualizing interaction for: HSC"

## [1] "Visualizing interaction for: LMPP"

## [1] "Visualizing interaction for: NK"

## [1] "Visualizing interaction for: pDC"

## [1] "Visualizing interaction for: Plasma"

## [1] "Visualizing interaction for: Pre_B"

## Group BMMC - visualising pathway: MIF 
## [1] "Creating hierarchy plot for specified pathways..."

## [1] "Creating circle plot for specified pathways..."

## [1] "Creating chord diagram for specified pathways..."

## [1] "Creating heatmap for specified pathways..."
## [1] "Subsetting data..."
## [1] "Identifying over-expressed genes..."
## [1] "Identifying over-expressed interactions..."
## [1] "Projecting data onto PPI network..."
## [1] "Computing communication probabilities..."
## triMean is used for calculating the average gene expression per cell group. 
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-08-19 22:08:21]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-08-19 22:11:58]"
## [1] "Filtering communications with minimum cell threshold..."
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  CD8_T 
## [1] "Computing communication probabilities for pathways..."
## [1] "Aggregating network data..."
## [1] "Visualizing interaction counts in a circle plot..."

## [1] "Visualizing interaction weights in a circle plot..."

## [1] "Visualizing individual rows of the interaction matrix in circle plots..."
## [1] "Visualizing interaction for: B"

## [1] "Visualizing interaction for: Basophils"

## [1] "Visualizing interaction for: CD14_Monocytes"

## [1] "Visualizing interaction for: CD16_Monocytes"

## [1] "Visualizing interaction for: CD4_T"

## [1] "Visualizing interaction for: CD8_T"

## [1] "Visualizing interaction for: cDC"

## [1] "Visualizing interaction for: CLP"

## [1] "Visualizing interaction for: Erythrocytes"

## [1] "Visualizing interaction for: GMP_Neutrophils"

## [1] "Visualizing interaction for: HSC"

## [1] "Visualizing interaction for: LMPP"

## [1] "Visualizing interaction for: NK"

## [1] "Visualizing interaction for: pDC"

## [1] "Visualizing interaction for: Plasma"

## [1] "Visualizing interaction for: Pre_B"

## Group CD34 - visualising pathway: MIF 
## [1] "Creating hierarchy plot for specified pathways..."

## [1] "Creating circle plot for specified pathways..."

## [1] "Creating chord diagram for specified pathways..."

## [1] "Creating heatmap for specified pathways..."

# recover the processed objects under their original names
cellchat_bmmc <- cellchat_list$BMMC
cellchat_cd34 <- cellchat_list$CD34

print("Process completed for all cellchat objects.")

## [1] "Process completed for all cellchat objects."

# sanity check: the inference must have persisted
stopifnot(!is.null(cellchat_bmmc@net$count), !is.null(cellchat_cd34@net$count))
cat("BMMC inferred pathways:", length(cellchat_bmmc@netP$pathways),
    "| CD34 inferred pathways:", length(cellchat_cd34@netP$pathways), "\n")

## BMMC inferred pathways: 23 | CD34 inferred pathways: 12

Common pathway

# Signalling pathways INFERRED FROM THE DATA.
# Note: cellchat@DB holds the reference CellChatDB catalogue, which is identical
# for both objects and would yield the full database rather than our results.
# The significant, data-derived pathways live in @netP$pathways.
pathways_bmmc <- cellchat_bmmc@netP$pathways
pathways_cd34 <- cellchat_cd34@netP$pathways

common_pathways <- intersect(pathways_bmmc, pathways_cd34)

cat("Significant pathways in BMMC:", length(pathways_bmmc), "\n")

## Significant pathways in BMMC: 23

cat("Significant pathways in CD34:", length(pathways_cd34), "\n")

## Significant pathways in CD34: 12

cat("Shared by both groups:       ", length(common_pathways), "\n\n")

## Shared by both groups:        10

print(common_pathways)

##  [1] "MIF"      "MHC-II"   "CD99"     "GALECTIN" "APP"      "RESISTIN"
##  [7] "SELL"     "MK"       "SELPLG"   "ITGB2"

cat("\nBMMC-specific:\n"); print(setdiff(pathways_bmmc, pathways_cd34))

## 
## BMMC-specific:

##  [1] "MHC-I"   "CLEC"    "ANNEXIN" "CD22"    "CD45"    "ICAM"    "LCK"    
##  [8] "SEMA4"   "IL16"    "CD23"    "FLT3"    "PARs"    "BAG"

cat("\nCD34-specific:\n"); print(setdiff(pathways_cd34, pathways_bmmc))

## 
## CD34-specific:

## [1] "NEGR" "CDH"

Interaction and Interaction Strength

# Aggregate number of interactions and interaction strength for each group.
# @net$count = inferred interaction counts, @net$weight = interaction strength.
interaction_summary <- data.frame(
  Group               = c("BMMC", "CD34"),
  Cell_Types          = c(nrow(cellchat_bmmc@net$count), nrow(cellchat_cd34@net$count)),
  Total_Interactions  = c(sum(cellchat_bmmc@net$count),  sum(cellchat_cd34@net$count)),
  Total_Strength      = c(sum(cellchat_bmmc@net$weight), sum(cellchat_cd34@net$weight)),
  Signalling_Pathways = c(length(pathways_bmmc), length(pathways_cd34))
)
print(interaction_summary)

##   Group Cell_Types Total_Interactions Total_Strength Signalling_Pathways
## 1  BMMC         16               2061      111.73806                  23
## 2  CD34         16                622       28.54916                  12

Per-pathway communication strength for the pathways shared by both groups. Strength is the summed communication probability over all sender-receiver pairs in that pathway.

pathway_strength <- function(cc, paths) {
  prob <- cc@netP$prob
  vapply(paths, function(p) {
    if (!is.null(prob) && p %in% dimnames(prob)[[3]]) sum(prob[, , p], na.rm = TRUE) else NA_real_
  }, numeric(1))
}

common_pathway_strength <- data.frame(
  Pathway       = common_pathways,
  BMMC_Strength = pathway_strength(cellchat_bmmc, common_pathways),
  CD34_Strength = pathway_strength(cellchat_cd34, common_pathways),
  row.names     = NULL
)
common_pathway_strength$Total <- with(common_pathway_strength, BMMC_Strength + CD34_Strength)
common_pathway_strength <- common_pathway_strength %>% arrange(desc(Total))

print(common_pathway_strength)

##     Pathway BMMC_Strength CD34_Strength     Total
## 1       MIF     20.778923     6.0525955 26.831518
## 2      CD99     11.589118     5.0553974 16.644516
## 3    MHC-II     13.115591     1.4963622 14.611953
## 4       APP      7.306862     4.8630772 12.169940
## 5  GALECTIN      7.683475     2.5053015 10.188777
## 6        MK      2.282923     5.1702259  7.453148
## 7      SELL      2.728488     2.9937032  5.722191
## 8  RESISTIN      3.472059     0.2227516  3.694811
## 9    SELPLG      1.707704     0.1771181  1.884822
## 10    ITGB2      1.263777     0.1359507  1.399727

11. Summary

We performed single-cell RNA sequencing on four bone marrow samples, including two Bone Marrow Mononuclear Cells (BMMC) and two CD34+ Enriched Bone Marrow Cells. The data was processed in Seurat, where quality control filtered out low-quality cells and genes. Doublets were removed with DoubletFinder, and Seurat’s integration method addressed batch effects. Dimensionality reduction via PCA and UMAP enabled clustering to identify cell populations, and cell types were annotated using SingleR and marker-based methods. Differential expression analysis highlighted key cell types, while pathway enrichment offered insights into biological processes. Monocle 3 was used for trajectory analysis focused on Common Lymphoid Progenitors (CLP), and CellChat identified cell-cell communication pathways.

Future improvements could involve SCTransform normalization, alternative clustering (e.g., Leiden), or using Scanpy and Harmony for flexibility and batch correction. Advanced machine learning could help identify rare cell types and additional pathways. Integrating multi-omics, such as spatial transcriptomics, ATAC-seq, and proteomics, would deepen insights into hematopoiesis, cellular interactions, and regulatory mechanisms.

Session Info

sessionInfo()

## R version 4.0.5 (2021-03-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /scratch/mdra00001/envs/single-cell/lib/libopenblasp-r0.3.27.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3          caret_6.0-93               
##  [3] lattice_0.20-45             pheatmap_1.0.12            
##  [5] Matrix_1.4-1                ROCR_1.0-11                
##  [7] KernSmooth_2.23-20          fields_14.1                
##  [9] viridis_0.6.2               viridisLite_0.4.1          
## [11] spam_2.9-1                  knitr_1.40                 
## [13] celldex_1.0.0               monocle3_1.0.0             
## [15] SeuratWrappers_0.3.0        SingleCellExperiment_1.12.0
## [17] CellChat_1.6.1              ggplot2_3.3.6              
## [19] igraph_1.3.4                enrichR_3.1                
## [21] SingleR_1.4.1               SummarizedExperiment_1.20.0
## [23] Biobase_2.50.0              GenomicRanges_1.42.0       
## [25] GenomeInfoDb_1.26.4         IRanges_2.24.1             
## [27] S4Vectors_0.28.1            BiocGenerics_0.36.0        
## [29] MatrixGenerics_1.2.1        matrixStats_0.62.0         
## [31] DoubletFinder_2.0.3         patchwork_1.1.2            
## [33] sp_1.5-0                    SeuratObject_4.1.1         
## [35] Seurat_4.0.1                spatstat.core_2.4-4        
## [37] rpart_4.1.16                nlme_3.1-159               
## [39] spatstat.random_2.2-0       spatstat.geom_2.4-0        
## [41] spatstat.data_2.2-0         dplyr_1.0.10               
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3                scattermore_0.8              
##   [3] ModelMetrics_1.2.2.2          coda_0.19-4                  
##   [5] tidyr_1.2.1                   bit64_4.0.5                  
##   [7] irlba_2.3.5                   DelayedArray_0.16.3          
##   [9] data.table_1.14.2             hardhat_1.2.0                
##  [11] RCurl_1.98-1.8                doParallel_1.0.17            
##  [13] generics_0.1.3                leidenbase_0.1.3             
##  [15] cowplot_1.1.1                 RSQLite_2.2.8                
##  [17] RANN_2.6.1                    proxy_0.4-27                 
##  [19] future_1.28.0                 bit_4.0.4                    
##  [21] lubridate_1.8.0               httpuv_1.6.6                 
##  [23] assertthat_0.2.1              gower_1.0.0                  
##  [25] xfun_0.33                     jquerylib_0.1.4              
##  [27] evaluate_0.16                 promises_1.2.0.1             
##  [29] fansi_1.0.3                   dbplyr_2.2.1                 
##  [31] DBI_1.1.3                     htmlwidgets_1.5.4            
##  [33] purrr_0.3.4                   ellipsis_0.3.2               
##  [35] RSpectra_0.16-1               ggpubr_0.4.0                 
##  [37] backports_1.4.1               gridBase_0.4-7               
##  [39] deldir_2.0-4                  sparseMatrixStats_1.2.1      
##  [41] vctrs_0.4.1                   ggalluvial_0.12.6            
##  [43] remotes_2.4.2                 Cairo_1.6-0                  
##  [45] abind_1.4-5                   cachem_1.0.6                 
##  [47] withr_2.5.0                   progressr_0.11.0             
##  [49] sctransform_0.3.4             sna_2.6                      
##  [51] goftest_1.2-3                 svglite_2.1.1                
##  [53] cluster_2.1.3                 ExperimentHub_1.16.0         
##  [55] cleanrmd_0.1.1                dotCall64_1.0-1              
##  [57] lazyeval_0.2.2                crayon_1.5.1                 
##  [59] recipes_1.0.1                 labeling_0.4.2               
##  [61] pkgconfig_2.0.3               nnet_7.3-17                  
##  [63] rlang_1.0.6                   globals_0.16.1               
##  [65] lifecycle_1.0.2               miniUI_0.1.1.1               
##  [67] registry_0.5-1                BiocFileCache_1.14.0         
##  [69] rsvd_1.0.5                    AnnotationHub_2.22.0         
##  [71] polyclip_1.10-0               lmtest_0.9-40                
##  [73] rngtools_1.5.2                carData_3.0-5                
##  [75] zoo_1.8-11                    ggridges_0.5.4               
##  [77] GlobalOptions_0.1.2           png_0.1-7                    
##  [79] rjson_0.2.21                  bitops_1.0-7                 
##  [81] pROC_1.18.0                   ggnetwork_0.5.10             
##  [83] blob_1.2.3                    DelayedMatrixStats_1.12.3    
##  [85] shape_1.4.6                   stringr_1.4.1                
##  [87] parallelly_1.32.1             rstatix_0.7.0                
##  [89] ggsignif_0.6.3                beachmat_2.6.4               
##  [91] scales_1.2.1                  memoise_2.0.1                
##  [93] magrittr_2.0.3                plyr_1.8.7                   
##  [95] ica_1.0-3                     zlibbioc_1.36.0              
##  [97] compiler_4.0.5                clue_0.3-60                  
##  [99] fitdistrplus_1.1-8            cli_3.4.1                    
## [101] XVector_0.30.0                listenv_0.8.0                
## [103] pbapply_1.5-0                 MASS_7.3-58.1                
## [105] mgcv_1.8-40                   tidyselect_1.1.2             
## [107] stringi_1.7.8                 highr_0.9                    
## [109] yaml_2.3.5                    BiocSingular_1.6.0           
## [111] ggrepel_0.9.1                 grid_4.0.5                   
## [113] sass_0.4.2                    tools_4.0.5                  
## [115] future.apply_1.9.1            circlize_0.4.15              
## [117] foreach_1.5.2                 gridExtra_2.3                
## [119] prodlim_2019.11.13            farver_2.1.1                 
## [121] Rtsne_0.16                    digest_0.6.29                
## [123] BiocManager_1.30.18           rgeos_0.5-9                  
## [125] lava_1.6.10                   FNN_1.1.3.1                  
## [127] shiny_1.7.2                   Rcpp_1.0.9                   
## [129] car_3.1-0                     broom_1.0.1                  
## [131] BiocVersion_3.12.0            later_1.2.0                  
## [133] RcppAnnoy_0.0.19              httr_1.4.4                   
## [135] AnnotationDbi_1.52.0          ComplexHeatmap_2.6.2         
## [137] colorspace_2.0-3              tensor_1.5                   
## [139] reticulate_1.30               splines_4.0.5                
## [141] uwot_0.1.14                   spatstat.utils_2.3-1         
## [143] plotly_4.10.0                 systemfonts_1.0.4            
## [145] xtable_1.8-4                  jsonlite_1.8.0               
## [147] timeDate_4021.104             ipred_0.9-13                 
## [149] R6_2.5.1                      pillar_1.8.1                 
## [151] htmltools_0.5.3               mime_0.12                    
## [153] NMF_0.28                      glue_1.6.2                   
## [155] fastmap_1.1.0                 BiocParallel_1.24.1          
## [157] BiocNeighbors_1.8.2           class_7.3-20                 
## [159] interactiveDisplayBase_1.28.0 codetools_0.2-18             
## [161] maps_3.4.0                    utf8_1.2.2                   
## [163] bslib_0.4.0                   spatstat.sparse_2.1-1        
## [165] tibble_3.1.8                  network_1.17.1               
## [167] curl_4.3.2                    leiden_0.4.3                 
## [169] limma_3.46.0                  survival_3.4-0               
## [171] rmarkdown_2.16                statnet.common_4.5.0         
## [173] munsell_0.5.0                 e1071_1.7-11                 
## [175] GetoptLong_1.0.5              GenomeInfoDbData_1.2.4       
## [177] iterators_1.0.14              reshape2_1.4.4               
## [179] gtable_0.3.1