Merged All samples with PBMC_10x and removed non CD4 T cells from Control and B cells from L4 regress nCount, percent.mt and rb and apply SCT

1. load libraries

Loading required package: SeuratObject

Loading required package: sp

Attaching package: 'SeuratObject'

The following objects are masked from 'package:base':

    intersect, t

── Installed datasets ──────────────────────────────── SeuratData v0.2.2.9001 ──

✔ pbmcref 1.0.0                         ✔ pbmcsca 3.0.0

────────────────────────────────────── Key ─────────────────────────────────────

✔ Dataset loaded successfully
❯ Dataset built with a newer version of Seurat than installed
❓ Unknown version of Seurat installed

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

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ ggplot2   3.5.1     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'magrittr'


The following object is masked from 'package:purrr':

    set_names


The following object is masked from 'package:tidyr':

    extract



Attaching package: 'dbplyr'


The following objects are masked from 'package:dplyr':

    ident, sql


Registered S3 method overwritten by 'SeuratDisk':
  method            from  
  as.sparse.H5Group Seurat



Attaching shinyBS

Loading required package: ggraph


Attaching package: 'ggraph'


The following object is masked from 'package:sp':

    geometry

2. Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
load("../0-Robj/SS_CD4_Tcells_PBMC10x_final_for_SCT_AzimuthAnnotation_and_Integration.robj")

# Azimuth l1
janitor::tabyl(filtered_seurat@meta.data, predicted.celltype.l1, cell_line)

# Azimuth l2
janitor::tabyl(filtered_seurat@meta.data, predicted.celltype.l2, cell_line)

# Azimuth l3
janitor::tabyl(filtered_seurat@meta.data, predicted.celltype.l3, cell_line)

All_samples_Merged <- filtered_seurat

All_samples_Merged

An object of class Seurat 
36752 features across 49388 samples within 5 assays 
Active assay: RNA (36601 features, 0 variable features)
 2 layers present: data, counts
 4 other assays present: ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 2 dimensional reductions calculated: integrated_dr, ref.umap

Summarizing Seurat Object

# Load necessary libraries
library(Seurat)

# Display basic metadata summary
head(All_samples_Merged@meta.data)

# Check if columns such as `orig.ident`, `nCount_RNA`, `nFeature_RNA`, `nUMI`, `ngene`, and any other necessary columns exist
required_columns <- c("orig.ident", "nCount_RNA", "nFeature_RNA", "nUMI", "ngene")
missing_columns <- setdiff(required_columns, colnames(All_samples_Merged@meta.data))

if (length(missing_columns) > 0) {
    cat("Missing columns:", paste(missing_columns, collapse = ", "), "\n")
} else {
    cat("All required columns are present.\n")
}

All required columns are present.

# Check cell counts and features
cat("Number of cells:", ncol(All_samples_Merged), "\n")

Number of cells: 49388 

cat("Number of features:", nrow(All_samples_Merged), "\n")

Number of features: 36601 

# Verify that each `orig.ident` label has the correct number of cells
cat("Cell counts per group:\n")

Cell counts per group:

print(table(All_samples_Merged$orig.ident))

     L1      L2      L3      L4      L5      L6      L7    PBMC PBMC10x 
   5825    5935    6428    6023    6022    5148    5331    5171    3505 

# Check that the cell IDs are unique (which ensures no issues from merging)
if (any(duplicated(colnames(All_samples_Merged)))) {
    cat("Warning: There are duplicated cell IDs.\n")
} else {
    cat("Cell IDs are unique.\n")
}

Cell IDs are unique.

# Check the assay consistency for RNA
DefaultAssay(All_samples_Merged) <- "RNA"

# Check dimensions of the RNA counts layer using the new method
cat("Dimensions of the RNA counts layer:", dim(GetAssayData(All_samples_Merged, layer = "counts")), "\n")

Dimensions of the RNA counts layer: 36601 49388 

cat("Dimensions of the RNA data layer:", dim(GetAssayData(All_samples_Merged, layer = "data")), "\n")

