Introduction

This practical session is designed to help you understand the importance of normalization, scaling, and regression in single-cell RNA-seq analysis. You will apply these steps using different methods and visualize the effects on the data.

Input data: Seurat object containing the filtered count matrix from the previous class. It’s called 05_TD3A_S5_Doublets.Filtered_12508.4035.RDS.

Output data: Seurat object after the normalization, scaling and regression. It’s called sobj_TD3A_normalized.rds.

Learning Objectives

  • Discuss why normalizing counts is necessary to compare cells
  • Be able to understand different normalization approaches
  • Be able to decide when to regress out a given variable, by evaluating the effects from any unwanted sources of variation

Motivation

Up to this point, we have filtered out low-quality cells, ambient RNA contamination, and doublets. The data is now a count matrix with cells as columns and genes as rows. These counts reflect the capture, reverse transcription, and sequencing steps of the scRNA-seq experiment, each introducing variability. This means the observed differences in gene expression may be influenced by sampling noise rather than true biological differences, resulting in uneven variance across the dataset.

Normalization addresses this by adjusting raw counts to minimize the effects of variable sampling, making the data more suitable for analysis, which often assumes uniform variance. Various normalization techniques exist, each tailored to support different downstream analyses.

A recent study by Ahlmann-Eltze and Huber (2023) benchmarked 22 normalization methods. It’s essential to select the normalization method based on the specific analysis goals. However, for this practical we will just explore the most well-known way to doing single-cell normalization.

Setup

We first load the packages that are needed :

library(Seurat)
library(ggplot2)
library(patchwork)
library(scater)
library(dplyr)

We then import the data that have been previously filtered (reminder : we removed low quality cells, removed ambient RNA and scored doublets) :

sobj <- readRDS("input/05_TD3A_S5_Doublets.Filtered_12508.4035.RDS")
sobj
## An object of class Seurat 
## 12508 features across 4035 samples within 1 assay 
## Active assay: RNA (12508 features, 0 variable features)
##  2 layers present: counts, data

For the sake of this practical session, we will use a smaller dataset which contains only 1000 cells. In order to do this, we can use the Seurat objet as a dataframe, where rows are genes, and columns are cells. In R language, this corresponds to : sobj[genes,cells].

We will take the first 1000 cells, by using the range 1:1000, like this :

# Select the cells to keep 
sobj_tiny <- sobj[,1:1000]
# Check that the object contains 1000 cells
sobj_tiny
## An object of class Seurat 
## 12508 features across 1000 samples within 1 assay 
## Active assay: RNA (12508 features, 0 variable features)
##  2 layers present: counts, data

Exercices 1 : explore NormalizeData function

Let’s start by exploring the NormalizeData function.

What is the input of this function ?

Answers

The input is your Seurat object (for us : sobj_tiny)

How many normalization methods are they ?

Can you spot which is the default one ?

Answers

There are 3 normalization approaches included in Seurat : LogNormalize, CLR and RC. The default one is LogNormalize.

Run the NormalizeData function on the sobj_tiny

You can now explore your new object : where are the normalized values stored ? Can you print the values of the first 4 genes in the 10 first cells ?

Answers 
sobj_tiny <- NormalizeData(sobj_tiny)

# The normalized values are stored in the slot 'data' : 
sobj_tiny@assays$RNA$data[1:4, 1:10]
## 4 x 10 sparse Matrix of class "dgCMatrix"
##                                                                          
## Mrpl15  . 2.365796 .        .        .        .        .        0.6103322
## Lypla1  . .        .        1.844929 .        1.627274 1.442936 0.6103322
## Gm37988 . .        .        .        .        .        .        .        
## Tcea1   . .        1.781829 .        1.708152 1.627274 2.010388 .        
##                          
## Mrpl15  .        .       
## Lypla1  .        1.052803
## Gm37988 .        .       
## Tcea1   1.701157 1.052803
# You can compare it to the slot 'count' (raw count) :
sobj_tiny@assays$RNA$count[1:4, 1:10]
## 4 x 10 sparse Matrix of class "dgCMatrix"
##                            
## Mrpl15  . 2 . . . . . 1 . .
## Lypla1  . . . 1 . 1 1 1 . 1
## Gm37988 . . . . . . . . . .
## Tcea1   . . 1 . 1 1 2 . 1 1

