I made a single chunk
for this time
# Install packages (only required once)
#install.packages(c("dplyr", "Seurat", "patchwork"))
# Load packages
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
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
## v5.6 adds significant improvements in speed and efficiency, particularly for large datasets.
## For more information, see https://satijalab.org/seurat/articles/announcements.
library(patchwork)
library(ggplot2)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "I:/OxfordBioDiscovery2025/2026/Single cell RNAseq/2. Single Cell RNAseq")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
pbmc
## An object of class Seurat
## 13714 features across 2700 samples within 1 assay
## Active assay: RNA (13714 features, 0 variable features)
## 1 layer present: counts
head(pbmc@meta.data)
## orig.ident nCount_RNA nFeature_RNA
## AAACATACAACCAC-1 pbmc3k 2419 779
## AAACATTGAGCTAC-1 pbmc3k 4903 1352
## AAACATTGATCAGC-1 pbmc3k 3147 1129
## AAACCGTGCTTCCG-1 pbmc3k 2639 960
## AAACCGTGTATGCG-1 pbmc3k 980 521
## AAACGCACTGGTAC-1 pbmc3k 2163 781
# Standard pre-processing workflow
#QC and selecting cells for further analysis
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
# FeatureScatter is typically used to visualize feature-feature relationships,
# But can be used for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
# From both these figures, we found that nCount_RNA range should be >200 and <2000
# percent.mt should be <5
# This selection will vary from experimental data, also, person to person.
pbmc_qc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
#Normalizing the data
pbmc_norm <- NormalizeData(pbmc_qc, normalization.method = "LogNormalize", scale.factor = 10000)
## Normalizing layer: counts
# Identification of highly variable features (feature selection)
# We next identify features with high cell-to-cell variation, meaning they are highly expressed in some cells and lowly expressed in others.
# Focusing on these genes in downstream analyses helps highlight biological signals in single-cell datasets.
pbmc_hv <- FindVariableFeatures(pbmc_norm, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
pbmc_hv
## An object of class Seurat
## 13714 features across 2638 samples within 1 assay
## Active assay: RNA (13714 features, 2000 variable features)
## 2 layers present: counts, data
# See, we have same number of cells
# But genes (features) has been reduced to 2000
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc_hv), 10)
top10
## [1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1"
## [9] "GNG11" "S100A8"
# Visualize variable features data
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc_hv)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
combined_plot <-plot1 + plot2
# Save high-resolution image
ggsave(
filename = "PBMC_Variable_Features.png",
plot = combined_plot,
width = 14,
height = 6,
dpi = 300
)
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
# Scaling the data
# Next, linear transformation (‘scaling’) prior to PCA.
# The ScaleData() function: Shifts the expression of each gene, so that the mean expression across cells is 0
# Scales the expression of each gene, so that the variance across cells is 1
# This step gives equal weight in downstream analyses, so that highly-expressed genes do not dominate
all.genes <- rownames(pbmc_hv)
pbmc_scale <- ScaleData(pbmc_hv, features = all.genes)
## Centering and scaling data matrix
# Perform linear dimensional reduction PCA
pbmc_pca <- RunPCA(pbmc_scale, features = VariableFeatures(object = pbmc_scale))
## PC_ 1
## Positive: CST3, TYROBP, LST1, AIF1, FTL, FTH1, LYZ, FCN1, S100A9, TYMP
## FCER1G, CFD, LGALS1, S100A8, CTSS, LGALS2, SERPINA1, IFITM3, SPI1, CFP
## PSAP, IFI30, SAT1, COTL1, S100A11, NPC2, GRN, LGALS3, GSTP1, PYCARD
## Negative: MALAT1, LTB, IL32, IL7R, CD2, B2M, ACAP1, CD27, STK17A, CTSW
## CD247, GIMAP5, AQP3, CCL5, SELL, TRAF3IP3, GZMA, MAL, CST7, ITM2A
## MYC, GIMAP7, HOPX, BEX2, LDLRAP1, GZMK, ETS1, ZAP70, TNFAIP8, RIC3
## PC_ 2
## Positive: NKG7, PRF1, CST7, GZMB, GZMA, FGFBP2, CTSW, GNLY, B2M, SPON2
## CCL4, GZMH, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX
## TTC38, APMAP, CTSC, S100A4, IGFBP7, ANXA1, ID2, IL32, XCL1, RHOC
## Negative: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1, HLA-DRA, LINC00926, CD79B, HLA-DRB1, CD74
## HLA-DMA, HLA-DPB1, HLA-DQA2, CD37, HLA-DRB5, HLA-DMB, HLA-DPA1, FCRLA, HVCN1, LTB
## BLNK, P2RX5, IGLL5, IRF8, SWAP70, ARHGAP24, FCGR2B, SMIM14, PPP1R14A, C16orf74
## PC_ 3
## Positive: PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU
## HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, PTCRA, CA2, ACRBP, MMD, TREML1
## NGFRAP1, F13A1, SEPT5, RUFY1, TSC22D1, MPP1, CMTM5, RP11-367G6.3, MYL9, GP1BA
## Negative: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1, HLA-DPA1, CD74, MS4A1, HLA-DRB1, HLA-DRA
## HLA-DRB5, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, CD37, HVCN1, FCRLA, IRF8
## PLAC8, BLNK, MALAT1, SMIM14, PLD4, P2RX5, IGLL5, LAT2, SWAP70, FCGR2B
## PC_ 4
## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1, CD74, HIST1H2AC, HLA-DPB1, PF4, SDPR
## TCL1A, HLA-DRB1, HLA-DPA1, HLA-DQA2, PPBP, HLA-DRA, LINC00926, GNG11, SPARC, HLA-DRB5
## GP9, AP001189.4, CA2, PTCRA, CD9, NRGN, RGS18, CLU, TUBB1, GZMB
## Negative: VIM, IL7R, S100A6, IL32, S100A8, S100A4, GIMAP7, S100A10, S100A9, MAL
## AQP3, CD2, CD14, FYB, LGALS2, GIMAP4, ANXA1, CD27, FCN1, RBP7
## LYZ, S100A11, GIMAP5, MS4A6A, S100A12, FOLR3, TRABD2A, AIF1, IL8, IFI6
## PC_ 5
## Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY, CCL4, CST7, PRF1, GZMA, SPON2
## GZMH, S100A9, LGALS2, CCL3, CTSW, XCL2, CD14, CLIC3, S100A12, RBP7
## CCL5, MS4A6A, GSTP1, FOLR3, IGFBP7, TYROBP, TTC38, AKR1C3, XCL1, HOPX
## Negative: LTB, IL7R, CKB, VIM, MS4A7, AQP3, CYTIP, RP11-290F20.3, SIGLEC10, HMOX1
## LILRB2, PTGES3, MAL, CD27, HN1, CD2, GDI2, CORO1B, ANXA5, TUBA1B
## FAM110A, ATP1A1, TRADD, PPA1, CCDC109B, ABRACL, CTD-2006K23.1, WARS, VMO1, FYB
# We will check first 15 PCA to check which PCA explain most of variation
DimHeatmap(pbmc_pca, dims = 1:15, cells = 500, balanced = TRUE)
# Determine the ‘dimensionality’ of the dataset
ElbowPlot(pbmc_pca)
# Cluster the cells
# ElbowPlot showed first 8 dimension reduced the data massively, so we don't need rest
pbmc_neigh <- FindNeighbors(pbmc_pca, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
pbmc_clus <- FindClusters(pbmc_neigh, resolution = 0.5) # resolution is theh most important feature. It can categorized to more cluster
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 95927
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8728
## Number of communities: 9
## Elapsed time: 0 seconds
# let's try with resolution =v3
#pbmc_clus <- FindClusters(pbmc_neigh, resolution = 3) # it increases number of cluster from 8 to 21
# Run non-linear dimensional reduction (UMAP / t-SNE)
pbmc_umap <- RunUMAP(pbmc_clus, dims = 1:10)
## 23:44:39 UMAP embedding parameters a = 0.9922 b = 1.112
## 23:44:39 Read 2638 rows and found 10 numeric columns
## 23:44:39 Using Annoy for neighbor search, n_neighbors = 30
## 23:44:39 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 23:44:40 Writing NN index file to temp file C:\Users\ahuqu\AppData\Local\Temp\Rtmp4UuG7u\file6b9c207468cb
## 23:44:40 Searching Annoy index using 1 thread, search_k = 3000
## 23:44:40 Annoy recall = 100%
## 23:44:41 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 23:44:41 Initializing from normalized Laplacian + noise (using RSpectra)
## 23:44:41 Commencing optimization for 500 epochs, with 105140 positive edges
## 23:44:41 Using rng type: pcg
## 23:44:50 Optimization finished
DimPlot(pbmc_umap, reduction = "umap")
# Finding differentially expressed features (cluster biomarkers)
# find all markers of cluster 2
cluster2.markers <- FindMarkers(pbmc_umap, ident.1 = 2)
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session
head(cluster2.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## IL32 2.894376e-90 1.3070772 0.947 0.465 3.969347e-86
## LTB 1.061535e-86 1.3312674 0.981 0.643 1.455789e-82
## CD3D 8.794641e-71 1.0597620 0.922 0.432 1.206097e-66
## IL7R 3.516098e-68 1.4377848 0.750 0.326 4.821977e-64
## LDHB 1.641513e-67 0.9911924 0.954 0.614 2.251171e-63
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc_umap, ident.1 = 5, ident.2 = c(0, 3))
head(cluster5.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## FCGR3A 8.246578e-205 6.794969 0.975 0.040 1.130936e-200
## IFITM3 1.677613e-195 6.192558 0.975 0.049 2.300678e-191
## CFD 2.401156e-193 6.015172 0.938 0.038 3.292945e-189
## CD68 2.900384e-191 5.530330 0.926 0.035 3.977587e-187
## RP11-290F20.3 2.513244e-186 6.297999 0.840 0.017 3.446663e-182
VlnPlot(pbmc_umap, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc_umap, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
## Warning: The `slot` argument of `VlnPlot()` is deprecated as of Seurat 5.0.0.
## ℹ Please use the `layer` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# FeaturePlot. If now marker genes for your desired cell type list them. It will show in which cluster our marker genes is present
FeaturePlot(pbmc_umap, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
saveRDS(pbmc, file = "I:/OxfordBioDiscovery2025/2026/Single cell RNAseq/2. Single Cell RNAseq/pbmc3k_final.rds")
# Assigning cluster names according to the cluster marker genes
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T",
"FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc_umap)
pbmc_umap <- RenameIdents(pbmc_umap, new.cluster.ids)
DimPlot(pbmc_umap, reduction = "umap", label = TRUE, pt.size = 0.5) +
NoLegend()