Dimensions of the RNA data layer: 36601 49388 

# Check the ADT assay (optional)
if ("ADT" %in% names(All_samples_Merged@assays)) {
    cat("ADT assay is present.\n")
    cat("Dimensions of the ADT counts layer:", dim(GetAssayData(All_samples_Merged, assay = "ADT", layer = "counts")), "\n")
} else {
    cat("ADT assay is not present.\n")
}

ADT assay is present.
Dimensions of the ADT counts layer: 56 49388 

3. QC

# Remove the percent.mito column
All_samples_Merged$percent.mito <- NULL


# Set identity classes to an existing column in meta data
Idents(object = All_samples_Merged) <- "cell_line"

All_samples_Merged[["percent.rb"]] <- PercentageFeatureSet(All_samples_Merged, 
                                                           pattern = "^RP[SL]")
# Convert 'percent.mt' to numeric, replacing "NaN" with 0
All_samples_Merged$percent.rb <- replace(as.numeric(All_samples_Merged$percent.rb), is.na(All_samples_Merged$percent.rb), 0)



# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
All_samples_Merged[["percent.mt"]] <- PercentageFeatureSet(All_samples_Merged, pattern = "^MT-")

# Convert 'percent.mt' to numeric, replacing "NaN" with 0
All_samples_Merged$percent.mt <- replace(as.numeric(All_samples_Merged$percent.mt), is.na(All_samples_Merged$percent.mt), 0)





VlnPlot(All_samples_Merged, features = c("nFeature_RNA", 
                                         "nCount_RNA", 
                                         "percent.mt",
                                         "percent.rb"), 
                            ncol = 4, pt.size = 0.1) & 
              theme(plot.title = element_text(size=10))

FeatureScatter(All_samples_Merged, feature1 = "percent.mt", 
                                  feature2 = "percent.rb")

VlnPlot(All_samples_Merged, features = c("nFeature_RNA", 
                                    "nCount_RNA", 
                                    "percent.mt"), 
                                      ncol = 3)

FeatureScatter(All_samples_Merged, 
               feature1 = "percent.mt", 
               feature2 = "percent.rb") +
        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'

##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.mt")+
  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'

Assign Cell-Cycle Scores

Running SCTransform on assay: RNA

Running SCTransform on layer: counts

vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.

Variance stabilizing transformation of count matrix of size 26179 by 49388

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

Found 487 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 26179 genes

Computing corrected count matrix for 26179 genes

Calculating gene attributes

Wall clock passed: Time difference of 8.04581 mins

Determine variable features

Getting residuals for block 1(of 10) for counts dataset

Getting residuals for block 2(of 10) for counts dataset

Getting residuals for block 3(of 10) for counts dataset

Getting residuals for block 4(of 10) for counts dataset

Getting residuals for block 5(of 10) for counts dataset

Getting residuals for block 6(of 10) for counts dataset

Getting residuals for block 7(of 10) for counts dataset

Getting residuals for block 8(of 10) for counts dataset

Getting residuals for block 9(of 10) for counts dataset

Getting residuals for block 10(of 10) for counts dataset

Finished calculating residuals for counts

Set default assay to SCT

Warning: The following features are not present in the object: MLF1IP, not
searching for symbol synonyms

Warning: The following features are not present in the object: FAM64A, HN1, not
searching for symbol synonyms

4. Normalize data

# Apply SCTransform
All_samples_Merged <- SCTransform(All_samples_Merged, 
                                  vars.to.regress = c("percent.rb","percent.mt", "nCount_RNA"), 
                                  do.scale=TRUE, 
                                  do.center=TRUE, 
                                  verbose = TRUE)

Running SCTransform on assay: RNA

Running SCTransform on layer: counts

vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.

Variance stabilizing transformation of count matrix of size 26179 by 49388

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

Found 487 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 26179 genes

Computing corrected count matrix for 26179 genes

Calculating gene attributes

Wall clock passed: Time difference of 6.396079 mins

Determine variable features

Regressing out percent.rb, percent.mt, nCount_RNA

