Merged All samples with PBMC_10x and removed non CD4 T cells from Control and B cells from L4 and ILC, NK, CD14 Mono 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-R_Objects/SS_CD4_Tcells_Azimuth_Annotated_PBMC10x_final_for_SCT_and_Integration.robj")

All_samples_Merged <- filtered_seurat

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: 49372 

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    6007    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 49372 

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

Dimensions of the RNA data layer: 36601 49372 

# 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 49372 

Azimuth Annotation

# InstallData("pbmcref")
# 
# # The RunAzimuth function can take a Seurat object as input
# All_samples_Merged <- RunAzimuth(All_samples_Merged, reference = "pbmcref")

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 26176 by 49372

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

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

Second step: Get residuals using fitted parameters for 26176 genes

Computing corrected count matrix for 26176 genes

Calculating gene attributes

Wall clock passed: Time difference of 8.277261 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 26176 by 49372

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

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

Second step: Get residuals using fitted parameters for 26176 genes

Computing corrected count matrix for 26176 genes

Calculating gene attributes

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

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

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: 49372
Number of edges: 1611982

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

Number of nodes: 49372
Number of edges: 1611982

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

Number of nodes: 49372
Number of edges: 1611982

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

Number of nodes: 49372
Number of edges: 1611982

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

Number of nodes: 49372
Number of edges: 1611982

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9190
Number of communities: 19
Elapsed time: 13 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 
6403 5977 5939 5931 5317 5215 5073 4107 3344 1813  196   57 

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    2    0    0    0    0
  B memory             8    6    1    0  116   86    0   31    0    4    0    0
  CD14 Mono            0    1    0    0    7    0    0    4    0    0    0    0
  CD4 CTL              0    0    0    0    0    0   12    0    0    0    0    1
  CD4 Naive            0    8    0    0    0    0  521    0 1479    0   33    1
  CD4 Proliferating 5448 2474 5388 2852 4144 4037    0 3257    6 1405    0    0
  CD4 TCM            871 3414  522  269  484  584 4478  109 1838   46  161   54
  CD4 TEM              0    1    0    0    0    0   62    0   21    0    0    0
  CD8 Proliferating    0    0    0    0    1    1    0    0    0    0    0    0
  CD8 TCM              0    1    0   16    0    0    0    0    0    0    0    0
  CD8 TEM              0    1    0    8    1    3    0    2    0    0    0    0
  cDC1                 0    0    0    0    0    6    0    2    0    0    0    0
  cDC2                 0    1    2    0   36    3    0   10    0    0    0    1
  dnT                  0    3    1    1    3    4    0    2    0    0    1    0
  HSPC                57   10    1    0  489  213    0  678    0  357    0    0
  NK Proliferating     4   40   23 2785   36  262    0   10    0    1    0    0
  Treg                15   14    1    0    0   14    0    0    0    0    1    0

8. clusTree

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

9. Azimuth Annotation

# InstallData("pbmcref")
# 
# # The RunAzimuth function can take a Seurat object as input
# All_samples_Merged <- RunAzimuth(All_samples_Merged, reference = "pbmcref")

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
  B intermediate       0    3    0    0    0    2    0    2    0    0    0    0
  B memory             8    6    1    0  116   86    0   31    0    4    0    0
  CD14 Mono            0    1    0    0    7    0    0    4    0    0    0    0
  CD4 CTL              0    0    0    0    0    0   12    0    0    0    0    1
  CD4 Naive            0    8    0    0    0    0  521    0 1479    0   33    1
  CD4 Proliferating 5448 2474 5388 2852 4144 4037    0 3257    6 1405    0    0
  CD4 TCM            871 3414  522  269  484  584 4478  109 1838   46  161   54
  CD4 TEM              0    1    0    0    0    0   62    0   21    0    0    0
  CD8 Proliferating    0    0    0    0    1    1    0    0    0    0    0    0
  CD8 TCM              0    1    0   16    0    0    0    0    0    0    0    0
  CD8 TEM              0    1    0    8    1    3    0    2    0    0    0    0
  cDC1                 0    0    0    0    0    6    0    2    0    0    0    0
  cDC2                 0    1    2    0   36    3    0   10    0    0    0    1
  dnT                  0    3    1    1    3    4    0    2    0    0    1    0
  HSPC                57   10    1    0  489  213    0  678    0  357    0    0
  NK Proliferating     4   40   23 2785   36  262    0   10    0    1    0    0
  Treg                15   14    1    0    0   14    0    0    0    0    1    0

