Advanced Usage

Introduction

This vignette covers advanced usage scenarios for scPharm, including:

Custom threshold calibration
Multi-drug analysis strategies
Combination therapy optimization
Integration with external tools
Performance tuning

library(ggplot2)
library(dplyr)

1. Threshold Calibration with Normal Tissue

Using scPharmGenNullDist

For accurate cell classification, calibrating thresholds using normal tissue samples is recommended:

library(scPharm)

# Load normal tissue reference
normal_seurat <- readRDS("normal_tissue.rds")

# Generate null distribution
null_dist <- scPharmGenNullDist(
  normal_seurat,
  cancer = "LUAD",
  drug = "Erlotinib",
  nmcs = 50,
  nfeatures = 200
)

# Extract calibrated thresholds
threshold_s <- null_dist$threshold.s
threshold_r <- null_dist$threshold.r

cat("Sensitive threshold:", threshold_s, "\n")
cat("Resistant threshold:", threshold_r, "\n")

Apply Calibrated Thresholds

# Load tumor sample
tumor_seurat <- readRDS("tumor_sample.rds")

# Run analysis with calibrated thresholds
result <- scPharmIdentify(
  tumor_seurat,
  type = "tissue",
  cancer = "LUAD",
  drug = "Erlotinib",
  threshold.s = threshold_s,
  threshold.r = threshold_r
)

Threshold Selection Strategy

Scenario	Recommendation
Strict classification	Higher `threshold.s`, lower `threshold.r`
Lenient classification	Lower `threshold.s`, higher `threshold.r`
Balanced	Use `scPharmGenNullDist()` defaults

2. Multi-Drug Analysis

Pan-Drug Screening

# Analyze all available drugs
result_all <- scPharmIdentify(
  seurat_obj,
  type = "cellline",
  cancer = "BRCA",
  drug = "all",
  cores = 8  # Parallel processing
)

# Get drug ranking
dr_scores <- scPharmDr(result_all)
print(head(dr_scores, 20))

Drug Class Analysis

# Load drug metadata
data(drug_info, package = "scPharm")

# Filter by drug class
chemo_drugs <- drug_info %>%
  filter(DRUG_TYPE == "chemotherapy") %>%
  pull(DRUG_NAME)

targeted_drugs <- drug_info %>%
  filter(DRUG_TYPE == "targeted") %>%
  pull(DRUG_NAME)

# Analyze specific drug classes
result_chemo <- scPharmIdentify(
  seurat_obj,
  cancer = "BRCA",
  drug = chemo_drugs
)

result_targeted <- scPharmIdentify(
  seurat_obj,
  cancer = "BRCA",
  drug = targeted_drugs
)

Multi-Cancer Analysis

# Analyze across multiple cancer contexts
cancers <- c("BRCA", "LUAD", "COREAD")

results <- lapply(cancers, function(cancer) {
  result <- scPharmIdentify(
    seurat_obj,
    type = "cellline",
    cancer = cancer,
    drug = "Paclitaxel"
  )
  
  # Extract NES values
  nes_col <- grep("scPharm_nes_", colnames(result@meta.data), value = TRUE)
  data.frame(
    cancer = cancer,
    mean_nes = mean(result@meta.data[[nes_col]], na.rm = TRUE),
    sd_nes = sd(result@meta.data[[nes_col]], na.rm = TRUE)
  )
})

# Combine results
multi_cancer_df <- do.call(rbind, results)
print(multi_cancer_df)

3. Combination Therapy Optimization

scPharmCombo Analysis

# Get drug ranking first
dr_scores <- scPharmDr(result)

# Identify combinations for top drugs
combos <- scPharmCombo(result, dr_scores, topN = 10)

# Examine results
for (drug_name in names(combos)) {
  cat("\n=== Drug:", drug_name, "===\n")
  print(head(combos[[drug_name]]$combination_scores))
}

Combination Selection Criteria

The combination score considers:

Complementarity: Do drugs target different resistant populations?
Coverage: What percentage of resistant cells become sensitive?
Synergy potential: Based on known drug interactions

Drug combination evaluation

4. Integration with Seurat Workflows

Pre-computed Embeddings

# Use existing Seurat clustering
seurat_obj <- FindVariableFeatures(seurat_obj)
seurat_obj <- ScaleData(seurat_obj)
seurat_obj <- RunPCA(seurat_obj)
seurat_obj <- FindNeighbors(seurat_obj)
seurat_obj <- FindClusters(seurat_obj)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:30)

# Run scPharm (MCA is computed independently)
result <- scPharmIdentify(seurat_obj, ...)

# Visualize with Seurat functions
DimPlot(result, group.by = "scPharm_label_Docetaxel")
FeaturePlot(result, features = "scPharm_nes_Docetaxel")

Cluster-Level Analysis

