How to create .h5ad for STCAT

1. load libraries

2. Load Data into Seurat

ss_Bor_2023 <- readRDS("../Borch_Combined/data/ss_Borcherding_Malignant_6_Normal_1_Integrated_object.rds")

3. Export to .h5ad (AnnData format)

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


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

# Step 3: Prepare Metadata (Optional, Recommended)
cell_metadata <- ss_Bor_2023[[]]  # 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)
writeH5AD(sce, file = "STCAT/ss_Bor_2023_raw_counts.h5ad")

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_Bor_2023_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_Bor_2023_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_Annotation_Borcherding/STCAT_annotations.csv", row.names = 1)

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

[1] "SS_P1_AAATGCCTCAACGGCC-1" "SS_P1_AATCGGTGTTCCAACA-1" "SS_P1_ACACCCTTCGCGTTTC-1" "SS_P1_ACAGCCGAGCTCCTCT-1"
[5] "SS_P1_ACTTACTAGCTAAGAT-1" "SS_P1_AGGCCACGTACCGTAT-1"

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

view(ss_Bor_2023@meta.data)

6. Read and merge with Seurat (STCAT annotation)

DimPlot(ss_Bor_2023_annotated, group.by = "Disease_state", label = TRUE, repel = T)

DimPlot(ss_Bor_2023_annotated, group.by = "orig.ident", label = TRUE, repel = T)

DimPlot(ss_Bor_2023_annotated, group.by = "integrated_snn_res.0.5", label = TRUE, repel = T)

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

NA
NA
NA
NA

7. Read and merge with Seurat (STCAT annotation)

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

 CD4 proliferation             CD4 Tc            CD4 Tcm            CD4 Tem          CD4 Temra            CD4 Tex 
               376               2353              21126               1320                103               1832 
           CD4 Tfh            CD4 Th1           CD4 Th17           CD4 Tisg             CD4 Tn    CD4 Tn adhesion 
               224                 63                725               1236               4420               3109 
          CD4 Treg CD4 Treg activated            CD4 Trm CD4 Trm cell-death CD4 Trm naive-like           CD4 Tstr 
              1729                132               1874               3252                621                434 
     Heterogeneous             None T                Tgd 
               128               5678                360 

table(ss_Bor_2023$orig.ident)         # should show L1–L7

Control   SS_P1   SS_P2   SS_P3   SS_P4   SS_P5   SS_P6 
   4437    3443   34179    3084     849    1582    3521 

# Contingency table of STCAT Prediction × Cell Line
table_annotation_orig_ident <- table(ss_Bor_2023$STCAT_Prediction, ss_Bor_2023$orig.ident)
table_annotation_cluster <- table(ss_Bor_2023$STCAT_Prediction, ss_Bor_2023$integrated_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_orig_ident, margin = 1)

# Heatmap with raw counts as labels
pheatmap(prop_table,
         display_numbers = table_annotation_orig_ident,
         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")