Centering and scaling data matrix

Getting residuals for block 1(of 10) for counts dataset

Getting residuals for block 2(of 10) for counts dataset

Getting residuals for block 3(of 10) for counts dataset

Getting residuals for block 4(of 10) for counts dataset

Getting residuals for block 5(of 10) for counts dataset

Getting residuals for block 6(of 10) for counts dataset

Getting residuals for block 7(of 10) for counts dataset

Getting residuals for block 8(of 10) for counts dataset

Getting residuals for block 9(of 10) for counts dataset

Getting residuals for block 10(of 10) for counts dataset

Regressing out percent.rb, percent.mt, nCount_RNA

Centering and scaling data matrix

Finished calculating residuals for counts

Set default assay to SCT

5. Perform PCA

Variables_genes <- All_samples_Merged@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)]

# Set the seed for clustering steps
set.seed(123)

# 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:  CD7, PRKCH, KIR3DL1, SEPTIN9, PTPRC, KIR2DL3, CLEC2B, CD52, ARHGAP15, CST7 
       MALAT1, CD3G, EPCAM, RPS27, ESYT2, XCL1, MATK, GZMM, LEF1, TRGV2 
       CD6, TC2N, MYO1E, KLRC1, KIR2DL4, KIR3DL2, KLRK1, CXCR3, LCK, PTPN6 
Negative:  NPM1, SEC11C, YBX3, VDAC1, MTHFD2, MTDH, CCT8, IL2RA, HDGFL3, RBM17 
       PRELID1, C12orf75, RAN, PRDX1, CCND2, HINT2, BATF3, MIR155HG, HSP90AB1, KRT7 
       SPATS2L, GAPDH, SRM, HSPD1, HTATIP2, CD74, CANX, PKM, MINDY3, SLC35F3 
PC_ 2 
Positive:  C12orf75, CYBA, HACD1, LY6E, SCCPDH, EGFL6, TNFRSF4, ATP5MC1, APRT, ENO1 
       BACE2, ARPC2, TIGIT, GGH, PTP4A3, SYT4, CCL17, SPINT2, CHCHD2, CORO1B 
       PON2, RPL27A, CTSC, COX6A1, GYPC, NME2, NET1, NME1, PLPP1, RHOC 
Negative:  PAGE5, RPL35A, RBPMS, CD74, NDUFV2, TENM3, LMNA, RPL22L1, CDKN2A, KIF2A 
       RPS3A, RPL11, PSMB9, ANXA5, PLD1, PPP2R2B, FAM241A, SPOCK1, B2M, VAMP5 
       STAT1, FAM50B, SH3KBP1, ERAP2, ZC2HC1A, GPX4, IFI27L2, RPS14, MSC-AS1, CTAG2 
PC_ 3 
Positive:  RPL30, RPL39, RPS27, RPS4Y1, MT-ND3, ETS1, MALAT1, BTG1, RPS29, TPT1 
       TCF7, RPL34, FYB1, ZBTB20, SELL, ANK3, SARAF, FAM107B, IL7R, LINC00861 
       TXNIP, CSGALNACT1, RIPOR2, PNRC1, PIK3IP1, TIGIT, EEF1A2, THEMIS, ATP8B4, LINC01934 
Negative:  PFN1, NME2, KIR3DL2, RPS15, EIF4A1, C1QBP, MIF, NDUFA4, CHCHD2, KIR2DL3 
       RPL19, ATP5MC3, ACTB, EIF5A, KIR3DL1, RPL27A, CLIC1, MT-CO2, CST7, HMGN2 
       DAD1, COX6A1, GAPDH, EPCAM, TRGV2, PSMB6, RPS2, RAB25, GGCT, WDR34 
PC_ 4 
Positive:  HSPE1, EIF5A, ATP5MC3, RPL34, RPS4Y1, MT-ND3, ODC1, CHCHD10, CYCS, CYC1 
       HSPD1, GCSH, RPL39, FKBP4, PPBP, PPID, RPS29, HSP90AA1, TCF7, FCER2 
       TOMM40, GSTP1, RPL30, CD7, FKBP11, DNAJC12, FAM162A, C1QBP, ATP5F1B, PRELID3B 
