#Load Seurat Object L4
load("../Documents/1-SS-STeps/4-Analysis_and_Robj_Marie/analyse juillet 2023/ObjetsR/L4_B.Robj")
L4 <- L4_B
L4
An object of class Seurat
36629 features across 6150 samples within 2 assays
Active assay: RNA (36601 features, 0 variable features)
2 layers present: counts, data
1 other assay present: ADT
# Set identity classes to an existing column in meta data
Idents(object = L4) <- "cell_line"
L4[["percent.rb"]] <- PercentageFeatureSet(L4, pattern = "^RP[SL]")
VlnPlot(L4, features = c("nFeature_RNA",
"nCount_RNA",
"percent.mito",
"percent.rb"),
ncol = 4, pt.size = 0.1) &
theme(plot.title = element_text(size=10))
FeatureScatter(L4, feature1 = "percent.mito",
feature2 = "percent.rb")
VlnPlot(L4, features = c("nFeature_RNA",
"nCount_RNA",
"percent.mito"),
ncol = 3)
FeatureScatter(L4,
feature1 = "percent.mito",
feature2 = "percent.rb") +
geom_smooth(method = 'lm')
FeatureScatter(L4,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA") +
geom_smooth(method = 'lm')
##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(L4,
feature1 = "nCount_RNA",
feature2 = "percent.mito")+
geom_smooth(method = 'lm')
FeatureScatter(L4,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA")+
geom_smooth(method = 'lm')
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 19823 by 6150
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 5000 cells
Found 142 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 19823 genes
Computing corrected count matrix for 19823 genes
Calculating gene attributes
Wall clock passed: Time difference of 24.85647 secs
Determine variable features
|
| | 0%
|
|============================== | 25%
|
|============================================================= | 50%
|
|============================================================================================ | 75%
|
|==========================================================================================================================| 100%
Place corrected count matrix in counts slot
Set default assay to SCT
Warning: The following features are not present in the object: MLF1IP, not searching for symbol synonymsWarning: The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms
# Apply SCTransform
L4 <- SCTransform(L4, vars.to.regress = c("percent.rb","percent.mito", "CC.Difference"),
do.scale=TRUE,
do.center=TRUE,
verbose = TRUE)
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 19823 by 6150
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 5000 cells
Found 142 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 19823 genes
Computing corrected count matrix for 19823 genes
Calculating gene attributes
Wall clock passed: Time difference of 21.62516 secs
Determine variable features
Regressing out percent.rb, percent.mito, CC.Difference
|
| | 0%
|
|= | 0%
|
|= | 1%
|
|== | 1%
|
|== | 2%
|
|=== | 2%
|
|=== | 3%
|
|==== | 3%
|
|==== | 4%
|
|===== | 4%
|
|====== | 5%
|
|======= | 5%
|
|======= | 6%
|
|======== | 6%
|
|======== | 7%
|
|========= | 7%
|
|========= | 8%
|
|========== | 8%
|
|========== | 9%
|
|=========== | 9%
|
|============ | 9%
|
|============ | 10%
|
|============= | 10%
|
|============= | 11%
|
|============== | 11%
|
|============== | 12%
|
|=============== | 12%
|
|=============== | 13%
|
|================ | 13%
|
|================ | 14%
|
|================= | 14%
|
|================== | 14%
|
|================== | 15%
|
|=================== | 15%
|
|=================== | 16%
|
|==================== | 16%
|
|==================== | 17%
|
|===================== | 17%
|
|===================== | 18%
|
|====================== | 18%
|
|======================= | 18%
|
|======================= | 19%
|
|======================== | 19%
|
|======================== | 20%
|
|========================= | 20%
|
|========================= | 21%
|
|========================== | 21%
|
|========================== | 22%
|
|=========================== | 22%
|
|=========================== | 23%
|
|============================ | 23%
|
|============================= | 23%
|
|============================= | 24%
|
|============================== | 24%
|
|============================== | 25%
|
|=============================== | 25%
|
|=============================== | 26%
|
|================================ | 26%
|
|================================ | 27%
|
|================================= | 27%
|
|================================== | 27%
|
|================================== | 28%
|
|=================================== | 28%
|
|=================================== | 29%
|
|==================================== | 29%
|
|==================================== | 30%
|
|===================================== | 30%
|
|===================================== | 31%
|
|====================================== | 31%
|
|====================================== | 32%
|
|======================================= | 32%
|
|======================================== | 32%
|
|======================================== | 33%
|
|========================================= | 33%
|
|========================================= | 34%
|
|========================================== | 34%
|
|========================================== | 35%
|
|=========================================== | 35%
|
|=========================================== | 36%
|
|============================================ | 36%
|
|============================================= | 36%
|
|============================================= | 37%
|
|============================================== | 37%
|
|============================================== | 38%
|
|=============================================== | 38%
|
|=============================================== | 39%
|
|================================================ | 39%
|
|================================================ | 40%
|
|================================================= | 40%
|
|================================================= | 41%
|
|================================================== | 41%
|
|=================================================== | 41%
|
|=================================================== | 42%
|
|==================================================== | 42%
|
|==================================================== | 43%
|
|===================================================== | 43%
|
|===================================================== | 44%
|
|====================================================== | 44%
|
|====================================================== | 45%
|
|======================================================= | 45%
|
|======================================================== | 46%
|
|========================================================= | 46%
|
|========================================================= | 47%
|
|========================================================== | 47%
|
|========================================================== | 48%
|
|=========================================================== | 48%
|
|=========================================================== | 49%
|
|============================================================ | 49%
|
|============================================================ | 50%
|
|============================================================= | 50%
|
|============================================================== | 50%
|
|============================================================== | 51%
|
|=============================================================== | 51%
|
|=============================================================== | 52%
|
|================================================================ | 52%
|
|================================================================ | 53%
|
|================================================================= | 53%
|
|================================================================= | 54%
|
|================================================================== | 54%
|
|================================================================== | 55%
|
|=================================================================== | 55%
|
|==================================================================== | 55%
|
|==================================================================== | 56%
|
|===================================================================== | 56%
|
|===================================================================== | 57%
|
|====================================================================== | 57%
|
|====================================================================== | 58%
|
|======================================================================= | 58%
|
|======================================================================= | 59%
|
|======================================================================== | 59%
|
|========================================================================= | 59%
|
|========================================================================= | 60%
|
|========================================================================== | 60%
|
|========================================================================== | 61%
|
|=========================================================================== | 61%
|
|=========================================================================== | 62%
|
|============================================================================ | 62%
|
|============================================================================ | 63%
|
|============================================================================= | 63%
|
|============================================================================= | 64%
|
|============================================================================== | 64%
|
|=============================================================================== | 64%
|
|=============================================================================== | 65%
|
|================================================================================ | 65%
|
|================================================================================ | 66%
|
|================================================================================= | 66%
|
|================================================================================= | 67%
|
|================================================================================== | 67%
|
|================================================================================== | 68%
|
|=================================================================================== | 68%
|
|==================================================================================== | 68%
|
|==================================================================================== | 69%
|
|===================================================================================== | 69%
|
|===================================================================================== | 70%
|
|====================================================================================== | 70%
|
|====================================================================================== | 71%
|
|======================================================================================= | 71%
|
|======================================================================================= | 72%
|
|======================================================================================== | 72%
|
|======================================================================================== | 73%
|
|========================================================================================= | 73%
|
|========================================================================================== | 73%
|
|========================================================================================== | 74%
|
|=========================================================================================== | 74%
|
|=========================================================================================== | 75%
|
|============================================================================================ | 75%
|
|============================================================================================ | 76%
|
|============================================================================================= | 76%
|
|============================================================================================= | 77%
|
|============================================================================================== | 77%
|
|=============================================================================================== | 77%
|
|=============================================================================================== | 78%
|
|================================================================================================ | 78%
|
|================================================================================================ | 79%
|
|================================================================================================= | 79%
|
|================================================================================================= | 80%
|
|================================================================================================== | 80%
|
|================================================================================================== | 81%
|
|=================================================================================================== | 81%
|
|=================================================================================================== | 82%
|
|==================================================================================================== | 82%
|
|===================================================================================================== | 82%
|
|===================================================================================================== | 83%
|
|====================================================================================================== | 83%
|
|====================================================================================================== | 84%
|
|======================================================================================================= | 84%
|
|======================================================================================================= | 85%
|
|======================================================================================================== | 85%
|
|======================================================================================================== | 86%
|
|========================================================================================================= | 86%
|
|========================================================================================================== | 86%
|
|========================================================================================================== | 87%
|
|=========================================================================================================== | 87%
|
|=========================================================================================================== | 88%
|
|============================================================================================================ | 88%
|
|============================================================================================================ | 89%
|
|============================================================================================================= | 89%
|
|============================================================================================================= | 90%
|
|============================================================================================================== | 90%
|
|============================================================================================================== | 91%
|
|=============================================================================================================== | 91%
|
|================================================================================================================ | 91%
|
|================================================================================================================ | 92%
|
|================================================================================================================= | 92%
|
|================================================================================================================= | 93%
|
|================================================================================================================== | 93%
|
|================================================================================================================== | 94%
|
|=================================================================================================================== | 94%
|
|=================================================================================================================== | 95%
|
|==================================================================================================================== | 95%
|
|===================================================================================================================== | 96%
|
|====================================================================================================================== | 96%
|
|====================================================================================================================== | 97%
|
|======================================================================================================================= | 97%
|
|======================================================================================================================= | 98%
|
|======================================================================================================================== | 98%
|
|======================================================================================================================== | 99%
|
|========================================================================================================================= | 99%
|
|========================================================================================================================= | 100%
|
|==========================================================================================================================| 100%
Centering and scaling data matrix
|
| | 0%
|
|============================== | 25%
|
|============================================================= | 50%
|
|============================================================================================ | 75%
|
|==========================================================================================================================| 100%
Place corrected count matrix in counts slot
Set default assay to SCT
Variables_genes <- L4@assays$SCT@var.features
# Exclude genes starting with "HLA-" AND "Xist" AND "TRBV, TRAV"
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
L4 <- RunPCA(L4,
features = Variables_genes_after_exclusion,
do.print = TRUE,
pcs.print = 1:5,
genes.print = 15,
npcs = 50)
PC_ 1
Positive: HSPD1, HSP90AB1, SRM, NPM1, NCL, SERBP1, HSPA9, CCT6A, RPS17, HNRNPAB
HSPE1, SET, C1QBP, YBX1, RANBP1, ATP5MC1, RBM17, ODC1, MTDH, RPL23
PRKDC, BATF3, CCT8, NME2, PAICS, NME1, VDAC1, TXNL4A, RAN, TRAP1
Negative: LGALS1, S100A6, S100A4, S100A11, B2M, LSP1, S100A10, SH3BGRL3, TMSB10, CRIP1
CYBA, MYL6, LAPTM5, IL32, CAPG, PRDX5, TAGLN2, ANXA2, FXYD5, IFITM2
TMSB4X, TIMP1, IFITM1, ANXA1, RHOC, VIM, CD63, IL9R, CARS, TNFRSF18
PC_ 2
Positive: RABGAP1L, FN1, ATP8A1, DENND4A, CAMK4, CABLES1, RORA, RNGTT, TNFRSF11A, APP
LINC01934, FNBP1, PCBP3, MAPK8, TMEM178B, ATP8B4, MCTP2, ZNF292, CCR7, CAMK1D
PDE4D, THEMIS, NEDD4L, OSBPL10, CDC42BPA, MACF1, AKAP13, RERE, MALAT1, ELMO1
Negative: CORO1A, CYCS, RPL35, RAN, DYNLL1, SNRPD1, UBE2S, RPL17, PSMB8, PA2G4
CDC20, CLIC1, PRELID1, HSPE1, NOP16, NPM1, TXN, ARHGDIB, SSBP1, PTTG1
ODC1, CCDC85B, MT1E, CCT2, PDCD5, IFITM2, TOMM40, HNRNPAB, CALM1, NOL7
PC_ 3
Positive: KCNQ5, AHNAK, RASGRP2, MKI67, PRUNE2, RRM2, LMNB1, ARHGEF6, SLC1A5, PXYLP1
CORO1A, DYNC1H1, RPS6KA5, HIST1H1B, FLNA, MYH9, ANTXR2, AAK1, KLF2, TPM4
ATAD2, MAD2L2, SMC4, GNAQ, CTDSPL, SPATA5, NSD2, UBAC1, UBASH3B, TUBA1B
Negative: RPL19, NME2, CA2, CCR7, DEGS1, TOMM20, MATN4, RPL27, SNRPE, CD74
SNHG16, SLC2A3, MINDY3, ARF1, ARID5A, JPT1, ANKRD37, PHB, PNRC1, PSMB3
NENF, TNFRSF11A, SEC14L1, RGS1, NME1, PKM, RPL22, MT-ND4, RPL38, RPS3A
PC_ 4
Positive: HSPA8, DDX21, CDC20, RAB11FIP1, CELF2, NCL, HSPA4, CD37, HSPH1, NAMPT
AK4, SLC2A3, HSPA5, JUN, GSPT1, NOP58, NOLC1, SERBP1, GTPBP4, ATP12A
NUDCD1, DDX17, PRR13, CORO1B, NOP16, SNHG3, DBN1, RSL1D1, DCUN1D5, CAPRIN1
Negative: HIST1H4C, TUBB, H2AFZ, TUBA1B, RRM2, RPL19, NME2, PPIA, GGH, H3F3A
STMN1, NUSAP1, NENF, TOP2A, EIF1, PCLAF, H1FX, HIST1H1D, HIST1H1C, HIST1H1B
EIF4EBP1, CDK1, CENPU, ZWINT, TYMS, PSAT1, H3F3B, ATAD2, BATF3, HIST1H1E
PC_ 5
Positive: TRAF3IP3, LIMD2, DANCR, RPL23, FRMPD4, JPT1, MT-ND1, BATF3, ARL6IP5, ZNF804A
CD52, SNRPE, AGBL1, ALAS1, NLGN1, TNFSF10, SLFN5, CD47, AC004687.1, LINC00271
ALOX5AP, TSHZ2, WFDC1, SH3BGRL3, GRIA4, CISD3, SLC20A2, RNGTT, RPL19, NDUFC2
Negative: GAPDH, LDHA, ENO1, MIF, TPI1, SLC2A3, DDIT4, PLIN2, ANKRD37, RGS1
PPP1R15A, BNIP3, HILPDA, BNIP3L, MXI1, DUSP4, DNAJB1, HSPA1B, MT-ND4, KLF6
C12orf75, HIST1H1E, RPL41, P4HA1, BHLHE40, N4BP3, TNFAIP3, JUN, AK4, GPI
# determine dimensionality of the data
ElbowPlot(L4, ndims =50)
NA
NA
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 L4$cell_line is a factor or character vector containing cell line names
data <- as.data.frame(table(L4$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 <- L4[["pca"]]@stdev / sum(L4[["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] 14
# TEST-2
# get significant PCs
stdv <- L4[["pca"]]@stdev
sum.stdv <- sum(L4[["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] 14
# 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()
NA
NA
NA
L4 <- FindNeighbors(L4,
dims = 1:min.pc,
verbose = FALSE)
# understanding resolution
L4 <- FindClusters(L4,
resolution = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6,
0.7,0.8, 0.9, 1, 1.1, 1.2))
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9423
Number of communities: 5
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9061
Number of communities: 6
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8818
Number of communities: 7
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8615
Number of communities: 9
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8458
Number of communities: 11
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8312
Number of communities: 10
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8157
Number of communities: 11
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8017
Number of communities: 11
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7899
Number of communities: 14
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7812
Number of communities: 14
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7720
Number of communities: 14
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 6150
Number of edges: 209562
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7631
Number of communities: 15
Elapsed time: 0 seconds
# non-linear dimensionality reduction --------------
L4 <- RunUMAP(L4,
dims = 1:min.pc,
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 Label Clusters function to help label
# individual clusters
DimPlot(L4,group.by = "cell_line",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.1",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.2",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.3",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.4",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.5",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.6",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.7",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.8",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.0.9",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.1",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.1.1",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4,
group.by = "SCT_snn_res.1.2",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
InstallData("pbmcref")
Warning: The following packages are already installed and will not be reinstalled: pbmcref
# The RunAzimuth function can take a Seurat object as input
L4 <- RunAzimuth(L4, reference = "pbmcref")
Warning: Overwriting miscellanous data for modelWarning: Adding a dimensional reduction (refUMAP) without the associated assay being presentWarning: Adding a dimensional reduction (refUMAP) without the associated assay being presentdetected 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 553 anchors
Finding integration vectors
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
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 ('-')
| | 0 % ~calculating
Integrating dataset 2 with reference dataset
Finding integration vectors
Integrating data
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s
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: Number of neighbors between query and reference is not equal to the number of neighbors within reference15:12:57 Read 6150 rows
15:12:57 Processing block 1 of 1
15:12:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
15:12:57 Initializing by weighted average of neighbor coordinates using 1 thread
15:12:57 Commencing optimization for 67 epochs, with 123000 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:12:57 Finished
Warning: No assay specified, setting assay as RNA by default.Projecting reference PCA onto query
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Projecting back the query cells into original PCA space
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Computing scores:
Finding neighbors of original query cells
Finding neighbors of transformed query cells
Computing query SNN
Determining bandwidth and computing transition probabilities
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Total elapsed time: 2.42626547813416
DimPlot(L4, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot (L4, group.by = "predicted.celltype.l2",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
#Load reference atlas and query data
ref <- readRDS(file = "CD4T_human_ref_v1.rds")
#Run Projection algorithm
query.projected <- Run.ProjecTILs(L4, ref = ref)
|
| | 0%[1] "Using assay SCT for query"
Pre-filtering cells with scGate...
### Detected a total of 4321 pure 'Target' cells (70.26% of total)
[1] "1829 out of 6150 ( 30% ) non-pure cells removed. Use filter.cells=FALSE to avoid pre-filtering"
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Aligning query to reference map for batch-correction..."
Warning: Layer counts isn't present in the assay object[[assay]]; returning NULLPreparing PCA embeddings for objects...
Warning: Number of dimensions changing from 50 to 20
| | 0 % ~calculating
|+++++++++++++++++++++++++ | 50% ~04s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
| | 0 % ~calculating
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s
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): GZMK, CCL20, MT1G, CD177, CGA, G0S2, IL17A, CCL4L2, XCL1, IGFL2, GZMH, TRDC, XCL2, IL1RL1, ACTG2, RASD1, KLRD1, ZBED2, CXCR6, TASL, CPA5, TMIGD2, H1-4, IL26, H2AZ1, GPR25, ELAPOR1, CCL3L3, POLR1F, EGR3, KIR3DL2, AHSP, CDKN2B, LIMS2, H1-2, CDKN2A, H1-0, WARS1, H1-3, ASCL2, DTHD1, H2BC11, GPX1, H2AC6, IRAG2, H1-10, MYO7A, FASLG, SCML1, CLEC7A, H3C10, NAP1L2, HS3ST3B1, FBLN7, FGFBP2, IL22, SLC28A3, PDLIM4, ZNF683, ECEL1, ARC, NLRP3, GFPT2, H4C3, RORC, GSTM2
|
|==========================================================================================================================| 100%
Creating slots functional.cluster and functional.cluster.conf in query object
#reference atlas
DimPlot(ref, label = T)
#Visualize projection
plot.projection(ref, query.projected, linesize = 0.5, pointsize = 0.5)
#Plot the predicted composition of the query in terms of reference T cell subtypes
plot.statepred.composition(ref, query.projected, metric = "Percent")
L4 <- ProjecTILs.classifier(query = L4, ref = ref)
|
| | 0%[1] "Using assay SCT for query"
Pre-filtering cells with scGate...
### Detected a total of 4321 pure 'Target' cells (70.26% of total)
[1] "1829 out of 6150 ( 30% ) non-pure cells removed. Use filter.cells=FALSE to avoid pre-filtering"
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
[1] "Aligning query to reference map for batch-correction..."
Warning: Layer counts isn't present in the assay object[[assay]]; returning NULLPreparing PCA embeddings for objects...
Warning: Number of dimensions changing from 50 to 20
| | 0 % ~calculating
|+++++++++++++++++++++++++ | 50% ~04s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
| | 0 % ~calculating
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=11s
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): GZMK, CCL20, MT1G, CD177, CGA, G0S2, IL17A, CCL4L2, XCL1, IGFL2, GZMH, TRDC, XCL2, IL1RL1, ACTG2, RASD1, KLRD1, ZBED2, CXCR6, TASL, CPA5, TMIGD2, H1-4, IL26, H2AZ1, GPR25, ELAPOR1, CCL3L3, POLR1F, EGR3, KIR3DL2, AHSP, CDKN2B, LIMS2, H1-2, CDKN2A, H1-0, WARS1, H1-3, ASCL2, DTHD1, H2BC11, GPX1, H2AC6, IRAG2, H1-10, MYO7A, FASLG, SCML1, CLEC7A, H3C10, NAP1L2, HS3ST3B1, FBLN7, FGFBP2, IL22, SLC28A3, PDLIM4, ZNF683, ECEL1, ARC, NLRP3, GFPT2, H4C3, RORC, GSTM2
|
|==========================================================================================================================| 100%
Creating slots functional.cluster and functional.cluster.conf in query object
DimPlot(L4, group.by = "functional.cluster",
reduction = "umap",
label.size = 3,
repel = T,
label = T)
NA
NA
#get reference datasets from celldex package
monaco.ref <- celldex::MonacoImmuneData()
snapshotDate(): 2023-10-24
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
hpca.ref <- celldex::HumanPrimaryCellAtlasData()
snapshotDate(): 2023-10-24
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
dice.ref <- celldex::DatabaseImmuneCellExpressionData()
snapshotDate(): 2023-10-24
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
bpe.ref <- celldex::BlueprintEncodeData()
snapshotDate(): 2023-10-24
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
#convert our Seurat object to single cell experiment (SCE)
sce <- as.SingleCellExperiment(DietSeurat(L4))
#using SingleR
monaco.main <- SingleR(test = sce,assay.type.test = 1,ref = monaco.ref,labels = monaco.ref$label.main)
monaco.fine <- SingleR(test = sce,assay.type.test = 1,ref = monaco.ref,labels = monaco.ref$label.fine)
hpca.main <- SingleR(test = sce,assay.type.test = 1,ref = hpca.ref,labels = hpca.ref$label.main)
hpca.fine <- SingleR(test = sce,assay.type.test = 1,ref = hpca.ref,labels = hpca.ref$label.fine)
dice.main <- SingleR(test = sce,assay.type.test = 1,ref = dice.ref,labels = dice.ref$label.main)
dice.fine <- SingleR(test = sce,assay.type.test = 1,ref = dice.ref,labels = dice.ref$label.fine)
bpe.main <- SingleR(test = sce,assay.type.test = 1,ref = bpe.ref,labels = bpe.ref$label.main)
bpe.fine <- SingleR(test = sce,assay.type.test = 1,ref = bpe.ref,labels = bpe.ref$label.fine)
#summary of general cell type annotations
#table(monaco.main$pruned.labels)
#table(hpca.main$pruned.labels)
#table(dice.main$pruned.labels)
#table(bpe.main$pruned.labels)
#The finer cell types annotations are you after, the harder they are to get reliably.
#This is where comparing many databases, as well as using individual markers from literature,
#would all be very valuable.
#table(monaco.fine$pruned.labels)
#table(hpca.fine$pruned.labels)
#table(dice.fine$pruned.labels)
#table(bpe.fine$pruned.labels)
#add the annotations to the Seurat object metadata
L4@meta.data$monaco.main <- monaco.main$pruned.labels
L4@meta.data$monaco.fine <- monaco.fine$pruned.labels
#
L4@meta.data$hpca.main <- hpca.main$pruned.labels
L4@meta.data$hpca.fine <- hpca.fine$pruned.labels
#
L4@meta.data$dice.main <- dice.main$pruned.labels
L4@meta.data$dice.fine <- dice.fine$pruned.labels
#
L4@meta.data$bpe.main <- bpe.main$pruned.labels
L4@meta.data$bpe.fine <- bpe.fine$pruned.labels
# Monaco Annotations
DimPlot(L4, group.by = "monaco.main",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(L4, group.by = "monaco.main",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4, group.by = "monaco.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(L4, group.by = "monaco.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
# HPCA Annotations
DimPlot(L4, group.by = "hpca.main",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4, group.by = "hpca.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(L4, group.by = "hpca.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
# DICE Annotations
DimPlot(L4, group.by = "dice.main",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4, group.by = "dice.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(L4, group.by = "dice.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
# BPE Annotations
DimPlot(L4, group.by = "bpe.main",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
DimPlot(L4, group.by = "bpe.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = F)
DimPlot(L4, group.by = "bpe.fine",
reduction = "umap",
label.size = 3,
repel = T,
label = T, label.box = T)
NA
NA
NA
library(clustree)
Loading required package: ggraph
Attaching package: 'ggraph'
The following object is masked from 'package:sp':
geometry
clustree(L4, prefix = "SCT_snn_res.")
save(L4, file = "L4_Analysis.Robj")
```