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 cell cycle 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 7.827107 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", "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 7.115858 mins

Determine variable features

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

# 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

# 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: 1625469

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

Number of nodes: 49372
Number of edges: 1625469

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9465
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: 1625469

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9377
Number of communities: 17
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: 1625469

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9294
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: 49372
Number of edges: 1625469

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9212
Number of communities: 20
Elapsed time: 14 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 
6404 5975 5930 5863 5278 5076 4938 4091 3341 1748  397  196   76   59 

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    0    2    1    0    0    1    0
  B memory             8    6    0    1  118    0   78   31    0    1    9    0
  CD14 Mono            0    0    0    0    7    0    0    4    0    0    0    0
  CD4 CTL              0    0    0    0    0   12    0    0    0    0    0    0
  CD4 Naive            0    8    0    0    0  521    0    0 1479    0    0   33
  CD4 Proliferating 5448 2473 2852 5320 4109    0 3891 3247    6 1348  249    0
  CD4 TCM            872 3415  268  514  484 4481  523  108 1835   43   64  161
  CD4 TEM              0    1    0    0    0   62    0    0   21    0    0    0
  CD8 Proliferating    0    0    0    0    1    0    1    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    0    3    2    0    0    0    0
  cDC1                 0    0    0    0    0    0    4    2    0    0    2    0
  cDC2                 0    0    0    2   35    0    3   10    0    1    0    0
  dnT                  0    3    1    1    3    0    1    2    0    0    3    1
  HSPC                57   10    0    1  486    0  202  674    0  354   21    0
  NK Proliferating     4   40 2785   23   34    0  228   10    0    1   36    0
  Treg                15   14    0    1    0    0    2    0    0    0   12    1
                   
                      12   13
  B intermediate       0    0
  B memory             0    0
  CD14 Mono            0    1
  CD4 CTL              0    1
  CD4 Naive            0    1
  CD4 Proliferating   68    0
  CD4 TCM              8   54
  CD4 TEM              0    0
  CD8 Proliferating    0    0
  CD8 TCM              0    0
  CD8 TEM              0    0
  cDC1                 0    0
  cDC2                 0    2
  dnT                  0    0
  HSPC                 0    0
  NK Proliferating     0    0
  Treg                 0    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    0    2    1    0    0    1    0
  B memory             8    6    0    1  118    0   78   31    0    1    9    0
  CD14 Mono            0    0    0    0    7    0    0    4    0    0    0    0
  CD4 CTL              0    0    0    0    0   12    0    0    0    0    0    0
  CD4 Naive            0    8    0    0    0  521    0    0 1479    0    0   33
  CD4 Proliferating 5448 2473 2852 5320 4109    0 3891 3247    6 1348  249    0
  CD4 TCM            872 3415  268  514  484 4481  523  108 1835   43   64  161
  CD4 TEM              0    1    0    0    0   62    0    0   21    0    0    0
  CD8 Proliferating    0    0    0    0    1    0    1    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    0    3    2    0    0    0    0
  cDC1                 0    0    0    0    0    0    4    2    0    0    2    0
  cDC2                 0    0    0    2   35    0    3   10    0    1    0    0
  dnT                  0    3    1    1    3    0    1    2    0    0    3    1
  HSPC                57   10    0    1  486    0  202  674    0  354   21    0
  NK Proliferating     4   40 2785   23   34    0  228   10    0    1   36    0
  Treg                15   14    0    1    0    0    2    0    0    0   12    1
                   
                      12   13
  B intermediate       0    0
  B memory             0    0
  CD14 Mono            0    1
  CD4 CTL              0    1
  CD4 Naive            0    1
  CD4 Proliferating   68    0
  CD4 TCM              8   54
  CD4 TEM              0    0
  CD8 Proliferating    0    0
  CD8 TCM              0    0
  CD8 TEM              0    0
  cDC1                 0    0
  cDC2                 0    2
  dnT                  0    0
  HSPC                 0    0
  NK Proliferating     0    0
  Treg                 0    0

Save the Seurat object as an Robj file

save(All_samples_Merged, file = "../0-R_Objects/CD4Tcells_SCTnormalized_regress_cellcycle_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
  theta =0.5,
  assay.use="SCT")

Transposing data matrix

Initializing state using k-means centroids initialization

Warning: Quick-TRANSfer stage steps exceeded maximum (= 2468600)

Harmony 1/10

Harmony 2/10

