SS_Robj_CellLines_PBMC_Merged_without_Correction_L1_Merged_first_HPC_PC-1:50-4
1. load libraries
2. Load Seurat Object
#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
load("All_T_cells_Merged_filtered_Mono_using_clusters.Robj")
All_samples_Merged <- filtered_data
All_samples_Merged## An object of class Seurat
## 62626 features across 46976 samples within 6 assays
## Active assay: SCT (25902 features, 3000 variable features)
## 3 layers present: counts, data, scale.data
## 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
## 4 dimensional reductions calculated: pca, umap, integrated_dr, ref.umap3. QC
Idents(object = All_samples_Merged) <- "cell_line"
VlnPlot(All_samples_Merged, features = c("nFeature_RNA",
"nCount_RNA",
"percent.mito"),
ncol = 3)
FeatureScatter(All_samples_Merged,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA") +
geom_smooth(method = 'lm')## `geom_smooth()` using formula = 'y ~ x'
##FeatureScatter is typically used to visualize feature-feature relationships ##for anything calculated by the object, ##i.e. columns in object metadata, PC scores etc.
FeatureScatter(All_samples_Merged,
feature1 = "nCount_RNA",
feature2 = "percent.mito")+
geom_smooth(method = 'lm')## `geom_smooth()` using formula = 'y ~ x'
FeatureScatter(All_samples_Merged,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA")+
geom_smooth(method = 'lm')## `geom_smooth()` using formula = 'y ~ x'
4. Normalize data
## Running SCTransform on assay: RNA## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 25901 by 46976## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells## Found 484 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 25901 genes## Computing corrected count matrix for 25901 genes## Calculating gene attributes## Wall clock passed: Time difference of 2.400643 mins## Determine variable features## Centering data matrix## Place corrected count matrix in counts slot## Warning: Different cells and/or features from existing assay SCT## Set default assay to SCT5. Perform PCA
Variables_genes <- All_samples_Merged@assays$SCT@var.features
# Exclude genes starting with "HLA-" or "Xist"
Variables_genes_after_exclusion <- Variables_genes[!grepl("^HLA-|^XIST|^TRBV|^TRAV", Variables_genes)]
# These are now standard steps in the Seurat workflow for visualization and clustering
All_samples_Merged <- RunPCA(All_samples_Merged,
features = Variables_genes_after_exclusion,
do.print = TRUE,
pcs.print = 1:5,
genes.print = 15,
npcs = 50)## PC_ 1
## Positive: CCL17, CD74, TNFRSF4, PPBP, IL2RA, LGALS1, CA2, MIR155HG, SEC11C, FABP5
## CD70, SYT4, C12orf75, EGFL6, MIIP, KRT7, BATF3, VIM, LGALS3, EEF1A2
## LTA, MT2A, MTHFD2, YBX3, GAPDH, FN1, CCL5, NPM1, RDH10, HDGFL3
## Negative: MALAT1, CD7, KRT1, LTB, GIMAP7, RPS4Y1, RPS27, CD52, CST7, IL7R
## GAS5, XCL1, RPLP1, GIMAP4, GIMAP5, TTC29, TCF7, KIR3DL1, CD5, PNRC1
## GZMM, RPL13, KLF2, RPS12, NKG7, SELL, RPL41, RPS15, RPS28, LEF1
## PC_ 2
## Positive: CCL17, CA2, TNFRSF4, SYT4, EGFL6, CCL5, CA10, C12orf75, KRT1, IGHE
## PRG4, GAS5, CFI, THY1, TIGIT, PLPP1, CYBA, STC1, RANBP17, MAP4K4
## HACD1, LTB, ONECUT2, BACE2, SEC11C, SPINT2, LAYN, IL13, HS3ST1, RXFP1
## Negative: PPBP, CD74, MT2A, PAGE5, LMNA, TENM3, CD70, STAT1, RPL22L1, GSTP1
## RBPMS, LGALS3, LGALS1, GAPDH, CCDC50, SPOCK1, PPP2R2B, MACROD2, FTL, IQCG
## FABP5, CTAG2, ANXA2, TMSB4X, SLC7A11-AS1, NDUFV2, STMN1, BASP1, CDKN2A, PIM2
## PC_ 3
## Positive: RPS4Y1, MALAT1, IL7R, B2M, PNRC1, BTG1, LINC00861, SELL, GIMAP7, TCF7
## GIMAP5, SARAF, CCR7, PIK3IP1, GIMAP4, PCED1B-AS1, RPS27, FYB1, PABPC1, ZFP36
## SESN3, TRBC1, FTH1, CCL17, TRBC2, YPEL3, BCL3, KLF2, GIMAP1, ANXA1
## Negative: XCL1, CD7, KIR3DL1, XCL2, KIR2DL3, MT1G, CST7, KLRC1, KIR2DL4, HIST1H4C
## ACTB, LTB, KRT81, ESYT2, KRT86, C1QBP, MATK, MYO1E, IFITM2, EPCAM
## ID3, TRGV2, KIR3DL2, CXCR3, NKG7, EIF5A, HIST1H1B, PTPN6, CHCHD2, GZMM
## PC_ 4
## Positive: KRT1, TTC29, PPBP, RPL13, RPLP1, GAS5, GZMA, RPS4X, LINC02752, EEF1A1
## CEBPD, WFDC1, AC069410.1, SP5, TBX4, IL4, RPS12, TNS4, RPS15, PLCB1
## EEF2, RPLP0, RPS28, IFNGR1, RPS2, CD74, HSPB1, RPS18, RPS8, EGLN3
## Negative: CD7, XCL1, RPS4Y1, XCL2, IL7R, MALAT1, S100A4, KIR3DL1, KIR2DL3, IFITM2
## LSP1, KLRC1, MT1G, TMSB4X, KIR2DL4, SARAF, LINC00861, TCF7, IFITM1, BTG1
## KRT81, PRKCH, FTH1, KLRK1, KRT86, MYO1E, RGCC, PNRC1, ZFP36, CAPG
## PC_ 5
## Positive: EEF1A2, TNFRSF4, IL2RA, WFDC1, PHLDA2, S100A4, FN1, MIIP, S100A11, HIST1H1C
## S100A6, PXYLP1, LGALS1, DUSP4, RDH10, GPAT3, TIGIT, TNFRSF18, CDKN1A, HOXC9
## AL136456.1, GATA3, CEP135, KRT1, HIST1H2BK, HACD1, TP73, TMEM163, MAP1B, LIME1
## Negative: CCL17, PPBP, MT2A, LTA, CD74, CA2, CD7, MIR155HG, XCL1, CCL5
## CA10, MGST3, STC1, XCL2, FCER2, KIR2DL3, RXFP1, IQCG, CFI, RANBP17
## RPL22L1, KIR3DL1, RPS4Y1, AL590550.1, STAT1, MALAT1, PAGE5, THY1, RAP1A, RYR2
Perform PCA TEST
library(ggplot2)
library(RColorBrewer)
# Assuming you have 10 different cell lines, generating a color palette with 10 colors
cell_line_colors <- brewer.pal(10, "Set3")
# Assuming All_samples_Merged$cell_line is a factor or character vector containing cell line names
data <- as.data.frame(table(All_samples_Merged$cell_line))
colnames(data) <- c("cell_line", "nUMI") # Change column name to nUMI
ncells <- ggplot(data, aes(x = cell_line, y = nUMI, fill = cell_line)) +
geom_col() +
theme_classic() +
geom_text(aes(label = nUMI),
position = position_dodge(width = 0.9),
vjust = -0.25) +
scale_fill_manual(values = cell_line_colors) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5)) + # Adjust the title position
ggtitle("Filtered cells per sample") +
xlab("Cell lines") + # Adjust x-axis label
ylab("Frequency") # Adjust y-axis label
print(ncells)
# TEST-1
# given that the output of RunPCA is "pca"
# replace "so" by the name of your seurat object
pct <- All_samples_Merged[["pca"]]@stdev / sum(All_samples_Merged[["pca"]]@stdev) * 100
cumu <- cumsum(pct) # Calculate cumulative percents for each PC
# Determine the difference between variation of PC and subsequent PC
co2 <- sort(which((pct[-length(pct)] - pct[-1]) > 0.1), decreasing = T)[1] + 1
# last point where change of % of variation is more than 0.1%. -> co2
co2## [1] 15# TEST-2
# get significant PCs
stdv <- All_samples_Merged[["pca"]]@stdev
sum.stdv <- sum(All_samples_Merged[["pca"]]@stdev)
percent.stdv <- (stdv / sum.stdv) * 100
cumulative <- cumsum(percent.stdv)
co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
co2 <- sort(which((percent.stdv[1:length(percent.stdv) - 1] -
percent.stdv[2:length(percent.stdv)]) > 0.1),
decreasing = T)[1] + 1
min.pc <- min(co1, co2)
min.pc## [1] 15# Create a dataframe with values
plot_df <- data.frame(pct = percent.stdv,
cumu = cumulative,
rank = 1:length(percent.stdv))
# Elbow plot to visualize
ggplot(plot_df, aes(cumulative, percent.stdv, label = rank, color = rank > min.pc)) +
geom_text() +
geom_vline(xintercept = 90, color = "grey") +
geom_hline(yintercept = min(percent.stdv[percent.stdv > 5]), color = "grey") +
theme_bw()
6. Clustering
All_samples_Merged <- FindNeighbors(All_samples_Merged,
dims = 1:50,
verbose = FALSE)
# understanding resolution
All_samples_Merged <- FindClusters(All_samples_Merged,
resolution = c(0.4,0.5, 0.6, 0.7,0.8, 0.9, 1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9546
## Number of communities: 19
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9464
## Number of communities: 23
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9390
## Number of communities: 26
## Elapsed time: 7 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9321
## Number of communities: 28
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9262
## Number of communities: 34
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9211
## Number of communities: 32
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 46976
## Number of edges: 1659128
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9161
## Number of communities: 35
## Elapsed time: 7 seconds# non-linear dimensionality reduction --------------
All_samples_Merged <- RunUMAP(All_samples_Merged,
dims = 1:50,
verbose = FALSE)## 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# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(All_samples_Merged,group.by = "cell_line",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
DimPlot(All_samples_Merged,
group.by = "SCT_snn_res.0.7",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
cluster_table <- table(Idents(All_samples_Merged))
barplot(cluster_table, main = "Number of Cells in Each Cluster",
xlab = "Cluster",
ylab = "Number of Cells",
col = rainbow(length(cluster_table)))
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 4051 3715 3399 3226 2734 2698 2393 2204 2193 2158 2147 1919 1889 1837 1736 1203
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 1035 927 889 657 645 506 501 444 384 344 255 220 130 118 107 93
## 32 33 34
## 92 85 427. Azimuth Annotation
## Warning: The following packages are already installed and will not be
## reinstalled: pbmcref# The RunAzimuth function can take a Seurat object as input
All_samples_Merged <- RunAzimuth(All_samples_Merged, reference = "pbmcref")## Warning: Overwriting miscellanous data for model## Warning: Adding a dimensional reduction (refUMAP) without the associated assay
## being present
## Warning: Adding a dimensional reduction (refUMAP) without the associated assay
## being present## detected inputs from HUMAN with id type Gene.name## reference rownames detected HUMAN with id type Gene.name## Normalizing query using reference SCT model## Warning: 113 features of the features specified were not present in both the reference query assays.
## Continuing with remaining 4887 features.## Projecting cell embeddings## Finding query neighbors## Finding neighborhoods## Finding anchors## Found 4803 anchors## Finding integration vectors## Finding integration vector weights## Predicting cell labels
## Predicting cell labels## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')## Predicting cell labels## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')##
## Integrating dataset 2 with reference dataset## Finding integration vectors## Integrating data## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from integrated_dr_ to integrateddr_## Computing nearest neighbors## Running UMAP projection## Warning in RunUMAP.default(object = neighborlist, reduction.model =
## reduction.model, : Number of neighbors between query and reference is not equal
## to the number of neighbors within reference## 11:59:27 Read 46976 rows## 11:59:27 Processing block 1 of 1## 11:59:28 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 11:59:28 Initializing by weighted average of neighbor coordinates using 1 thread
## 11:59:28 Commencing optimization for 67 epochs, with 939520 positive edges
## 11:59:32 Finished## Warning: No assay specified, setting assay as RNA by default.## Projecting reference PCA onto query
## Finding integration vector weights
## Projecting back the query cells into original PCA space
## Finding integration vector weights
## Computing scores:
## Finding neighbors of original query cells
## Finding neighbors of transformed query cells
## Computing query SNN
## Determining bandwidth and computing transition probabilities
## Total elapsed time: 21.1084017753601DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
8. Cell type annotation using ProjectTils
#Load reference atlas and query data
ref <- readRDS(file = "../../8-Cell_Lines_Test/CD4T_human_ref_v1.rds")
#Run Projection algorithm
query.projected <- Run.ProjecTILs(All_samples_Merged, ref = ref)## | | | 0%[1] "Using assay SCT for query"## Pre-filtering cells with scGate...##
## ### Detected a total of 26819 pure 'Target' cells (57.09% of total)## [1] "20157 out of 46976 ( 43% ) non-pure cells removed. Use filter.cells=FALSE to avoid pre-filtering"
## [1] "Aligning query to reference map for batch-correction..."## Warning: Layer counts isn't present in the assay object[[assay]]; returning
## NULL## Preparing PCA embeddings for objects...## Warning: Number of dimensions changing from 50 to 20##
## Projecting corrected query onto Reference PCA space
##
## Projecting corrected query onto Reference UMAP space## Warning: Not all features provided are in this Assay object, removing the
## following feature(s): CD177, TASL, H1-4, H2AZ1, ELAPOR1, CCL3L3, POLR1F, AHSP,
## H1-2, H1-0, WARS1, H1-3, H2BC11, GPX1, H2AC6, IRAG2, H1-10, H3C10, IL22, ECEL1,
## H4C3## | |======================================================================| 100%## Creating slots functional.cluster and functional.cluster.conf in query object
#Plot the predicted composition of the query in terms of reference T cell subtypes
plot.statepred.composition(ref, query.projected, metric = "Percent")
## | | | 0%[1] "Using assay SCT for query"## Pre-filtering cells with scGate...
##
## ### Detected a total of 26819 pure 'Target' cells (57.09% of total)## [1] "20157 out of 46976 ( 43% ) non-pure cells removed. Use filter.cells=FALSE to avoid pre-filtering"
## [1] "Aligning query to reference map for batch-correction..."## Warning: Layer counts isn't present in the assay object[[assay]]; returning
## NULL## Preparing PCA embeddings for objects...## Warning: Number of dimensions changing from 50 to 20##
## Projecting corrected query onto Reference PCA space
##
## Projecting corrected query onto Reference UMAP space## Warning: Not all features provided are in this Assay object, removing the
## following feature(s): CD177, TASL, H1-4, H2AZ1, ELAPOR1, CCL3L3, POLR1F, AHSP,
## H1-2, H1-0, WARS1, H1-3, H2BC11, GPX1, H2AC6, IRAG2, H1-10, H3C10, IL22, ECEL1,
## H4C3## | |======================================================================| 100%## Creating slots functional.cluster and functional.cluster.conf in query objectDimPlot(All_samples_Merged, group.by = "functional.cluster",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
clusTree
## Loading required package: ggraph##
## Attaching package: 'ggraph'## The following object is masked from 'package:sp':
##
## geometry
Azimuth Visualization
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## 0 1 2 3 4 5 6 7 8 9 10 11
## ASDC 45 0 0 0 0 26 7 0 0 8 0 0
## B intermediate 4 0 0 0 3 0 0 0 0 0 0 0
## B memory 1 0 0 1 0 0 0 0 0 0 0 0
## CD4 CTL 0 0 0 0 0 0 0 0 0 0 0 6
## CD4 Naive 0 0 0 741 0 0 0 0 0 0 0 1
## CD4 Proliferating 5268 2901 2412 0 3669 3867 2887 2857 1341 1080 1180 0
## CD4 TCM 445 291 2810 4198 61 326 66 4 7 171 6 50
## CD4 TEM 0 0 11 15 0 0 0 0 0 0 0 33
## CD8 Naive 1 0 0 48 1 147 15 0 0 41 0 326
## CD8 Proliferating 0 0 0 0 0 0 0 0 0 0 0 0
## CD8 TCM 0 14 434 34 0 0 0 0 0 0 0 218
## CD8 TEM 0 1 0 1 0 0 0 0 0 0 0 206
## cDC2 471 0 0 0 827 38 206 71 92 449 6 0
## dnT 1 1 3 23 0 0 0 0 0 0 0 5
## gdT 0 0 0 0 0 0 0 0 0 0 0 13
## HSPC 40 0 0 3 180 0 645 424 396 0 1 0
## ILC 0 0 0 1 0 0 0 0 0 0 0 0
## MAIT 0 0 0 5 0 0 0 0 0 0 0 53
## NK 0 0 0 0 0 0 0 0 0 0 0 91
## NK Proliferating 3 2725 17 0 207 37 12 6 0 21 5 0
## Platelet 0 0 0 1 0 0 0 0 0 0 0 0
## Treg 1 0 1 158 0 0 0 0 0 1 0 0
##
## 12 13 14 15 16 17 18
## ASDC 0 0 0 0 2 1 0
## B intermediate 5 36 0 0 0 0 2
## B memory 0 2 0 0 0 0 0
## CD4 CTL 0 0 0 0 0 0 0
## CD4 Naive 0 5 0 0 0 0 0
## CD4 Proliferating 355 123 217 39 73 72 12
## CD4 TCM 2 121 1 74 4 4 9
## CD4 TEM 0 0 0 0 0 0 0
## CD8 Naive 0 3 1 0 0 0 0
## CD8 Proliferating 0 1 0 0 0 0 0
## CD8 TCM 0 1 0 17 0 0 0
## CD8 TEM 0 0 0 0 0 0 0
## cDC2 74 15 0 0 26 7 0
## dnT 0 21 0 0 0 0 0
## gdT 0 0 0 0 0 0 0
## HSPC 45 3 0 0 2 0 3
## ILC 0 0 0 0 0 0 0
## MAIT 0 0 0 0 0 0 0
## NK 0 1 0 0 0 0 0
## NK Proliferating 9 24 0 0 0 1 3
## Platelet 0 1 0 0 0 0 0
## Treg 0 20 0 0 0 0 13