Merged All samples with PBMC_10x and removed non CD4 T cells from Control 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-imp_Robj/SS_CD4_Tcells_Azimuth_Annotated_PBMC10x_excluding_nonCD4_cells_from_Control.robj")

All_samples_Merged <- seurat_filtered

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

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    6150    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 49515 

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

Dimensions of the RNA data layer: 36601 49515 

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

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 26195 by 49515

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

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

Second step: Get residuals using fitted parameters for 26195 genes

Computing corrected count matrix for 26195 genes

Calculating gene attributes

Wall clock passed: Time difference of 8.012294 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", "CC.Difference"), 
                                  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 26195 by 49515

Model formula is y ~ log_umi

Get Negative Binomial regression parameters per gene

Using 2000 genes, 5000 cells

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

Second step: Get residuals using fitted parameters for 26195 genes

Computing corrected count matrix for 26195 genes

Calculating gene attributes

Wall clock passed: Time difference of 6.716735 mins

Determine variable features

Regressing out percent.rb, percent.mt, CC.Difference

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, CC.Difference

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


# 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:  MALAT1, RPS27, PTPRC, SARAF, RPS4Y1, TCF7, RPL39, RPS29, RPL34, FYB1 
       ITM2B, GIMAP7, PNRC1, LEF1, BTG1, B2M, FOXP1, GIMAP5, IL7R, EVL 
       CD52, RIPOR2, LINC00861, SELL, YPEL3, CD3E, PIK3IP1, ITK, ZFP36L2, ANKRD12 
Negative:  PPIA, NPM1, RAN, NME2, GAPDH, PRDX1, PRELID1, RPS2, TUBA1B, HSP90AB1 
       HSPD1, RPLP0, TPI1, HMGA1, VDAC1, H2AFZ, SRM, TXN, ATP5F1B, CYC1 
       UBE2S, ENO1, PPP1R14B, NME1, RANBP1, CCT8, JPT1, SNRPD1, MIF, HNRNPAB 
PC_ 2 
Positive:  KIR3DL1, CST7, EPCAM, TRGV2, KIR2DL3, RPL27A, MATK, DAD1, KIR3DL2, XCL1 
       C1QBP, KLRC1, KIR2DL4, PFN1, MYO1E, RAB25, GZMM, CXCR3, CD7, NDUFA4 
       EIF4A1, ESYT2, ZBTB16, KLRK1, XCL2, KRT86, RCBTB2, SPINT2, PTPN6, TRGV4 
Negative:  RPL30, FAM107B, SEC11C, ELL2, CD74, LRRFIP1, MTHFD2, YBX3, VIM, RAD21 
       MT-ND3, SMAP2, CFLAR, TAP1, CD2, ANXA1, TMEM173, HDGFL3, IL2RA, CD70 
       AHNAK, HTATIP2, BATF3, MTLN, EEF1A1, CD58, C3orf14, RPL35A, RPL39, CCR7 
PC_ 3 
Positive:  TNFRSF4, C12orf75, HACD1, EGFL6, TIGIT, BACE2, ARPC2, SYT4, NET1, PXYLP1 
       LY6E, GRIA4, GGH, PARK7, MAP1B, UBE2D2, CCL17, KRT7, PON2, ADGRB3 
       PTP4A3, ACTN1, SCCPDH, MCTP2, CYBA, EEF1A2, PLEKHH2, SMIM3, MSC, HDGFL3 
Negative:  PAGE5, NDUFV2, RPL35A, RBPMS, CDKN2A, RPL22L1, CD74, STMN1, KIF2A, LMNA 
       TENM3, PSMB2, PSMB9, ANXA2, GPX4, FAM50B, PPP2R2B, GAPDH, ANXA5, IFI27L2 
       RPL11, VAMP5, FAM241A, NEURL1, PPBP, SH3KBP1, RPS3A, PLD1, HEBP2, TYMS 