Harmony converged after 2 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  harmony_5
L1_AAACCTGAGGGCTTCC-1   2.8110708  0.5260470 -2.029509  -5.414869  0.4810349
L1_AAACCTGGTGCAGGTA-1 -10.8277517 -1.3080442  2.621653 -12.684867 -6.6769900
L1_AAACCTGGTTAAAGTG-1 -12.5325103 -2.3480316 -2.091807  -7.248564 -2.5818361
L1_AAACCTGTCAGGTAAA-1   0.7168548 -2.0496518 -7.199748  -1.961662  2.3792218
L1_AAACCTGTCCCTGACT-1   1.8990115 -0.8495511 -1.194518  -2.758270 -1.3860372
L1_AAACCTGTCCTTCAAT-1 -14.4141663 -1.6009173  5.471255  -6.792593 -7.5389136
                      harmony_6  harmony_7  harmony_8 harmony_9 harmony_10
L1_AAACCTGAGGGCTTCC-1  1.061357 -2.5527498  3.1985310 -2.746724  -6.830058
L1_AAACCTGGTGCAGGTA-1  9.410060 -2.7334282 -1.3222509  1.055235  -1.004011
L1_AAACCTGGTTAAAGTG-1  5.746505 -1.7399981  4.0411710  4.279450   3.112614
L1_AAACCTGTCAGGTAAA-1  1.533169  0.5831356  0.3504292  2.237259   3.447279
L1_AAACCTGTCCCTGACT-1  1.237426 -2.0802556  1.9108124  0.174732  -3.821618
L1_AAACCTGTCCTTCAAT-1  5.175538 -5.4316561  4.4416909  2.553290  -1.915657
                      harmony_11  harmony_12  harmony_13 harmony_14 harmony_15
L1_AAACCTGAGGGCTTCC-1 -0.9714769 -0.05527739  0.15581513   1.732362 -0.2018541
L1_AAACCTGGTGCAGGTA-1  3.0184137 -1.61621018 -0.04082989  -2.164035 -1.8986730
L1_AAACCTGGTTAAAGTG-1 -1.7221516  2.55422857 -2.60542844  -5.640573 -4.7569498
L1_AAACCTGTCAGGTAAA-1 -0.2740579  2.57582176 -1.07542398  -2.386170 -1.0541965
L1_AAACCTGTCCCTGACT-1 -0.3247629  0.54034581 -0.04814952   1.265659  2.0557533
L1_AAACCTGTCCTTCAAT-1 -0.1971245 -1.58663476 -0.52186800  -1.627686 -0.3703728
                      harmony_16 harmony_17 harmony_18 harmony_19 harmony_20
L1_AAACCTGAGGGCTTCC-1  1.4182175  0.8927475 -0.2090148  2.7118999  -1.618971
L1_AAACCTGGTGCAGGTA-1 -1.1801555  0.2123179  2.9383929 -1.9859795  -6.354939
L1_AAACCTGGTTAAAGTG-1 -0.4801443 -3.3607603  1.4361882 -0.9222678  -1.423420
L1_AAACCTGTCAGGTAAA-1 -1.0764376  1.3806085  0.9536258  0.5033541   2.127277
L1_AAACCTGTCCCTGACT-1  4.0194691 -4.2731920 -1.9344298  3.3722244  -0.404554
L1_AAACCTGTCCTTCAAT-1  0.1813054 -1.3051187  1.8163951 -0.7331759  -1.153906
                      harmony_21 harmony_22 harmony_23 harmony_24 harmony_25
L1_AAACCTGAGGGCTTCC-1   3.840658   5.463864 -4.9235944 -0.6539436  2.3889750
L1_AAACCTGGTGCAGGTA-1   4.356782   4.226991 -0.7697845  0.9406667  0.9692749
L1_AAACCTGGTTAAAGTG-1   3.676776   2.112633 -0.9075745  2.1831911  0.9240647
L1_AAACCTGTCAGGTAAA-1  -2.483335  -1.448657  1.9807517 -0.6691265 -0.8103514
L1_AAACCTGTCCCTGACT-1   5.803254   4.323493 -2.1814562  0.4847289 -0.2041591
L1_AAACCTGTCCTTCAAT-1   2.230887   4.614193 -2.5862105 -2.5612540  0.3535346
                       harmony_26 harmony_27 harmony_28 harmony_29 harmony_30
L1_AAACCTGAGGGCTTCC-1 -0.49859576  1.5008159  0.4244375  0.8604441 -1.5174028
L1_AAACCTGGTGCAGGTA-1 -0.07987602  0.1661470  0.8491917 -1.2835925 -0.3087157
L1_AAACCTGGTTAAAGTG-1 -1.84695160  1.7279675 -0.4244429  2.3041313 -2.0890832
L1_AAACCTGTCAGGTAAA-1 -1.59549939  0.2004768  0.8547986 -2.0092601 -1.2147589
L1_AAACCTGTCCCTGACT-1  0.97805498  0.1558526 -0.8726602  1.1839138 -0.6738418
L1_AAACCTGTCCTTCAAT-1 -0.22986679  2.2838265  0.3776235 -1.9399405 -2.5364713
                      harmony_31 harmony_32 harmony_33 harmony_34 harmony_35