Exercises 2 : let’s test some different normalization strategies

Method 1: LogNormalize

LogNormalize method or “scaling normalization” is the simplest and most commonly used class of normalization approaches. This involves dividing all counts for each cell by a cell-specific scaling factor, often called a “size factor” (Anders and Huber 2010). The assumption here is that any cell-specific bias (e.g., in capture or amplification efficiency) affects all genes equally via scaling of the expected mean count for that cell. The size factor for each cell represents the estimate of the relative bias in that cell, so division of its counts by its size factor should remove that bias. The resulting “normalized expression values” can then be used for downstream analyses such as clustering and dimensionality reduction. (From https://bioconductor.org/books/3.14/OSCA.basic/normalization.html)

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

We can now visualize the effect of LogNormalize with some violin plots :

VlnPlot(sobj_tiny,
        features = c("Srm", "Apoe", "Top2a"),
        layer = "data")

Can you explain to your neighbor what does this plot mean ?

Violin plots are sometimes hard to read. Here is a function to build box plots out of Seurat objects :

BoxPlotFromSeurat <- function(seurat_obj, features, log = TRUE, layer = "data", plot.title = "Box Plot") {
  # Extract the data for the features of interest
  df <- Seurat::FetchData(seurat_obj, vars = features, layer = layer)
  
  # Add metadata for grouping if group.by is specified
  df <- cbind(df, seurat_obj[[]])
  
  # Reshape the data for ggplot2
  df_long <- tidyr::gather(df, key = "feature", value = "expression", group)
  df_long <- df %>% group_by_(group.by) %>%
    tidyr::pivot_longer(cols = features)

  # Create box plots using ggplot2
  if (log == TRUE) {
    p <- ggplot(df_long, ggplot2::aes(x = name, y = log10(value))) +
    geom_boxplot(outlier.shape = NA) +
    theme_minimal() +
    labs(x = "Feature", y = "Expression Level (log10)") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  } else {
    p <- ggplot(df_long, ggplot2::aes(x = name, y = value)) +
    geom_boxplot(outlier.shape = NA) +
    theme_minimal() +
    labs(x = "Feature", y = "Expression Level") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  }
  return(p)
}

# Example usage
# BoxPlotFromSeurat(sobj_tiny, features = c("Srm", "Apoe", "Top2a"))

Can you do the same Violin plots but with the raw counts ?

You can try also with the boxplots.

Answers 
VlnPlot(sobj_tiny,
        features = c("Srm", "Apoe", "Top2a"), 
        layer = "count")

Can you plot the distribution of the values ?

Try to do an histogram, with 1) the raw counts values and 2) the log normalized values. You can use the R hist() function. What do you observe ? You will need to extract your data with a command such as :

raw_counts <- as.vector(as.matrix(sobj_tiny@assays$RNA$count))
Answers 
# Raw values
raw_counts <- as.vector(as.matrix(sobj_tiny@assays$RNA$count))
hist(raw_counts)

# Log Norm values
logNorm_counts <- as.vector(as.matrix(sobj_tiny@assays$RNA$data))
hist(logNorm_counts)

# If you want to visualize the counts without the zero values, add these lines :
raw_counts <- raw_counts[raw_counts > 0]
logNorm_counts <- logNorm_counts[logNorm_counts > 0]

# And then again : 
hist(raw_counts)

hist(logNorm_counts)

# You can play with the "breaks" parameters to have a better visualization :
hist(raw_counts, breaks = 100)

hist(logNorm_counts, breaks = 100)

We can also plot the distribution of all genes in our dataset, and visualize it with some boxplots.

# Prepare the data 
raw_df <- data.frame(Expression = raw_counts, Type = "Raw Counts")
norm_df <- data.frame(Expression = logNorm_counts, Type = "Normalized Counts")
plot_data <- rbind(raw_df, norm_df)

