SS_Robj_CellLines_PBMC_Merged_without_Correction_L1_Merged_first_HPC_PC-1:50-3
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.379441 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", 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: TRAV38-2DV8, CD7, KRT1, MALAT1, LTB, GIMAP7, TRBV20-1, RPS27, RPS4Y1, CD52
## CST7, XCL1, GAS5, RPLP1, IL7R, TTC29, KIR3DL1, GIMAP4, GIMAP5, CD5
## GZMM, TCF7, RPL13, RPS12, KLF2, NKG7, PNRC1, GZMA, RPS15, RPS28
## Negative: CCL17, CD74, TNFRSF4, PPBP, IL2RA, LGALS1, CA2, MIR155HG, SEC11C, FABP5
## CD70, SYT4, C12orf75, EGFL6, MIIP, VIM, BATF3, KRT7, LGALS3, EEF1A2
## MT2A, LTA, MTHFD2, YBX3, GAPDH, FN1, CCL5, RDH10, NPM1, HDGFL3
## PC_ 2
## Positive: CCL17, CA2, TNFRSF4, SYT4, EGFL6, C12orf75, CA10, CCL5, KRT1, TRAV38-2DV8
## IGHE, PRG4, PLPP1, CFI, THY1, GAS5, TIGIT, CYBA, STC1, LTB
## HACD1, RANBP17, MAP4K4, SPINT2, ONECUT2, BACE2, SEC11C, LAYN, IL13, RXFP1
## Negative: PPBP, CD74, MT2A, PAGE5, LMNA, TENM3, CD70, STAT1, TRAV17, RPL22L1
## GSTP1, LGALS3, RBPMS, LGALS1, SPOCK1, PPP2R2B, CCDC50, MACROD2, TRAV9-2, IQCG
## GAPDH, FTL, CTAG2, FABP5, ANXA2, NDUFV2, PIM2, SLC7A11-AS1, TNFSF10, TMSB4X
## PC_ 3
## Positive: RPS4Y1, MALAT1, IL7R, PNRC1, BTG1, LINC00861, B2M, GIMAP7, SELL, TCF7
## CCL17, GIMAP5, SARAF, CCR7, PIK3IP1, GIMAP4, RPS27, PCED1B-AS1, PABPC1, ZFP36
## FYB1, SESN3, FTH1, TRBC1, TRBC2, YPEL3, KLF2, GIMAP1, BCL3, ITM2B
## Negative: XCL1, TRAV38-2DV8, KIR3DL1, CD7, XCL2, KIR2DL3, MT1G, CST7, KLRC1, ACTB
## HIST1H4C, KIR2DL4, LTB, KRT81, ESYT2, C1QBP, KRT86, MATK, IFITM2, MYO1E
## ID3, EPCAM, KIR3DL2, TRGV2, NKG7, EIF5A, CXCR3, PFN1, CHCHD2, HSP90AA1
## PC_ 4
## Positive: CD7, XCL1, XCL2, RPS4Y1, KIR3DL1, IL7R, MALAT1, KIR2DL3, IFITM2, S100A4
## MT1G, KLRC1, KIR2DL4, LSP1, TMSB4X, KRT81, IFITM1, KRT86, KLRK1, SARAF
## PRKCH, MYO1E, LINC00861, TCF7, FTH1, BTG1, MATK, CAPG, PTPN6, ESYT2
## Negative: KRT1, TTC29, TRBV20-1, PPBP, RPL13, RPLP1, GAS5, GZMA, RPS4X, LINC02752
## EEF1A1, CEBPD, WFDC1, AC069410.1, SP5, TBX4, RPS12, IL4, CD74, PLCB1
## EEF2, TNS4, RPS15, RPLP0, RPS28, IFNGR1, EGLN3, HSPB1, RPS2, RPS18
## PC_ 5
## Positive: CCL17, PPBP, MT2A, LTA, CA2, CD74, MIR155HG, CD7, XCL1, CCL5
## CA10, MGST3, STC1, FCER2, XCL2, RXFP1, IQCG, KIR2DL3, CFI, RPL22L1
## RANBP17, KIR3DL1, AL590550.1, RPS4Y1, THY1, PAGE5, RAP1A, STAT1, IGHE, RYR2
## Negative: EEF1A2, TNFRSF4, IL2RA, WFDC1, PHLDA2, S100A4, FN1, MIIP, S100A11, HIST1H1C
## TRAV4, S100A6, PXYLP1, LGALS1, DUSP4, RDH10, GPAT3, TIGIT, TNFRSF18, CDKN1A
## AL136456.1, HOXC9, GATA3, CEP135, HIST1H2BK, HACD1, TP73, TMEM163, MAP1B, LIME1
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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9549
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9468
## Number of communities: 24
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9396
## Number of communities: 26
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9328
## Number of communities: 29
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9271
## Number of communities: 31
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9220
## Number of communities: 35
## 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: 1663835
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9172
## 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
## 3416 3360 3313 2902 2869 2666 2561 2363 2240 2192 2121 2082 1957 1838 1729 1541
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
## 1229 919 833 589 581 507 502 443 439 421 334 242 217 128 117 106
## 32 33 34
## 91 86 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## 10:52:47 Read 46976 rows## 10:52:47 Processing block 1 of 1## 10:52:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 10:52:48 Initializing by weighted average of neighbor coordinates using 1 thread
## 10:52:48 Commencing optimization for 67 epochs, with 939520 positive edges
## 10:52:52 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.3896441459656DimPlot(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: 4 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: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## 0 1 2 3 4 5 6 7 8 9 10 11
## ASDC 47 0 0 8 0 0 26 7 0 0 0 0
## B intermediate 4 0 0 0 3 0 0 0 0 0 0 5
## B memory 1 0 0 0 0 1 0 0 0 0 0 0
## CD4 CTL 0 0 0 0 0 0 0 0 0 0 6 0
## CD4 Naive 0 0 0 0 0 667 0 0 0 0 5 0
## CD4 Proliferating 5295 2901 2412 3950 3667 0 3876 2886 1339 1173 0 353
## CD4 TCM 442 291 2813 175 63 3906 326 65 7 6 31 3
## CD4 TEM 0 0 11 0 0 27 0 0 0 0 15 0
## CD8 Naive 1 0 0 41 1 41 147 15 0 0 310 0
## CD8 Proliferating 0 0 0 0 0 0 0 0 0 0 0 0
## CD8 TCM 0 14 434 0 0 39 0 0 0 0 199 0
## CD8 TEM 0 1 0 0 0 1 0 0 0 0 202 0
## cDC2 484 0 0 516 832 0 38 206 94 6 0 73
## dnT 1 1 3 0 0 23 0 0 0 0 6 0
## gdT 0 0 0 0 0 0 0 0 0 0 13 0
## HSPC 42 0 0 424 181 3 0 645 396 1 0 45
## ILC 0 0 0 0 0 1 0 0 0 0 0 0
## MAIT 0 0 0 0 0 28 0 0 0 0 30 0
## NK 0 0 0 0 0 0 0 0 0 0 91 0
## NK Proliferating 3 2725 17 29 207 0 37 11 1 5 0 9
## Platelet 0 0 0 0 0 1 0 0 0 0 0 0
## Treg 1 0 1 1 0 146 0 0 0 0 0 0
##
## 12 13 14 15 16 17 18
## ASDC 0 0 0 0 0 1 0
## B intermediate 0 34 0 2 0 0 2
## B memory 0 1 0 1 0 0 0
## CD4 CTL 0 0 0 0 0 0 0
## CD4 Naive 70 0 0 5 0 0 0
## CD4 Proliferating 0 106 214 57 39 73 12
## CD4 TCM 312 39 1 85 72 4 9
## CD4 TEM 6 0 0 0 0 0 0
## CD8 Naive 23 0 1 3 0 0 0
## CD8 Proliferating 0 1 0 0 0 0 0
## CD8 TCM 14 0 0 1 17 0 0
## CD8 TEM 4 0 0 0 0 0 0
## cDC2 0 8 0 18 0 7 0
## dnT 0 7 0 13 0 0 0
## gdT 0 0 0 0 0 0 0
## HSPC 0 2 0 0 0 0 3
## ILC 0 0 0 0 0 0 0
## MAIT 0 0 0 0 0 0 0
## NK 0 0 0 1 0 0 0
## NK Proliferating 0 22 0 0 0 1 3
## Platelet 0 1 0 0 0 0 0
## Treg 13 12 0 7 0 0 13