Negative:  RPS4X, GAS5, EGLN3, KRT1, LINC02752, WFDC1, TTC29, TBX4, RPLP1, RPL13 
       AC069410.1, IFNGR1, IL32, PLCB1, TNS4, SP5, S100A11, FAM9C, SEMA4A, IL4 
       S100A4, S100A6, NKG7, LINC00469, VIM, HSPB1, CEBPD, RPLP0, VIPR2, SOCS1 
PC_ 5 
Positive:  TMSB4X, LGALS1, TMSB10, S100A11, S100A4, S100A6, COTL1, IFITM2, LSP1, TP73 
       TAGLN2, TMEM163, GPAT3, LIME1, HOXC9, CRIP1, GAS2L1, LAPTM5, DUSP4, TNFRSF18 
       GPAT2, EMP3, IFITM1, EEF1A2, MIIP, QPRT, PRDX5, CARHSP1, ACTB, RBM38 
Negative:  CCL17, MIR155HG, MAP4K4, LRBA, PRKCA, RUNX1, MYO1D, RXFP1, IMMP2L, CA10 
       CFI, DOCK10, CA2, FRMD4A, AL590550.1, NFIB, THY1, EZH2, LTA, SNTB1 
       SLC35F3, RANBP17, HS3ST1, IGHE, CCL5, NME2, AKAP12, DENND4A, AC100801.1, MGST3 

# determine dimensionality of the data
ElbowPlot(All_samples_Merged, ndims = 50)

6. 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] 17

# 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] 17

# 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()

7. Clustering

# Set the seed for clustering steps
set.seed(123)

All_samples_Merged <- FindNeighbors(All_samples_Merged, 
                                dims = 1:min.pc, 
                                verbose = FALSE)

# understanding resolution
All_samples_Merged <- FindClusters(All_samples_Merged, 
                                    resolution = c(0.4, 0.5, 0.6, 0.7,0.8))

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49388
Number of edges: 1638914

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9557
Number of communities: 14
Elapsed time: 16 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49388
Number of edges: 1638914

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9462
Number of communities: 15
Elapsed time: 18 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49388
Number of edges: 1638914

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9372
Number of communities: 18
Elapsed time: 17 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49388
Number of edges: 1638914

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9287
Number of communities: 18
Elapsed time: 17 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 49388
Number of edges: 1638914

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9210
Number of communities: 22
Elapsed time: 17 seconds

UMAP Visualization

# Set the seed for clustering steps
set.seed(123)

# non-linear dimensionality reduction --------------
All_samples_Merged <- RunUMAP(All_samples_Merged, 
                          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(All_samples_Merged,group.by = "cell_line", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.4", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.5", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.6", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.7", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.8", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

# Set identity classes to an existing column in meta data
Idents(object = All_samples_Merged) <- "SCT_snn_res.0.4"

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)))

