Celltype Prediction using different References-18-09-2025

load libraries————————————

1. Load and Subset Normal CD4 T Cells

ctrl <- readRDS("../0-Seurat_RDS_OBJECT_FINAL/Seurat_object_Final_changes/All_samples_Merged_with_STCAT_Annotation_final-5-09-2025.rds")

Idents(ctrl) <- "seurat_clusters"
table(ctrl@active.ident)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
6789 5275 4663 4661 4086 3634 3536 3409 3338 3273 3212 1675 1063  691 

2. Reference data

reference <- scPred::pbmc_1
reference

An object of class Seurat 
32838 features across 3500 samples within 1 assay 
Active assay: RNA (32838 features, 0 variable features)
 2 layers present: counts, data

Reference data

# Run SCTransform (normalization + variance stabilization)
reference <- SCTransform(reference, verbose = FALSE)

# Find variable features is done inside SCTransform, no need to run separately

# Run PCA on SCT assay
reference <- RunPCA(reference, assay = "SCT", verbose = FALSE, reduction.name = "pca")


# Run UMAP on PCA
reference <- RunUMAP(reference, dims = 1:30)

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

DimPlot(reference, group.by = "cell_type", label = TRUE, repel = TRUE) + NoAxes()

Run all steps of the analysis for the ctrl sample as well

DimPlot(ctrl, group.by = "seurat_clusters", label = TRUE, repel = TRUE) + NoAxes()

3. Label transfer

# Set default assay to SCT
DefaultAssay(reference) <- "SCT"

# 1. Make sure SCT is the default assay in your query (ctrl)
DefaultAssay(ctrl) <- "SCT"

# 2. (Optional but recommended) Run SCTransform if not already done
# If you already normalized with SCTransform, you can skip this.
# ctrl <- SCTransform(ctrl, verbose = FALSE)

# 3. Find anchors between reference and query
transfer.anchors <- FindTransferAnchors(
    reference = reference,
    query = ctrl,
    normalization.method = "SCT",   # MUST be SCT
    reference.reduction = "pca",   # Azimuth references use spca
    dims = 1:30
)

# 4. Transfer annotations (cell types, etc.)
predictions <- TransferData(
    anchorset = transfer.anchors,
    refdata = reference$cell_type,  # column from reference metadata
    dims = 1:30
)

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# 5. Add predictions back into your query object
ctrl <- AddMetaData(ctrl, metadata = predictions)

Dimplot-myumap

DimPlot(ctrl, group.by = "predicted.id", label = T, repel = T) + NoAxes()

barplot Distributions

ggplot(ctrl@meta.data, aes(x = seurat_clusters, fill = predicted.id)) +
    geom_bar() +
    theme_classic()

4. SinlgeR

  # Set default assay to SCT
  DefaultAssay(ctrl) <- "SCT"
  
  # Convert to SingleCellExperiment quietly
  sce <- as.SingleCellExperiment(ctrl)

Immune cell reference

 immune = celldex::DatabaseImmuneCellExpressionData()
singler.immune <- SingleR(test = sce, ref = immune, assay.type.test=1,
    labels = immune$label.fine)

head(singler.immune)

HPCA reference

hpca <- celldex::HumanPrimaryCellAtlasData()
singler.hpca <- SingleR(test = sce, ref = hpca, assay.type.test=1,
    labels = hpca$label.fine)

head(singler.hpca)

Compare results:

ctrl$singler.immune = singler.immune$pruned.labels
ctrl$singler.hpca = singler.hpca$pruned.labels


DimPlot(ctrl, group.by = c("singler.hpca", "singler.immune"), ncol = 2, label = T, label.box = T, repel = T)

Compare results:

ctrl$singler.immune = singler.immune$pruned.labels
ctrl$singler.hpca = singler.hpca$pruned.labels


DimPlot(ctrl, group.by = c("singler.hpca", "singler.immune"), ncol = 2, label = F, label.box = F, repel = T)

5. Azimuth

# options(future.globals.maxSize = 1e9)
# 
# library(SeuratData)
# 
# DefaultAssay(ctrl) <- "SCT"
# 
# ctrl <- RunAzimuth(ctrl, reference = "pbmcref")

Azimuth Dimplot

DimPlot(ctrl, group.by = "predicted.celltype.l1", label = T, repel = T) + NoAxes()

DimPlot(ctrl, group.by = "predicted.celltype.l2", label = T, repel = T) + NoAxes()

DimPlot(ctrl, group.by = "predicted.celltype.l3", label = T, repel = T) + NoAxes()

6. ProjectTils Annotation_CD4T_human_ref_v1

#reference atlas
 DimPlot(ref, label = T)

#Visualize projection
 plot.projection(ref, query.projected, linesize = 0.5, pointsize = 0.5)


 #Plot the predicted composition of the query in terms of reference T cell subtypes
 plot.statepred.composition(ref, query.projected, metric = "Percent")

# Now you can plot
DimPlot(ctrl, group.by = "functional.cluster", 
        reduction = "umap",
        label.size = 3,
        repel = TRUE,
        label = FALSE)

NA
NA

7. Compare results

crossTab(ctrl, "predicted.id", "singler.hpca")

crossTab(ctrl, "predicted.id", "singler.immune")

crossTab(ctrl, "singler.hpca", "singler.immune")

crossTab(ctrl, "predicted.id", "predicted.celltype.l1")

crossTab(ctrl, "predicted.celltype.l2", "functional.cluster")

crossTab(ctrl, "Prediction", "functional.cluster")

NA

Compare results

wrap_plots(
    DimPlot(ctrl, label = T, group.by = "predicted.id") + NoAxes() + ggtitle("LabelTransfer"),
    DimPlot(ctrl, label = F, group.by = "singler.hpca" ) + NoAxes() + ggtitle("SingleR HPCA"),
    DimPlot(ctrl, label = F, group.by = "singler.immune") + NoAxes() + ggtitle("SingleR Ref"),
    DimPlot(ctrl, label = T, group.by = "predicted.celltype.l1") + NoAxes() + ggtitle("Azimuth l1"),
    ncol = 2
)

Compare results

wrap_plots(
    DimPlot(ctrl, label = T, group.by = "predicted.id") + NoAxes() + ggtitle("LabelTransfer"),
    DimPlot(ctrl, label = F, group.by = "singler.hpca") + NoAxes() + ggtitle("SingleR HPCA"),
    DimPlot(ctrl, label = F, group.by = "singler.immune") + NoAxes() + ggtitle("SingleR Ref"),
    DimPlot(ctrl, label = F, group.by = "predicted.celltype.l2") + NoAxes() + ggtitle("LabelTransfer"),
    DimPlot(ctrl, label = T, group.by = "functional.cluster") + NoAxes() + ggtitle("ProjecTils CD4"),
    DimPlot(ctrl, label = F, group.by = "Prediction") + NoAxes() + ggtitle("STCAT"),
    ncol = 2
)

8. save all the annotation in object into RDS

# Save the Seurat object with all annotations
saveRDS(ctrl, file = "All_samples_Merged_with_all_annotations_ATLAS-18-09-2025.rds")