PC_ 4 
Positive:  RPS4X, EGLN3, GAS5, WFDC1, KRT1, LINC02752, TTC29, TBX4, RPLP1, IL32 
       IFNGR1, AC069410.1, PLCB1, S100A11, TNS4, SP5, RPL13, FAM9C, SEMA4A, S100A4 
       IL4, NKG7, S100A6, LINC00469, VIM, PTGIS, HSPB1, CEBPD, VIPR2, NPTX1 
Negative:  EIF5A, HSPE1, ATP5MC3, ODC1, CHCHD10, PPID, MT-ND3, CYCS, PPBP, CYC1 
       RPL34, GCSH, FCER2, RPS4Y1, DNAJC12, FKBP4, FAM162A, GSTP1, CD7, FKBP11 
       HSPD1, RPS29, MTDH, RPL39, HSP90B1, TCF7, SET, PSMB6, SLC7A11-AS1, C1QBP 
PC_ 5 
Positive:  TMSB4X, TMSB10, LGALS1, S100A11, COTL1, S100A4, TP73, LSP1, IFITM2, S100A6 
       TMEM163, GPAT3, TAGLN2, HOXC9, LIME1, GPAT2, DUSP4, GAS2L1, EEF1A2, TNFRSF18 
       IFITM1, QPRT, MIIP, RBM38, LAPTM5, CEP135, PGAM1, CRIP1, PXYLP1, PRDX5 
Negative:  CCL17, MIR155HG, MAP4K4, MYO1D, RXFP1, LRBA, CA10, CFI, CA2, RUNX1 
       PRKCA, THY1, FRMD4A, AL590550.1, EZH2, IGHE, IMMP2L, HS3ST1, NFIB, LTA 
       SNTB1, RANBP17, MGST3, AKAP12, SLC35F3, AC100801.1, DOCK10, CCL5, EPB41L2, ONECUT2 

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

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

# 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

All_samples_Merged <- FindNeighbors(All_samples_Merged, 
                                dims = 1:22, 
                                verbose = FALSE)

# understanding resolution
All_samples_Merged <- FindClusters(All_samples_Merged, 
                                    resolution = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,0.8, 0.9, 1,1.2,1.5,2))

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9591
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: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

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

