A stripped down version of Example #1, modified to allow key parameters to be easily changed in order to look at the stability of the final results.
Tests:
nFeature_RNA is cut at 2100 rather than the default of 2500 which seems to include one outlier at ca. 2400
Quite similar, 9 communities, UMAP just a little different but the same messages from data
nFeature_RNA is cut at 3000
Quite similar, 9 communities, UMAP just a little different
nFeature_RNA is cut at 4000 (not cut)
Now interesting
9 communities are generated but the cluster integer order is different and so labels need to be readjusted
Cells NK cluster is extended
Seems like there are some better gene expression markers than IL7R e.g.s LEF1 and PRKCQ-AS1
Seurat v3 and v5 UMAP have different default labels - v3 uppercase and v5 are lowercase. Useful for identification of versions but also the UMAP cluster rotations in 2D are different, together with some cluster separation and cluster shape aspects. Suspect Seurat is simply using a refined version of UMAP in v5.0.0.
ChangeLog
Updated to Seurat v5
Includes
Environment Load and Check
this code section is packaged as an include for reuse across all examples
it uses the HTML details tag directly to wrap code blocks and output as drop-down sections
Show Environment
library(patchwork)library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
library(Seurat)
Loading required package: SeuratObject
Loading required package: sp
'SeuratObject' was built under R 4.3.0 but the current version is
4.3.2; it is recomended that you reinstall 'SeuratObject' as the ABI
for R may have changed
Attaching package: 'SeuratObject'
The following object is masked from 'package:base':
intersect
Is there a typo in the message above? Application Programming Interface, API != ABI
Rolling back to R 4.3.0 was not possible with the current version of Seurat
the indication was that Seurat requires a version of base Matrix that is not present in R 4.3.0
# which Seurat?packageVersion("Seurat")
[1] '5.0.0'
# which R?version[['version.string']]
[1] "R version 4.3.2 (2023-10-31)"
# presto was installed # For a (much!) faster implementation of the # Wilcoxon Rank Sum TestpackageVersion('presto')
[1] '1.0.0'
# check python is available via reticulateimport sysprint(sys.version.split(" ")[0])
3.12.0
# shell checkpython3-V
Python 3.12.0
# shell checkquarto-v
1.4.489
Functions
Show Functions
# Useful for code development.# Save the object at a point and reload it into the R console # i.e. for developing alternative reports # without having to run the pipeline right from the start# which can be slow## NB: Files produced by saveRDS (or serialized to a file connection) # are not suitable as an interchange format between machines# # For that use hdf5 or transfer data and code to reproduce the result saveRDS_overwrite <-function(file_path) {if (file.exists(file_path)) {file.remove(file_path) } saveRDS(pbmc, file = file_path)}
n_pcs_chosen <-10# clusters# create the K-nearest neighbor (KNN) graph of cellspbmc <-FindNeighbors(pbmc, dims =1:n_pcs_chosen)
Computing nearest neighbor graph
Computing SNN
# define clusters according to resolutionpbmc <-FindClusters(pbmc, resolution =0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 2643
Number of edges: 96640
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8753
Number of communities: 9
Elapsed time: 0 seconds
# Look at cluster IDs of the first 5 cells# In this case we have 9 levels (0 - 8)# The structure is the relation between cell barcode and the cluster (community) head(Idents(pbmc), 5)
# Each cell that survived filtering above is represented length(Idents(pbmc))
[1] 2643
pbmc
An object of class Seurat
13714 features across 2643 samples within 1 assay
Active assay: RNA (13714 features, 2000 variable features)
3 layers present: counts, data, scale.data
1 dimensional reduction calculated: pca
# UMAPpbmc <-RunUMAP(pbmc, dims =1:n_pcs_chosen)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
16:26:11 UMAP embedding parameters a = 0.9922 b = 1.112
16:26:11 Read 2643 rows and found 10 numeric columns
16:26:11 Using Annoy for neighbor search, n_neighbors = 30
16:26:11 Building Annoy index with metric = cosine, n_trees = 50
# save the objectfile_path <-paste0("./seurat_object_checkpoints/pbmc_sw1",EGN,".rds")saveRDS_overwrite(file_path)# to restore# pbmc <- readRDS(file_path)# find markers for every cluster compared to all remaining cells, # report only the positive onespbmc.markers <-FindAllMarkers(pbmc, only.pos =TRUE, min.pct =0.25, logfc.threshold =0.25)
# heatmaps # note that 'wt' specifies the variable to use for ordering # we get the best markers in terms of size effectpbmc.markers %>%group_by(cluster) %>%top_n(n =5, wt = avg_log2FC) -> top_nhead(top_n,10)
# heat map shows that cluster 1 and 2 are not easily distiguished# by just a few genes others are DoHeatmap(pbmc, features = top_n$gene)
# + NoLegend()# consider the canonical cluster 0 markers# "IL7R" is not particularly good marker for cluster 0 VlnPlot(pbmc, features =c("IL7R", "CCR7"))
FeaturePlot(pbmc, features =c("IL7R", "CCR7"))
# Try the top 5# Seems like there are some better markers than IL7R e.g.s LEF1 and PRKCQ-AS1VlnPlot(pbmc, features =c("LDHB", "CCR7","LEF1","PRKCQ-AS1","LDLRAP1"))
FeaturePlot(pbmc, features =c("LDHB", "CCR7","LEF1","PRKCQ-AS1","LDLRAP1"))
Compare results
# cell type to cluster#new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", # "B", "CD8 T", "FCGR3A+ Mono",# "NK", "DC", "Platelet")# the unlabelled clustersDimPlot(pbmc, reduction ="umap")