Save the Seurat object as an Robj file

#save(All_samples_Merged, file = "../0-R_Objects/CD4Tcells_SCTnormalized_done_on_HPC.robj")

11. Harmony Integration

# Load required libraries
library(Seurat)
library(harmony)

Loading required package: Rcpp

library(ggplot2)

# Run Harmony, adjusting for batch effect using "cell_line" or another grouping variable
All_samples_Merged <- RunHarmony(
  All_samples_Merged,
  group.by.vars = c("cell_line"),  # Replace with the metadata column specifying batch or cell line
  assay.use="SCT")

Transposing data matrix

Initializing state using k-means centroids initialization

Harmony 1/10

Harmony 2/10

Harmony 3/10

Harmony converged after 3 iterations

# Check results in harmony embeddings
harmony_embeddings <- Embeddings(All_samples_Merged, reduction = "harmony")
head(harmony_embeddings)

                       harmony_1  harmony_2   harmony_3    harmony_4
L1_AAACCTGAGGGCTTCC-1   1.948641  1.0484874 -0.00806881  -5.81445998
L1_AAACCTGGTGCAGGTA-1  -9.909179 -0.6214129  2.65110042 -10.38620322
L1_AAACCTGGTTAAAGTG-1  -9.933443 -1.3016055 -2.97527749  -6.22359137
L1_AAACCTGTCAGGTAAA-1   3.357574 -1.0140749 -6.82128648  -0.06995136
L1_AAACCTGTCCCTGACT-1   2.312686  0.7096574  1.68673081  -2.59094358
L1_AAACCTGTCCTTCAAT-1 -12.462611  0.9807047  9.26965937  -3.24161085
                        harmony_5 harmony_6 harmony_7  harmony_8  harmony_9
L1_AAACCTGAGGGCTTCC-1   2.5424795 0.3628653 -1.380085 -1.7012642 2.70955231
L1_AAACCTGGTGCAGGTA-1  -6.1565395 7.6147219 -4.462583  0.4974874 0.03019239
L1_AAACCTGGTTAAAGTG-1   0.1939102 3.6047253 -1.446832 -2.9123242 2.76017343
L1_AAACCTGTCAGGTAAA-1   2.6122609 1.0670470 -1.139194  0.5914508 0.09361857
L1_AAACCTGTCCCTGACT-1  -0.7888451 0.1279331 -1.067497 -1.2229540 2.27768388
L1_AAACCTGTCCTTCAAT-1 -13.3783475 3.1114342 -9.032892 -7.8000125 5.35600316
                      harmony_10 harmony_11 harmony_12 harmony_13 harmony_14
L1_AAACCTGAGGGCTTCC-1  5.5639822   3.356909 -0.4091590 -0.9810571 -0.2569345
L1_AAACCTGGTGCAGGTA-1  1.1001777  -1.470097  1.9795566 -2.0074855 -0.1735669
L1_AAACCTGGTTAAAGTG-1 -2.0128436  -3.019327 -3.1762648  1.8248146  2.0013316
L1_AAACCTGTCAGGTAAA-1 -2.1448504  -4.378139 -1.7137991  2.0660880  0.9943382
L1_AAACCTGTCCCTGACT-1  1.9488831   2.436063 -0.3790963  0.5073282  0.4521224
L1_AAACCTGTCCTTCAAT-1  0.3802278   4.246651 -1.8281061 -2.9236830  0.2534402
                      harmony_15  harmony_16 harmony_17 harmony_18 harmony_19
