How to create .h5ad for STCAT of Herrera data

1. load libraries

2. Load Data into Seurat

ss_Herrera <- readRDS("Herrara_All_samples.rds")

3. Export to .h5ad (AnnData format)

# Step 1: Rejoin Layers to Get Unified Raw Counts
ss_Herrera[["RNA"]] <- JoinLayers(ss_Herrera[["RNA"]])
raw_counts <- ss_Herrera[["RNA"]]$counts


# Step 2: Transpose to Cells × Genes
# Transpose: genes × cells → cells × genes
counts_transposed <- t(as.matrix(raw_counts))

Warning: sparse->dense coercion: allocating vector of size 7.6 GiB

# Step 3: Prepare Metadata (Optional, Recommended)
cell_metadata <- ss_Herrera[[]]  # All cell metadata

# Step 4: Convert to SingleCellExperiment
library(SingleCellExperiment)
sce <- SingleCellExperiment(
  assays = list(counts = raw_counts),       # Keep raw_counts (genes × cells!)
  colData = cell_metadata                   # Attach cell metadata
)

#Step 5: Export to .h5ad
library(zellkonverter)

Registered S3 method overwritten by 'zellkonverter':
  method                                             from      
  py_to_r.pandas.core.arrays.categorical.Categorical reticulate

writeH5AD(sce, file = "STCAT/ss_Herrera_raw_counts.h5ad")

ℹ Using the 'counts' assay as the X matrix

4. Python Script to run (make sure you have right directory)

# import warnings
# warnings.filterwarnings("ignore")
# import STCAT
# import scanpy as sc
# 
# 
# # Load your .h5ad data
# adata = sc.read("/home/bioinfo/2025_NewHarmony_Integrated_Files/Year3_Analysis/Biomarkers_Validation_with_Public_Data/Borch_Combined/STCAT/STCAT_Annotation_Borcherding/ss_Herrera_raw_counts.h5ad")
# print(adata.var_names[:10])
# 
# 
# # Run STCAT
# results = STCAT.STCAT(adata)
# 
# # Save results
# results.write("/home/bioinfo/2025_NewHarmony_Integrated_Files/Year3_Analysis/Biomarkers_Validation_with_Public_Data/Borch_Combined/STCAT/STCAT_Annotation_Borcherding/STCAT_Annotation_ss_Herrera_raw_counts.h5ad")
# 
# 
# print(adata)
# 

5. Read and merge with Seurat (STCAT annotation)

# Load Seurat and read your annotated table
library(Seurat)
library(dplyr)

# Read your annotation table
annots <- read.csv("STCAT/STCAT_annotations.csv", row.names = 1)

# Confirm the rownames are cell barcodes
head(rownames(annots))  # e.g., AAACCTGGT...

[1] "SS1_Blood_H20_AAGGTTCCATATGAGA" "SS1_Blood_H20_CACTCCAGTGCACCAC" "SS1_Blood_H20_AACTCCCTCCGAACGC"
[4] "SS1_Blood_H20_AACTCTTAGACCGGAT" "SS1_Blood_H20_AAGGCAGAGTCAAGCG" "SS1_Blood_H20_AAGTCTGGTCTAAAGA"

# Match to your Seurat object (assuming your Seurat object is named `ss_Herrera`)
ss_Herrera$STCAT_Prediction <- annots[Cells(ss_Herrera), "Prediction"]
ss_Herrera$STCAT_Cluster <- annots[Cells(ss_Herrera), "Cluster"]
ss_Herrera$STCAT_Uncertainty <- annots[Cells(ss_Herrera), "Uncertainty.score"]

view(ss_Herrera@meta.data)

6. Read and merge with Seurat (STCAT annotation)

DimPlot(ss_Herrera, group.by = "tissue", label = TRUE, repel = T)

DimPlot(ss_Herrera, group.by = "sample_id", label = TRUE, repel = T)

DimPlot(ss_Herrera, group.by = "RNA_snn_res.0.5", label = TRUE, repel = T)

DimPlot(ss_Herrera, group.by = "STCAT_Prediction", label = TRUE, repel = T)

NA
NA
NA
NA

7. Read and merge with Seurat (STCAT annotation)

table(ss_Herrera$STCAT_Prediction)  # should show all predicted labels

      CD4 activated   CD4 proliferation              CD4 Tc             CD4 Tcm             CD4 Tem           CD4 Temra 
                 40                 108                 361                4935                 272                  62 
            CD4 Tex             CD4 Tfh            CD4 Th17            CD4 Tisg              CD4 Tn     CD4 Tn adhesion 
                906                  31                 137                 702                3927                 163 
           CD4 Treg CD4 Treg naive-like             CD4 Trm  CD4 Trm cell-death   CD8 proliferation      CD8 senescence 
                228                 142                 615                1134                 134                   4 
             CD8 Tc             CD8 Tcm             CD8 Tem           CD8 Temra             CD8 Tex              CD8 Tn 
                743                  19                 669                1264                   2                1381 
            CD8 Trm  CD8 Trm naive-like            CD8 Tstr     Double Negative                MAIT              None T 
                142                 476                  26                  32                 385                9622 
        other CD8 T                 Tgd 
                  9                1531 

table(ss_Herrera$sample_id)         # should show L1–L7

Healthy_Blood  Healthy_Skin     MF1_Blood      MF1_Skin     SS1_Blood      SS1_Skin     SS2_Blood      SS2_Skin 
         4386            33          4880           293          3631           374          6078           244 
    SS3_Blood      SS3_Skin     SS4_Blood      SS4_Skin     SS5_Blood     SS6_Blood 
         2109           272          2072           320          2832          2821 

table(ss_Herrera$tissue) 

Blood  Skin 
28809  1536 

# Contingency table 
table_annotation_tissue <- table(ss_Herrera$STCAT_Prediction, ss_Herrera$tissue)
table_annotation_sample_id <- table(ss_Herrera$STCAT_Prediction, ss_Herrera$sample_id)
table_annotation_cluster <- table(ss_Herrera$STCAT_Prediction, ss_Herrera$RNA_snn_res.0.5)

library(pheatmap)

# Proportion table (color scale)
prop_table <- prop.table(table_annotation_cluster, margin = 1)

# Heatmap with raw counts as labels
pheatmap(prop_table,
         display_numbers = table_annotation_cluster,
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         color = colorRampPalette(c("lightyellow", "yellow", "lightgreen", "white"))(100),
         fontsize_row = 10,
         fontsize_col = 10,
         main = "STCAT Prediction vs Clusters")

library(pheatmap)

# Proportion table (color scale)
prop_table <- prop.table(table_annotation_sample_id, margin = 1)

# Heatmap with raw counts as labels
pheatmap(prop_table,
         display_numbers = table_annotation_sample_id,
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         color = colorRampPalette(c("lightyellow", "yellow", "lightgreen", "white"))(100),
         fontsize_row = 10,
         fontsize_col = 10,
         main = "STCAT Prediction vs Patients")

# Proportion table (color scale)
prop_table <- prop.table(table_annotation_tissue, margin = 1)

# Heatmap with raw counts as labels
pheatmap(prop_table,
         display_numbers = table_annotation_tissue,
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         color = colorRampPalette(c("lightyellow", "yellow", "lightgreen", "white"))(100),
         fontsize_row = 10,
         fontsize_col = 10,
         main = "STCAT Prediction vs Patients")

Save the Seurat object as an RDS

saveRDS(ss_Herrera, file = "STCAT/All_Patients_Integrated_ss_Herrera_annotated.RDS")