# Create the boxplot
boxplot_plot <- ggplot(plot_data, aes(x = Type, y = Expression, fill = Type)) +
  geom_boxplot(outlier.shape = NA) +
  scale_y_log10() +
  ggtitle("Boxplots of Gene Expression Across All Genes (Before and After Normalization)") +
  xlab("Normalization Status") +
  ylab("Expression (Log10 Scale)") +
  theme_minimal()

# Print the boxplot
print(boxplot_plot)

What can you observe when comparing the 2 boxplots ?

How do we know that we effectively normalized the data by reading the boxplots ?

Until now, we explored the variability among the genes count values. We can do the same with the cells :

# Prepare the data 
df <- as.data.frame(sobj_tiny[[]])
df$cell_index <- rownames(df)

# Plots before and after normalization
ggplot(df, aes(x = cell_index, y = nCount_RNA)) + geom_bar(stat="identity")

ggplot(df, aes(x = cell_index, y = log10_nCount_RNA)) + geom_bar(stat="identity")

What can you conclude about the sequencing depths in the cells ?

Finally, we want to see the impact of the normalization on a highly variable gene. We chose to focus on “Tmsb10”.

Can you plot the total number of Tmsb10 expression versus the total number of UMI ?

You can use the FeatureScatter function for this.

Answers 
FeatureScatter(sobj_tiny,
               feature1 = "nCount_RNA",
               feature2 = "Tmsb10",
               slot = "count")

Compare it with the normalized values

What can you observe ?

Answers 
FeatureScatter(sobj_tiny,
               feature1 = "nCount_RNA",
               feature2 = "Tmsb10",
               slot = "data")

Discussion : LogNormalize reduces the impact of sequencing depth differences but may not handle highly variable genes effectively.

Method 2: SCTransform

SCTransform is a normalization method designed for single-cell RNA-seq data. It uses a variance-stabilizing transformation (VST) to model and remove technical noise, making it particularly effective for datasets with high variability and technical noise, such as dropouts and differences in sequencing depth. SCTransform is available in the Seurat package and is often used as an alternative to traditional log normalization.

Key Features of SCTransform :

  • Variance stabilization: SCTransform aims to stabilize the variance across genes, making it less affected by highly variable or highly expressed genes.

  • Data modeling: it uses a generalized linear model (GLM) to model the counts as a function of sequencing depth, accounting for technical noise.

  • No manual scaling needed: SCTransform automatically standardizes the data, so additional scaling is not required before PCA.

# Apply SCTransform
sobj_tiny <- SCTransform(sobj_tiny,
                         verbose = FALSE)

# Visualize the effect of SCTransform
VlnPlot(sobj_tiny,
        features = c("Srm", "Apoe", "Top2a"),
        log = TRUE)

Try to answers all the same questions we asked above for SCTransform

Compare the results between log normalization and SCTransform

Here is an example of the results on the Violin Plots :

# Normalize data using LogNormalize (for comparison)
sobj_tiny_log <- NormalizeData(sobj_tiny, normalization.method = "LogNormalize")
sobj_tiny_log <- ScaleData(sobj_tiny_log)

# Normalize data using SCTransform
sobj_tiny_sct <- SCTransform(sobj_tiny, verbose = FALSE)

# Compare the distributions of gene expression using violin plots
vln_log <- VlnPlot(sobj_tiny_log, features = c("Srm", "Apoe", "Top2a"), log = TRUE)
vln_sct <- VlnPlot(sobj_tiny_sct, features = c("Srm", "Apoe", "Top2a"))

# Display the plots side by side
vln_log

vln_sct

Discussion: SCTransform handles variability better, making it suitable for noisy datasets.

Exercises 3 : Scaling

Until now, we have talked a lot about the normalization but what about the scaling ??? Well, it is just an additional step that we use to standardize the expression values across all cells.

A few words on scaling

In single-cell RNA-seq, scaling typically refers to standardizing the expression of each gene across all cells (not across genes !). The goal is to make the expression of each gene comparable across all cells, ensuring that highly expressed genes with large variance don’t dominate the analysis.