L1_AAACCTGAGGGCTTCC-1  -2.964799  0.09820109  1.2908460  1.1812929  0.3942512
L1_AAACCTGGTGCAGGTA-1   0.648601  0.68182775 -1.4236026  1.1082348  2.4290357
L1_AAACCTGGTTAAAGTG-1   3.622343  3.01020267 -0.4361311 -2.1964401  0.2251690
L1_AAACCTGTCAGGTAAA-1   1.158282 -0.15300011 -0.2161842  0.2028314 -0.3280452
L1_AAACCTGTCCCTGACT-1  -1.122470 -1.38757877  4.1277006 -3.5861213 -1.3269795
L1_AAACCTGTCCTTCAAT-1   2.081071 -0.12642277  0.0329925 -0.1628263  0.0139672
                      harmony_20 harmony_21 harmony_22 harmony_23 harmony_24
L1_AAACCTGAGGGCTTCC-1  2.1410368  1.4172562   4.148122 -2.6884267  4.1542966
L1_AAACCTGGTGCAGGTA-1 -2.8486478  2.9504390   5.915975 -1.6108688  0.2497149
L1_AAACCTGGTTAAAGTG-1 -0.3587958 -1.0043108   1.752189 -1.6077171  0.2970351
L1_AAACCTGTCAGGTAAA-1  0.6377344 -1.5956720  -2.625786  0.4315349 -2.0538848
L1_AAACCTGTCCCTGACT-1  2.4283149 -1.3512176   4.361970 -2.1512464  1.0155496
L1_AAACCTGTCCTTCAAT-1 -2.2630667 -0.0616899   4.644279 -1.1644565 -0.8620327
                      harmony_25 harmony_26 harmony_27  harmony_28 harmony_29
L1_AAACCTGAGGGCTTCC-1 -1.9753098 -0.6855722  0.1651495 -0.61223199  0.7309983
L1_AAACCTGGTGCAGGTA-1 -0.1450759 -2.0074265  0.6578202  1.16190690 -0.1538706
L1_AAACCTGGTTAAAGTG-1 -0.3928556 -1.8189002 -0.9134228 -1.44268285  0.8661139
L1_AAACCTGTCAGGTAAA-1 -0.6663157  0.5972026 -1.0971138  0.07522274  0.2697583
L1_AAACCTGTCCCTGACT-1  0.3343284  1.3049786  0.3305861 -0.43349992 -0.4754056
L1_AAACCTGTCCTTCAAT-1 -2.5028190 -0.9493836 -0.4230632  0.11739141  1.1017206
                       harmony_30 harmony_31 harmony_32 harmony_33 harmony_34
L1_AAACCTGAGGGCTTCC-1  1.07025712 -0.9370386  1.2876791  0.1932222  1.6984845
L1_AAACCTGGTGCAGGTA-1 -0.08871048  0.5716273 -0.9012281 -1.3277878 -4.2008927
L1_AAACCTGGTTAAAGTG-1  2.11187936 -0.9554580  3.2039421  2.1747311  2.6182738
L1_AAACCTGTCAGGTAAA-1  0.20415517 -0.9194567 -1.1883659 -0.2724933 -0.3280742
L1_AAACCTGTCCCTGACT-1  0.23015569 -0.3927227  0.5519720  0.9935055  2.8882997
L1_AAACCTGTCCTTCAAT-1  1.26190254 -2.8976018 -0.6587909  2.9940796 -0.9998876
                      harmony_35 harmony_36 harmony_37 harmony_38  harmony_39
L1_AAACCTGAGGGCTTCC-1  1.6291505 -0.4622979  0.9918008 -0.1095008 -0.09044753
L1_AAACCTGGTGCAGGTA-1 -1.3991633  0.9025145  0.5690480  0.7315200  1.04802111
L1_AAACCTGGTTAAAGTG-1  0.7861378 -0.7397459 -0.8360783 -0.8497256 -0.95682510
L1_AAACCTGTCAGGTAAA-1  0.2906318 -1.2332554  2.2572042  1.1422478  0.54278056
L1_AAACCTGTCCCTGACT-1  0.7024852  0.9226384  0.8747567 -1.3661143 -1.46986865
L1_AAACCTGTCCTTCAAT-1 -0.7107164  0.5236052  4.3380333  0.9728332 -1.39132047
                      harmony_40 harmony_41 harmony_42 harmony_43 harmony_44
