Celltype Prediction using Herrera data

1 load libraries

2 Load and Subset Herrera Data

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

   2    3    0   12    6   17    5    4    1   15   14   18    9   13    7    8   11   16   10 
3241 1006  185  124  404   44  223 2278   92   48  183   10   84  150  230  923  727    4  829 

2.1 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

2.2 Reference data analysis

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

Registered S3 method overwritten by 'rmarkdown':
  method         from
  print.paged_df     

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

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

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

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

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

2.4 Label transfer

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

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

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

2.5 Dimplot-myumap

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

2.6 barplot Distributions

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

3 SinlgeR

 # make sure you are on Seurat >= 5
DefaultAssay(ctrl) <- "RNA"

# merge all RNA layers into one
ctrl <- JoinLayers(ctrl, assay = "RNA")

# now convert to SCE
sce <- as.SingleCellExperiment(ctrl, assay = "RNA")

3.1 Immune cell reference

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

head(singler.immune)

DataFrame with 6 rows and 4 columns
                                                          scores                 labels delta.next          pruned.labels
                                                        <matrix>            <character>  <numeric>            <character>
SS_SS1_Blood_H20_AAACGGGTCTGTACGA 0.126734:0.123792:0.127111:...     T cells, CD4+, Th2 0.01020391     T cells, CD4+, Th2
SS_SS1_Blood_H20_AAAGCAAAGTCCCACG 0.146714:0.150715:0.153418:...     T cells, CD4+, Th2 0.00787036     T cells, CD4+, Th2
SS_SS1_Blood_H20_AAAGTAGAGTTTCCTT 0.147983:0.135853:0.135882:...     T cells, CD4+, Th2 0.06289714     T cells, CD4+, Th2
SS_SS1_Blood_H20_AAAGTAGTCCACGACG 0.127739:0.115666:0.116248:...     T cells, CD4+, Th2 0.00174197     T cells, CD4+, Th2
SS_SS1_Blood_H20_AAATGCCCACCCATGG 0.140546:0.143729:0.155259:... T cells, CD4+, memor.. 0.00138212 T cells, CD4+, memor..
SS_SS1_Blood_H20_AACACGTAGCCCAATT 0.137978:0.121588:0.138889:...     T cells, CD4+, Th2 0.00429020     T cells, CD4+, Th2

3.2 HPCA reference

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

head(singler.hpca)

DataFrame with 6 rows and 4 columns
                                                          scores                 labels delta.next          pruned.labels
                                                        <matrix>            <character>  <numeric>            <character>
SS_SS1_Blood_H20_AAACGGGTCTGTACGA 0.158345:0.252571:0.222786:... T_cell:CD4+_central_..  0.1610613 T_cell:CD4+_central_..
SS_SS1_Blood_H20_AAAGCAAAGTCCCACG 0.167377:0.305519:0.274721:... T_cell:CD4+_central_..  0.1069279 T_cell:CD4+_central_..
SS_SS1_Blood_H20_AAAGTAGAGTTTCCTT 0.182532:0.312291:0.278981:... T_cell:CD4+_central_..  0.0303539 T_cell:CD4+_central_..
SS_SS1_Blood_H20_AAAGTAGTCCACGACG 0.168245:0.274507:0.247855:... T_cell:CD4+_effector..  0.0999581 T_cell:CD4+_effector..
SS_SS1_Blood_H20_AAATGCCCACCCATGG 0.182977:0.305724:0.273373:... T_cell:CD4+_central_..  0.0980566 T_cell:CD4+_central_..
SS_SS1_Blood_H20_AACACGTAGCCCAATT 0.175543:0.282729:0.252079:... T_cell:CD4+_effector..  0.0399778 T_cell:CD4+_effector..

3.3 Compare results:

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

3.4 Compare results:

DimPlot(ctrl, group.by = "singler.hpca", label = T, repel = T) + NoAxes()

DimPlot(ctrl, group.by = "singler.immune", label = T, repel = T) + NoAxes()

NA
NA

3.5 Compare results:

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

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

4 Azimuth

options(future.globals.maxSize = 1e9)

library(SeuratData)

DefaultAssay(ctrl) <- "SCT"

ctrl <- RunAzimuth(ctrl, reference = "pbmcref")

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

  |                                                  | 0 % ~calculating  
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  

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

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

5 ProjectTils Annotation_CD4T_human_ref_v1

6 Compare results from all Annotation methods

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

NA

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

6.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")
  
)

7 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-Herrera-03-02-2026.rds")