It typically involves transforming each gene’s expression values so that they have a mean of 0 and a standard deviation of 1 (z-score transformation). This ensures that each gene contributes equally to downstream analyses, such as PCA, clustering, and UMAP.

If you want a comparison, it’s like converting temperatures from Celsius and Fahrenheit to a common scale so we can compare them easily. Moreover, we scale to ensure that highly expressed genes don’t dominate the analysis (e.g., in PCA or clustering). Without scaling, a small set of genes with very high variance can bias the results.

Please note that not all genes need to be scaled (e.g., marker genes of interest). Be careful also that scaling should be done after normalization, not before.

How it works ?

For each gene, we take the expression values across all cells and :

1- Center them by subtracting the mean expression of the gene

2- Scale them by dividing by the standard deviation of the gene

\[z = \frac{x - \text{mean}}{\text{standard deviation}}\]

Where, \(x\) is the expression value for a particular gene in a specific cell. Mean and standard deviation are calculated across all cells for that gene.

What it is not ? Scaling is not about making the expression of each cell comparable across all genes. That’s what normalization (e.g., LogNormalize, SCTransform) does : adjusting for sequencing depth differences across cells.

How to scale the data in Seurat ? With the function ScaleData.

Start by exploring the ScaleData function.

What are the arguments you think will be useful for us ?

Scale the data of the tiny Seurat object

Answers 
sobj_tiny@active.assay <- "RNA"
sobj_tiny <- ScaleData(sobj_tiny)

Store the new scaled data in a matrix called scaled_data

Answers 
scaled_data <- as.matrix(sobj_tiny@assays$RNA@layers$scale.data)

After scaling, the expression values for each gene will have a mean of 0 and a standard deviation of 1.

Can you check that you also find those values ?

You will probably need the R functions rowMeans and rowSds. If you don’t know them, have a look at their help page.

Answers 
# Calculate the mean and standard deviation for each gene without using apply()
gene_means <- rowMeans(scaled_data)
gene_sds <- rowSds(scaled_data)

# Check the results
summary(gene_means)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.02466 -0.00290  0.00000 -0.00453  0.00000  0.00000
summary(gene_sds)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.9327  1.0000  0.8859  1.0000  1.0000

Advanced : use scaled data for differential expression analysis

Let’s see what happens when we use scaled data for DE analysis. You can run the analyses with the functions FindMarkers or FindAllMarkers. For this exercise, it will be more useful if you first run a clustering analyses in order to get identified groups of cells.

Example :