L1_AAACCTGAGGGCTTCC-1  0.4395848  0.2294327   0.732884 -0.2133815 -0.1088725
L1_AAACCTGGTGCAGGTA-1 -0.2547631  2.1443439   1.523305 -1.0858964 -1.0545710
L1_AAACCTGGTTAAAGTG-1 -0.6769785  0.3808731  -2.526305  1.3170597 -0.7941223
L1_AAACCTGTCAGGTAAA-1 -0.2889691 -1.4221899   0.874801  2.1957545  0.4364966
L1_AAACCTGTCCCTGACT-1  1.4626448 -1.7592194  -1.799442 -0.4650283  0.3306556
L1_AAACCTGTCCTTCAAT-1  0.4297156 -1.0515038  -3.547262  3.1270887 -1.8461041
                       harmony_45  harmony_46 harmony_47 harmony_48 harmony_49
L1_AAACCTGAGGGCTTCC-1 -0.32592691  1.97329800  0.2392751  0.7179419  0.3961060
L1_AAACCTGGTGCAGGTA-1  0.45834064 -1.42438576  0.2309867 -1.6689774  0.3301350
L1_AAACCTGGTTAAAGTG-1  0.08314617 -0.65869247 -1.5812362 -1.0555953  2.3985817
L1_AAACCTGTCAGGTAAA-1  0.22720771 -0.11615639 -0.8078383 -1.2600196  0.8054183
L1_AAACCTGTCCCTGACT-1  0.40463735  0.75701944 -0.3071644  1.4050215 -1.4059841
L1_AAACCTGTCCTTCAAT-1 -0.44516266 -0.08467537 -1.2774570 -0.3671591  0.0802340
                       harmony_50
L1_AAACCTGAGGGCTTCC-1 -1.43227758
L1_AAACCTGGTGCAGGTA-1 -0.34804648
L1_AAACCTGGTTAAAGTG-1  1.86972487
L1_AAACCTGTCAGGTAAA-1  0.01128852
L1_AAACCTGTCCCTGACT-1 -0.52700408
L1_AAACCTGTCCTTCAAT-1  2.67437059

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

# Run UMAP on Harmony embeddings
All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:16)

13:04:49 UMAP embedding parameters a = 0.9922 b = 1.112

13:04:49 Read 49372 rows and found 16 numeric columns

13:04:49 Using Annoy for neighbor search, n_neighbors = 30

13:04:49 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

**************************************************|
13:04:55 Writing NN index file to temp file /tmp/Rtmp8yjtoZ/file180ec3b102544
13:04:55 Searching Annoy index using 1 thread, search_k = 3000
13:05:14 Annoy recall = 100%
13:05:16 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
13:05:20 Initializing from normalized Laplacian + noise (using RSpectra)
13:05:23 Commencing optimization for 200 epochs, with 2099108 positive edges
13:06:28 Optimization finished

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

# Optionally, find neighbors and clusters (if you plan to do clustering analysis)
All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:16)

Computing nearest neighbor graph
Computing SNN

All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)  # Adjust resolution as needed

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

Number of nodes: 49372
Number of edges: 1435467

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8548
Number of communities: 15
Elapsed time: 14 seconds

1 singletons identified. 14 final clusters.

# Visualize UMAP
DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line", label = TRUE, pt.size = 0.5) +
    ggtitle("UMAP of Harmony-Integrated Data")

# Visualize UMAP with batch/cell line information
DimPlot(All_samples_Merged, reduction = "umap", group.by = "cell_line", label = TRUE, pt.size = 0.5) +
    ggtitle("UMAP - Colored by Cell Line (After Harmony Integration)")

# Visualize UMAP with clusters
DimPlot(All_samples_Merged, reduction = "umap", group.by = "seurat_clusters", label = TRUE, pt.size = 0.5) +
    ggtitle("UMAP - Clustered Data (After Harmony Integration)")

# Visualize specific cell types or other metadata
DimPlot(All_samples_Merged, reduction = "umap", group.by = "predicted.celltype.l2", label = TRUE, pt.size = 0.5) +
    ggtitle("UMAP - Cell Types After Harmony Integration")

12.Save the Seurat object as an Robj file

#save(All_samples_Merged, file = "../../../0-IMP-OBJECTS/All_Samples_Merged_with_10x_Azitmuth_Annotated_SCT_HPC_without_harmony_integration.robj")