L1_AAACCTGAGGGCTTCC-1  2.4931168  1.2598178  1.4299188  0.9984621  -3.274602
L1_AAACCTGGTGCAGGTA-1 -2.0953837 -3.5361026 -2.4394641 -1.8543248   2.233153
L1_AAACCTGGTTAAAGTG-1  4.3203859  0.8505935  0.6683941  1.0755620  -1.239794
L1_AAACCTGTCAGGTAAA-1 -1.3730414 -1.5866399 -1.4500221  0.3013226  -1.476863
L1_AAACCTGTCCCTGACT-1  2.3268213  2.1669031  2.5589161  1.5222361  -1.881328
L1_AAACCTGTCCTTCAAT-1  0.9452486 -1.7511897 -2.2235795  1.5247648  -3.331400
                       harmony_36 harmony_37 harmony_38  harmony_39  harmony_40
L1_AAACCTGAGGGCTTCC-1  0.73833256 0.90280566  0.7604341 -0.03901279  0.21012290
L1_AAACCTGGTGCAGGTA-1  1.07019382 0.07081572 -0.3472956 -0.03092420 -2.88190834
L1_AAACCTGGTTAAAGTG-1  0.07938292 0.78780908  1.2989845 -1.76037099  0.01909956
L1_AAACCTGTCAGGTAAA-1 -0.54561715 0.67121861 -0.6762170 -1.48964267  0.28825101
L1_AAACCTGTCCCTGACT-1  1.87907707 1.92788476  1.4872688  0.31696953  2.08077128
L1_AAACCTGTCCTTCAAT-1  3.89226889 1.57740339  1.2876424 -0.41865757 -1.31192084
                      harmony_41 harmony_42 harmony_43  harmony_44 harmony_45
L1_AAACCTGAGGGCTTCC-1   0.489528 -1.0329314  0.2194132  0.10138125 -1.7691260
L1_AAACCTGGTGCAGGTA-1  -1.219079  1.3811536  1.4949549  1.02236295  0.2819950
L1_AAACCTGGTTAAAGTG-1   2.228459 -2.4615098  1.7832945 -0.05588167  2.2745851
L1_AAACCTGTCAGGTAAA-1  -1.023605 -2.1711387 -1.0596324  0.79648989  0.8016935
L1_AAACCTGTCCCTGACT-1   2.483687  0.0620254 -0.6379408  0.58823324 -1.2029052
L1_AAACCTGTCCTTCAAT-1   3.225511 -3.2063560 -0.2389498  0.87613710 -1.5828924
                      harmony_46 harmony_47 harmony_48 harmony_49  harmony_50
L1_AAACCTGAGGGCTTCC-1 -0.9828127 -0.5140391 -0.6583532  0.9544030 -0.02661797
L1_AAACCTGGTGCAGGTA-1  0.7282084  1.3802107  2.1433726  1.0424374 -1.16203040
L1_AAACCTGGTTAAAGTG-1 -2.7526590 -2.1933917  2.1782184 -0.7480607 -2.21273992
L1_AAACCTGTCAGGTAAA-1  0.5060130  0.8837644  1.9178650  0.2556520 -1.34849293
L1_AAACCTGTCCCTGACT-1 -0.7984049  0.8086429 -1.3205108 -0.4491508  0.13190215
L1_AAACCTGTCCTTCAAT-1 -2.3086164  0.3514620  1.4098153 -0.8431389 -2.71995300

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

19:54:26 UMAP embedding parameters a = 0.9922 b = 1.112

19:54:26 Read 49372 rows and found 16 numeric columns

19:54:26 Using Annoy for neighbor search, n_neighbors = 30

19:54:26 Building Annoy index with metric = cosine, n_trees = 50

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

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

**************************************************|
19:54:31 Writing NN index file to temp file /tmp/RtmpvgbiFo/filecd707e69cc04
19:54:31 Searching Annoy index using 1 thread, search_k = 3000
19:54:48 Annoy recall = 100%
19:54:50 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
19:54:55 Initializing from normalized Laplacian + noise (using RSpectra)
19:54:58 Commencing optimization for 200 epochs, with 2055060 positive edges
19:56:00 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, reduction = "harmony", resolution = 0.5)  # Adjust resolution as needed

Warning: The following arguments are not used: reduction

Warning: The following arguments are not used: reduction

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

Number of nodes: 49372
Number of edges: 1554209

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

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