sobj_tiny <- NormalizeData(sobj_tiny)
sobj_tiny <- FindVariableFeatures(sobj_tiny)
sobj_tiny <- ScaleData(sobj_tiny)
sobj_tiny <- RunPCA(sobj_tiny, dims = 1:20)
sobj_tiny <- FindNeighbors(sobj_tiny, dims = 1:20)
sobj_tiny <- FindClusters(sobj_tiny, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1000
## Number of edges: 41202
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7202
## Number of communities: 7
## Elapsed time: 0 seconds
sobj_tiny <- RunUMAP(sobj_tiny, dims = 1:20)
DimPlot(sobj_tiny, reduction = "umap")

You can also run a first differential expression analyses to be able to compare after your results :

# Perform DE analyses
sobj_tiny.markers <- FindAllMarkers(sobj_tiny, only.pos = TRUE)

# Extract top 10 genes per cluster
top10 <- sobj_tiny.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup()

# Display heatmap of DE results
DoHeatmap(sobj_tiny, features = top10$gene) + NoLegend() + theme(text = element_text(size = 7))
## Warning in DoHeatmap(sobj_tiny, features = top10$gene): The following features
## were omitted as they were not found in the scale.data slot for the RNA assay:
## Gm11725, Nln, Notch3, Dennd3, Pdk1, Slc43a1, Tle2, Smoc1, Cep85, Tdrd5,
## BB031773, Plxdc1, Patj

# Display top DE results
head(top10, n = 30)
## # A tibble: 30 × 7
##            p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
##            <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 0.00000000185       1.68 0.242 0.114 0.0000231 0       Patj    
##  2 0.00000000657       1.05 0.374 0.247 0.0000821 0       Plxdc1  
##  3 0.0000000760        1.60 0.206 0.096 0.000951  0       Slc6a19 
##  4 0.000000922         1.48 0.16  0.068 0.0115    0       BB031773
##  5 0.00000100          1.41 0.204 0.1   0.0125    0       Tdrd5   
##  6 0.00000290          1.04 0.364 0.278 0.0362    0       Cep85   
##  7 0.00000474          2.02 0.107 0.036 0.0593    0       Smoc1   
##  8 0.00000877          1.76 0.115 0.043 0.110     0       Tle2    
##  9 0.0000318           1.06 0.318 0.236 0.398     0       Slc43a1 
## 10 0.0000360           1.23 0.262 0.181 0.451     0       Pdk1    
## # ℹ 20 more rows
Answers 
# Attempt DE analysis using scaled data
sobj_tiny_scaled.markers <- FindAllMarkers(sobj_tiny, only.pos = TRUE, assay = "RNA", slot = "scale.data")
## Warning: No DE genes identified
## Warning: The following tests were not performed:
## Warning: When testing 0 versus all:
##  subscript out of bounds
## Warning: When testing 1 versus all:
##  subscript out of bounds
## Warning: When testing 2 versus all:
##  subscript out of bounds
## Warning: When testing 3 versus all:
##  subscript out of bounds
## Warning: When testing 4 versus all:
##  subscript out of bounds
## Warning: When testing 5 versus all:
##  subscript out of bounds
## Warning: When testing 6 versus all:
##  subscript out of bounds
# Not working !! Scaled data is designed for clustering, dimensionality reduction, and visualizations.

Can you try also with raw counts (not normalized) ? What do you observe ?

Answers 
markers <- FindAllMarkers(sobj_tiny, only.pos = TRUE, slot = "counts")

top10 <- markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup()

# Display heatmap of DE results
DoHeatmap(sobj_tiny, features = top10$gene) + NoLegend() + theme(text = element_text(size = 7))

# Display top DE results
head(markers)
##                p_val avg_log2FC pct.1 pct.2    p_val_adj cluster    gene
## mt-Co2  1.493196e-21  0.1506044 0.997 0.995 1.867690e-17       0  mt-Co2
## Eif1    4.006820e-18  0.1422588 0.977 0.992 5.011731e-14       0    Eif1
## Cd3d    1.441059e-17  0.1042215 0.936 0.979 1.802476e-13       0    Cd3d
## Serf2   1.783625e-15  0.1231241 0.952 0.980 2.230958e-11       0   Serf2
## Laptm5  2.012218e-15  0.1146634 0.944 0.959 2.516882e-11       0  Laptm5
## mt-Atp8 3.073331e-15  0.1966736 1.000 0.995 3.844122e-11       0 mt-Atp8

Discussion : * Why might the DE results be non-significant or unexpected ? * What information is lost during scaling that is critical for DE analysis ?

Advanced: change the values of the features option

This value is set by default to 2000. It represents the number of highly variable genes to take into account. You can explore the impact it has on the PCA. Hint : use the RunPCA function from Seurat.

To answer this question, you will need a function to reinitialize your Seurat object as follows :

# Function to reinitialize the Seurat object
reinitialize_sobj <- function(seurat_obj) {
  
  # Copy the object
  seurat_copy <- seurat_obj
  # Remove variable features
  seurat_copy@assays$RNA@meta.data <- data.frame()
  # Remove normalized and scales data
  seurat_copy@assays$RNA@layers$data <- NULL
  seurat_copy@assays$RNA@layers$scale.data <- NULL
  # Remove PCA and UMAP
  seurat_copy@reductions <- list()
  # Remove SCT
  seurat_copy@assays$SCT <- NULL
  # Set active assay to RNA
  seurat_copy@active.assay <- "RNA"
  
  return(seurat_copy)
}
Answers 
# Function to perform analysis with a specified number of variable features
perform_analysis <- function(seurat_obj, num_features) {
  # Create a new Seurat object from raw counts
  seurat_copy <- reinitialize_sobj(seurat_obj)
  
  # Normalize the data
  seurat_copy <- NormalizeData(seurat_copy)
  
  # Set the variable features (up to the specified number)
  seurat_copy <- FindVariableFeatures(seurat_copy,
                                     selection.method = "vst",
                                     nfeatures = num_features)
  
  # Identify the 10 most highly variable genes
  top_features <- head(VariableFeatures(seurat_copy), 20)
  print(top_features)
  
  # Scale the data, and run PCA
  seurat_copy <- ScaleData(seurat_copy)
  seurat_copy <- RunPCA(seurat_copy,
                        features = VariableFeatures(seurat_copy),
                        verbose = FALSE)
  
  # Generate PCA plot
  pca_plot <- DimPlot(seurat_copy, reduction = "pca") + ggtitle(paste(num_features, "Variable Features"))
  
  # Generate DimHeatmap for the first two principal components
  DimHeatmap(seurat_copy,
             dims = 1:2,
             balanced = TRUE,
             reduction = "pca")
  
  # Cluster the cells and run UMAP
  seurat_copy <- FindNeighbors(seurat_copy, dims = 1:20)
  seurat_copy <- FindClusters(seurat_copy, resolution = 0.5)
  seurat_copy <- RunUMAP(seurat_copy, dims = 1:20)
  umap_plot <- DimPlot(seurat_copy, reduction = "umap")
  
  # Find markers
  seurat_copy.markers <- FindAllMarkers(seurat_copy, only.pos = TRUE)
  top10 <- seurat_copy.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup()
  heatmap_markers_plot <- DoHeatmap(seurat_copy, features = top10$gene) + NoLegend()
  
  # Generate Variable Feature Plot with labeled top 20 features
  vf_plot <- VariableFeaturePlot(seurat_copy)
  vf_plot <- LabelPoints(plot = vf_plot, points = top_features, repel = TRUE) +
    ggtitle(paste("Top Variable Features:", num_features))
  
  # Print your new Seurat object
  print(seurat_copy)
  
  # Return a list of plots
  return(list(
    SeuratObject = seurat_copy,
    PCA = pca_plot,
    VariableFeatures = vf_plot,
    UMAP_plot = umap_plot,
    Heatmap_markers = heatmap_markers_plot))
}

# Analyze with different numbers of variable features
results_100 <- perform_analysis(sobj_tiny, 100)
##  [1] "Hist1h2ao" "Pclaf"     "Hist1h2ap" "Rrm2"      "Top2a"     "Hist1h2ae"
##  [7] "Hmgb2"     "Ms4a4b"    "Tuba1b"    "Itm2a"     "Birc5"     "Trbv16"   
## [13] "Ifi27l2a"  "Trbv29"    "Trbv15"    "Trbv5"     "H2-K1"     "Nkg7"     
## [19] "Ube2c"     "Hist1h1b"
## Warning in irlba(A = t(x = object), nv = npcs, ...): You're computing too large
## a percentage of total singular values, use a standard svd instead.
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1000
## Number of edges: 34398
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8030
## Number of communities: 14
## Elapsed time: 0 seconds
## Warning in DoHeatmap(seurat_copy, features = top10$gene): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: Tyrobp, Vcam1, Ear2, C1qb, Lgals3, C1qc, Cd74, Fcer1g, C1qa, Rflnb, Cd7,
## Gimap4, H2-Q7, Gm2682, S1pr1, Itgae, Tmc4, Gm20517, Trdj1, Ephb6, Trbj2-1,
## Cyp11a1, Trbj1-5, Traj2, Gm49355, Sdr42e1, Comp, Gabrd, Sdcbp2, Paqr8,
## C920008G01Rik, E2f8, Clspn, Rad51, Tmie, Pcsk1, Grik1, Gm30373, Zc2hc1a, Fbln1,
## Trbv13-2, Lypd1, Cd69, Cd40lg, Tigit, Nab2, Egr2, 0610005C13Rik, Sbsn, Uckl1os,
## 1700029I15Rik, Catspere2, 4732491K20Rik, Gm29019, Trav3-4, Gm44731, Zfp316,
## Gm43062, Gm15417, Tyro3, Gm38037, Myo5b, Slc9a3r2, Slit3, Trbj2-2, Gm17349,
## Gm34552, Gm30906, Trav6-2, Ubap1l, A830035O19Rik, Gm12840, Isl2, Star, Eif2s3y,
## Gm26526, Trav6-4, Mgam, Trim80, Cabyr, Tes, Casp9, Pam, Gm45623, Gm13920,
## 2810414N06Rik, Gm13919, Tpbgl, Trbj2-4, Fbln2, Ddx3y, 5830405F06Rik, Gm33104,
## Trbv24, Pde4dip, Tle2, Trbv23, Trbv21

## An object of class Seurat 
## 12508 features across 1000 samples within 1 assay 
## Active assay: RNA (12508 features, 100 variable features)
##  3 layers present: counts, data, scale.data
##  2 dimensional reductions calculated: pca, umap
results_2000 <- perform_analysis(sobj_tiny, 2000)
##  [1] "Hist1h2ao" "Pclaf"     "Hist1h2ap" "Rrm2"      "Top2a"     "Hist1h2ae"
##  [7] "Hmgb2"     "Ms4a4b"    "Tuba1b"    "Itm2a"     "Birc5"     "Trbv16"   
## [13] "Ifi27l2a"  "Trbv29"    "Trbv15"    "Trbv5"     "H2-K1"     "Nkg7"     
## [19] "Ube2c"     "Hist1h1b"
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1000
## Number of edges: 41202
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7202
## Number of communities: 7
## Elapsed time: 0 seconds
## Warning in DoHeatmap(seurat_copy, features = top10$gene): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: Gm11725, Nln, Notch3, Dennd3, Pdk1, Slc43a1, Tle2, Smoc1, Cep85, Tdrd5,
## BB031773, Plxdc1, Patj

## An object of class Seurat 
## 12508 features across 1000 samples within 1 assay 
## Active assay: RNA (12508 features, 2000 variable features)
##  3 layers present: counts, data, scale.data
##  2 dimensional reductions calculated: pca, umap
results_5000 <- perform_analysis(sobj_tiny, 5000)
##  [1] "Hist1h2ao" "Pclaf"     "Hist1h2ap" "Rrm2"      "Top2a"     "Hist1h2ae"
##  [7] "Hmgb2"     "Ms4a4b"    "Tuba1b"    "Itm2a"     "Birc5"     "Trbv16"   
## [13] "Ifi27l2a"  "Trbv29"    "Trbv15"    "Trbv5"     "H2-K1"     "Nkg7"     
## [19] "Ube2c"     "Hist1h1b"
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1000
## Number of edges: 40708
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7142
## Number of communities: 6
## Elapsed time: 0 seconds
## Warning in DoHeatmap(seurat_copy, features = top10$gene): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: St8sia1, Ccnd3, Tiparp, Gm44175, Cep85, Plxdc1, Gtf2h4, Arl5c

## An object of class Seurat 
## 12508 features across 1000 samples within 1 assay 
## Active assay: RNA (12508 features, 5000 variable features)
##  3 layers present: counts, data, scale.data
##  2 dimensional reductions calculated: pca, umap
# Display PCA plots for comparison
print(results_100$PCA)

print(results_2000$PCA)

print(results_5000$PCA)

# Display Variable Feature Plots
print(results_100$VariableFeatures)
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: ggrepel: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

print(results_2000$VariableFeatures)
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## ggrepel: 4 unlabeled data points (too many overlaps). Consider increasing max.overlaps

print(results_5000$VariableFeatures)
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning: ggrepel: 2 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

# Display UMAPs
print(results_100$UMAP_plot)

print(results_2000$UMAP_plot)