Number of nodes: 49515
Number of edges: 1697162

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8724
Number of communities: 38
Elapsed time: 13 seconds

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

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

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.3",
        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)

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

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

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

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

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.2", 
        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.7"

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   14   15 
6259 5929 5435 4945 3890 3347 3301 3237 2447 2217 1927 1919 1897  789  517  506 
  16   17   18   19 
 401  283  197   72 

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

                   
                       0    1    2    3    4    5    6    7    8    9   10   11
  B intermediate       0    0    0    3    0    1    1    0    0    0    0   55
  B memory             9    2    0   86  117    0   30    0    0    0    3   38
  B naive              0    0    0    0    0    0    0    0    0    0    0   41
  CD14 Mono            0    2    0    0    7    0    4    0    0    0    0   13
  CD4 CTL              0    0    0    0    0   11    0    0    0    0    0    2
  CD4 Naive            0    0    0    0    0  521    0    0 1479    0    0    9
  CD4 Proliferating 5450 5386 2852 4116 4064    0 3214 1530    5  930 1454    9
  CD4 TCM            878  523  267  588  482 4471  107 2323 1841  990   48  152
  CD4 TEM              0    0    0    0    0   61    0    1   22    0    0    0
  CD8 Proliferating    0    0    0    1    1    0    0    0    0    0    0    0
  CD8 TCM              0    0   16    0    0    0    0    1    0    0    0    0
  CD8 TEM              0    0    8    3    1    0    2    0    0    1    0    0
  cDC1                 0    0    0    6    0    0    2    0    0    0    0    0
  cDC2                 0    2    0    3   36    0    9    0    0    0    1    2
  dnT                  0    1    1    4    3    0    2    1    0    0    0    1
  HSPC                57    1    0  216  486    0  672    0    0    0  363   10
  ILC                  0    0    0    0    0    0    0    0    0    0    0    1
  NK                   0    0    0    0    0    0    0    0    0    0    0    1
  NK Proliferating     6   23 2785  263   34    0    9   24    0   15    2    0
  Treg                15    1    0   15    0    0    0    0    0    1    0   12
                   
                      12
  B intermediate       0
  B memory             0
  B naive              0
  CD14 Mono            0
  CD4 CTL              0
  CD4 Naive           33
  CD4 Proliferating    1
  CD4 TCM            160
  CD4 TEM              0
  CD8 Proliferating    0
  CD8 TCM              0
  CD8 TEM              0
  cDC1                 0
  cDC2                 0
  dnT                  2
  HSPC                 0
  ILC                  0
  NK                   0
  NK Proliferating     0
  Treg                 1

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 = T, label.box = T)

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

                   
                       0    1    2    3    4    5    6    7    8    9   10   11
  B intermediate       0    0    0    3    0    1    1    0    0    0    0   55
  B memory             9    2    0   86  117    0   30    0    0    0    3   38
  B naive              0    0    0    0    0    0    0    0    0    0    0   41
  CD14 Mono            0    2    0    0    7    0    4    0    0    0    0   13
  CD4 CTL              0    0    0    0    0   11    0    0    0    0    0    2
  CD4 Naive            0    0    0    0    0  521    0    0 1479    0    0    9
  CD4 Proliferating 5450 5386 2852 4116 4064    0 3214 1530    5  930 1454    9
  CD4 TCM            878  523  267  588  482 4471  107 2323 1841  990   48  152
  CD4 TEM              0    0    0    0    0   61    0    1   22    0    0    0
  CD8 Proliferating    0    0    0    1    1    0    0    0    0    0    0    0
  CD8 TCM              0    0   16    0    0    0    0    1    0    0    0    0
  CD8 TEM              0    0    8    3    1    0    2    0    0    1    0    0
  cDC1                 0    0    0    6    0    0    2    0    0    0    0    0
  cDC2                 0    2    0    3   36    0    9    0    0    0    1    2
  dnT                  0    1    1    4    3    0    2    1    0    0    0    1
  HSPC                57    1    0  216  486    0  672    0    0    0  363   10
  ILC                  0    0    0    0    0    0    0    0    0    0    0    1
  NK                   0    0    0    0    0    0    0    0    0    0    0    1
  NK Proliferating     6   23 2785  263   34    0    9   24    0   15    2    0
  Treg                15    1    0   15    0    0    0    0    0    1    0   12
                   
                      12
  B intermediate       0
  B memory             0
  B naive              0
  CD14 Mono            0
  CD4 CTL              0
  CD4 Naive           33
  CD4 Proliferating    1
  CD4 TCM            160
  CD4 TEM              0
  CD8 Proliferating    0
  CD8 TCM              0
  CD8 TEM              0
  cDC1                 0
  cDC2                 0
  dnT                  2
  HSPC                 0
  ILC                  0
  NK                   0
  NK Proliferating     0
  Treg                 1

11.Harmony Integration

# 
# # Load required libraries
# library(Seurat)
# library(harmony)
# library(ggplot2)
# 
# # Run Harmony, adjusting for batch effect using "cell_line" or another grouping variable
# All_samples_Merged <- RunHarmony(
#   object = All_samples_Merged,
#   group.by.vars = "cell_line",  # Replace with the metadata column specifying batch or cell line
#   dims.use = 1:22  # Use the same dimensions as PCA
# )
# 
# # Check results in harmony embeddings
# harmony_embeddings <- Embeddings(All_samples_Merged, reduction = "harmony")
# head(harmony_embeddings)
# 
# 
# # Run UMAP on Harmony embeddings
# All_samples_Merged <- RunUMAP(All_samples_Merged, reduction = "harmony", dims = 1:22)
# 
# # Optionally, find neighbors and clusters (if you plan to do clustering analysis)
# All_samples_Merged <- FindNeighbors(All_samples_Merged, reduction = "harmony", dims = 1:22)
# All_samples_Merged <- FindClusters(All_samples_Merged, resolution = 0.5)  # Adjust resolution as needed
# 
# # 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")