# Aggregate NES by cluster
cluster_summary <- result@meta.data %>%
  group_by(seurat_clusters) %>%
  summarise(
    n_cells = n(),
    mean_nes = mean(scPharm_nes_Docetaxel, na.rm = TRUE),
    pct_sensitive = mean(scPharm_label_Docetaxel == "sensitive") * 100,
    pct_resistant = mean(scPharm_label_Docetaxel == "resistant") * 100
  )

print(cluster_summary)

5. Performance Tuning

Memory Management

For large datasets, consider:

# Process in chunks
cell_chunks <- split(colnames(seurat_obj), 
                     ceiling(seq_along(colnames(seurat_obj)) / 5000))

results <- lapply(cell_chunks, function(cells) {
  subset_obj <- subset(seurat_obj, cells = cells)
  scPharmIdentify(subset_obj, ...)
})

# Merge results
final_result <- merge(results[[1]], results[-1])

Parallel Processing

# Detect available cores
n_cores <- parallel::detectCores() - 1

# Run with multiple cores
result <- scPharmIdentify(
  seurat_obj,
  drug = "all",
  cores = n_cores
)

Parameter Optimization

Parameter	Effect on Speed	Effect on Accuracy
`nmcs` ↑	Slower	Higher (to a point)
`nfeatures` ↑	Minimal	Higher specificity
`cores` ↑	Faster	None

6. Custom Drug Signatures

Using External Gene Sets

# Load custom drug sensitivity signatures
custom_signatures <- list(
  DrugA_sensitive = c("GENE1", "GENE2", "GENE3", ...),
  DrugA_resistant = c("GENE4", "GENE5", "GENE6", ...)
)

# Note: This requires modifying internal functions
# Contact maintainer for custom signature integration

7. Quality Control

Pre-analysis Checks

# Check gene coverage
data(bulkdata, package = "scPharm")
gene_overlap <- intersect(rownames(seurat_obj), rownames(bulkdata))
cat("Gene overlap:", length(gene_overlap), "/", nrow(bulkdata), "\n")

# Minimum recommended: 5000 genes
if (length(gene_overlap) < 5000) {
  warning("Low gene overlap may affect accuracy")
}

# Check cell numbers
if (ncol(seurat_obj) < 100) {
  warning("Small cell numbers may lead to unstable estimates")
}

Post-analysis Validation

# Check NES distribution
nes_values <- result@meta.data$scPharm_nes_Docetaxel

# Should be roughly normal/bimodal
hist(nes_values, breaks = 50, main = "NES Distribution")

# Check classification balance
table(result@meta.data$scPharm_label_Docetaxel)

8. Troubleshooting

Common Issues

Issue	Possible Cause	Solution
All cells classified as “other”	Thresholds too strict	Lower threshold_s, raise threshold_r
No tumor cells detected	CNV detection failed	Provide tumor.cells manually
Memory error	Dataset too large	Process in chunks
Slow performance	Too many drugs	Use parallel processing

Debug Mode

# Verbose output
options(scPharm.verbose = TRUE)

# Step-by-step execution
result <- scPharmIdentify(
  seurat_obj,
  type = "cellline",
  cancer = "BRCA",
  drug = "Docetaxel"
)

9. Exporting Results

To Data Frame

# Extract all scPharm columns
scpharm_cols <- grep("^scPharm_|^cell\\.label$", 
                     colnames(result@meta.data), value = TRUE)

export_df <- result@meta.data[, scpharm_cols]
export_df$cell_id <- rownames(export_df)

write.csv(export_df, "scpharm_results.csv", row.names = FALSE)

To Seurat Object

# Save annotated Seurat object
saveRDS(result, "annotated_seurat.rds")

Session Info

sessionInfo()
#> R version 4.5.3 (2026-03-11)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dplyr_1.2.1    ggplot2_4.0.3  rmarkdown_2.31
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.51         rlang_1.2.0       
#>  [5] xfun_0.57          otel_0.2.0         generics_0.1.4     S7_0.2.2          
#>  [9] jsonlite_2.0.0     labeling_0.4.3     glue_1.8.1         buildtools_1.0.0  
#> [13] htmltools_0.5.9    maketools_1.3.2    sys_3.4.3          sass_0.4.10       
#> [17] scales_1.4.0       grid_4.5.3         tibble_3.3.1       evaluate_1.0.5    
#> [21] jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5   
#> [25] compiler_4.5.3     RColorBrewer_1.1-3 pkgconfig_2.0.3    farver_2.1.2      
#> [29] digest_0.6.39      R6_2.6.1           tidyselect_1.2.1   pillar_1.11.1     
#> [33] magrittr_2.0.5     bslib_0.10.0       withr_3.0.2        tools_4.5.3       
#> [37] gtable_0.3.6       cachem_1.1.0