print(cluster_table)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
6405 5978 5930 5871 5272 5108 5075 4102 3342 1751  216  196   72   70 

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.4)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11
  B intermediate       0    3    0    0    0    2    0    1    0    0    1    0
  B memory             8    7    0    1  116   79    0   32    0    2    7    0
  CD14 Mono            0    0    0    2    6    0    0    5    0    0    0    0
  CD4 CTL              0    0    0    0    0    0   12    0    0    0    0    0
  CD4 Naive            0    8    0    0    0    0  521    0 1479    0    0   33
  CD4 Proliferating 5448 2475 2852 5323 4101 3974    0 3256    6 1353  158    0
  CD4 TCM            873 3414  268  517  486  568 4480  109 1835   42   18  161
  CD4 TEM              0    1    0    0    0    0   61    0   22    0    0    0
  CD8 Proliferating    0    0    0    0    1    1    0    0    0    0    0    0
  CD8 TCM              0    1   16    0    0    0    0    0    0    0    0    0
  CD8 TEM              0    1    8    0    1    3    0    2    0    0    0    0
  cDC1                 0    0    0    0    0    5    0    2    0    0    1    0
  cDC2                 0    0    0    2   36    3    0   10    0    0    0    0
  dnT                  0    3    1    1    3    4    0    2    0    0    0    1
  HSPC                57   10    0    1  486  204    0  673    0  353   21    0
  ILC                  0    1    0    0    0    0    0    0    0    0    0    0
  NK                   0    0    0    0    0    0    0    0    0    0    0    0
  NK Proliferating     4   40 2785   23   36  252    0   10    0    1   10    0
  Treg                15   14    0    1    0   13    1    0    0    0    0    1
                   
                      12   13
  B intermediate       0    0
  B memory             0    0
  CD14 Mono           13    0
  CD4 CTL              1    0
  CD4 Naive            1    0
  CD4 Proliferating    0   65
  CD4 TCM             54    5
  CD4 TEM              0    0
  CD8 Proliferating    0    0
  CD8 TCM              0    0
  CD8 TEM              0    0
  cDC1                 0    0
  cDC2                 2    0
  dnT                  0    0
  HSPC                 0    0
  ILC                  0    0
  NK                   1    0
  NK Proliferating     0    0
  Treg                 0    0

8. clusTree

# clustree(All_samples_Merged, prefix = "SCT_snn_res.")

9. Azimuth Annotation

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
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 5126 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

18:42:57 Read 49388 rows

18:42:57 Processing block 1 of 1

18:42:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
18:42:57 Initializing by weighted average of neighbor coordinates using 1 thread
18:42:57 Commencing optimization for 67 epochs, with 987760 positive edges
18:43:14 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: 39.0926122665405

10. 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)

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 = F)

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.4)

                   
                       0    1    2    3    4    5    6    7    8    9   10   11
  ASDC                24    0    0    9    6    2    0    1    0    0    0    0
  B intermediate       6    8    0    0    1    4    0    0    0    0    1    0
  B memory             4   17    0    0    4    2    7    0    0    0    0    0
  B naive              4   64    0    0    1    2   22    0    1    1    0    0
  CD14 Mono            0    0    0    0    0    0    0    0    0    0    0    0
  CD16 Mono            0    0    0    0    0    0    0    0    0    0    0    0
  CD4 CTL              0    0    0    0    0    0   11    0    0    0    0    0
  CD4 Naive            0    1    0    0    0    0  416    0 1447    0    0   35
  CD4 Proliferating 5263 2444 2893 5306 3909 3803    0 2998    5 1248  151    0
  CD4 TCM            645 3327  234  508  268  123 4552   83 1871   13    4  161
  CD4 TEM              0    2    0    0    0    0   66    0   18    0    0    0
  CD8 Proliferating    1    0    0    0    0    2    0    0    0    0    0    0
  CD8 TCM              0   73    0    0    0    0    0    0    0    0    0    0
  CD8 TEM              0    1   70    0    0    1    1    0    0    0    0    0
  cDC2               332    7    0   39  464  753    0  198    0   76   31    0
  HSPC               123    9    0    2  601  222    0  816    0  412   21    0
  ILC                  0    1    0    0    0    0    0    0    0    0    0    0
  NK Proliferating     3   24 2733    7   18  194    0    6    0    1    8    0
                   
                      12   13
  ASDC                 0    1
  B intermediate       0    0
  B memory             0    0
  B naive              4    0
  CD14 Mono           14    0
  CD16 Mono           15    0
  CD4 CTL              1    0
  CD4 Naive            1    0
  CD4 Proliferating    0   64
  CD4 TCM             31    0
  CD4 TEM              0    0
  CD8 Proliferating    0    0
  CD8 TCM              0    0
  CD8 TEM              5    0
  cDC2                 1    5
  HSPC                 0    0
  ILC                  0    0
  NK Proliferating     0    0

Save the Seurat object as an Robj file

save(All_samples_Merged, file = "../0-Robj/CD4Tcells_SCTnormalized_and_Azimuth_Annotation_done_on_HPC.robj")