--- title: "4. Apply pdactrace to your own multi-omics cohort" author: "Jibeom Ko" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{4. Apply pdactrace to your own multi-omics cohort} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") suppressPackageStartupMessages({ library(pdactrace) library(data.table) }) set.seed(42) ``` # Goal The bundled v0.3.0 reference atlas covers four PDAC RNA cohorts, two PDAC tissue-protein cohorts, ten scRNA cohorts, and three serum-DIA cohorts. v0.4.0 adds an additive, backwards-compatible API so a new user can run the same frozen 3-axis + 2-gate audit framework on their own multi-omics cohort, without retraining and without modifying the bundled atlas. This vignette walks through the full five-step workflow on a synthetic example so the framework's input/output schema is concrete. # Step 1 — Load the bundled toy cohort The package ships with a tiny synthetic cohort (`toy_counts` / `toy_protein` / `toy_coldata`) so the user-cohort workflow can be exercised end-to-end without an external download. 50 genes / proteins x 24 samples balanced across the four canonical stages and two synthetic cohorts; the first 10 RNA rows and first 8 protein rows carry an injected `Early_Burst_Up` signal. ```{r toy-rna} data(toy_counts) data(toy_coldata) data(toy_protein) sample_meta <- toy_coldata counts <- toy_counts prot_intensity <- toy_protein up_genes <- rownames(counts)[1:10] # injected stage-progressive RNA dim(counts); table(sample_meta$stage, sample_meta$cohort) ``` # Step 2 — Fit the stage-aware RNA model ```{r rna-fit, eval = requireNamespace("DESeq2", quietly = TRUE)} rna_fit <- fit_stage_de(counts, stage = sample_meta$stage, cohort = sample_meta$cohort) rna_pat <- classify_trajectory(rna_fit, rho_cutoff = 0.85) head(rna_pat[!is.na(rna_pattern), .(gene_symbol, rna_pattern, rna_pattern_rho, lrt_padj)]) ``` # Step 3 — Fit the protein limma model ```{r prot-fit, eval = requireNamespace("limma", quietly = TRUE)} prot_fit <- fit_stage_de_protein(prot_intensity, stage = sample_meta$stage, cohort = sample_meta$cohort) # `sig_only = FALSE` so the toy cohort (where the protein signal is # small relative to noise) still surfaces best-template calls. prot_pat <- classify_protein_trajectory(prot_fit, rho_cutoff = 0.85, sig_only = FALSE) if ("prot_pattern" %in% names(prot_pat)) { head(prot_pat[!is.na(prot_pattern), .(gene_symbol, prot_pattern, prot_pattern_rho)]) } else { message("No protein patterns surfaced on the toy cohort.") } ``` # Step 4 — Assemble the user evidence table ```{r assemble, eval = exists("rna_pat") && exists("prot_pat")} user_ev <- assemble_user_evidence( rna_fit = rna_pat, prot_fit = prot_pat, signal_peptide = up_genes) str(user_ev, max.level = 1) ``` # Step 5 — Run the audit framework on user evidence ```{r audit, eval = exists("user_ev")} user_audit <- compute_audit_score(evidence = user_ev) table(user_audit$audit_class) head(user_audit[order(-audit_score)], 10) ``` The same `evaluate_anchor_enrichment()` interface accepts an optional `anchors` table so a user can evaluate their own anchor list: ```{r enrich, eval = exists("user_ev")} user_anchors <- data.table( gene = up_genes[1:5], include_primary_eval = TRUE, include_secondary_eval = TRUE) evaluate_anchor_enrichment(top_n = 20, tier = "secondary", evidence = user_ev, anchors = user_anchors) ``` # Step 6 — Bioconductor-native input via `SummarizedExperiment` (v0.99.2) Steps 2 and 3 above show the matrix interface. `fit_stage_de()` and `fit_stage_de_protein()` are also S4 generics with a `SummarizedExperiment` method, so the same workflow can be driven from a Bioconductor container without rebuilding intermediate vectors: ```{r se-input, eval = requireNamespace("SummarizedExperiment", quietly = TRUE)} suppressPackageStartupMessages(library(SummarizedExperiment)) rna_se <- SummarizedExperiment( assays = list(counts = counts), colData = DataFrame(stage = sample_meta$stage, cohort = sample_meta$cohort), rowData = DataFrame(gene_symbol = rownames(counts))) rna_fit_se <- fit_stage_de(rna_se, stage_col = "stage", cohort_col = "cohort") identical(dim(rna_fit_se), dim(rna_fit)) ``` The matrix and SE entry points are functionally equivalent — the S4 method extracts `assay(se, assay_name)` and `colData(se)` and calls the same fitting kernel. The same is true of `fit_stage_de_protein(se, stage_col, cohort_col, assay_name = "intensity")`. # Step 7 — Single-patient trajectory alignment (v0.4.1) `align_patient_profile()` is the single-sample analog of `score_trajectory()`. It takes one patient's per-gene log2(tumor / matched-normal) vector and reports how that profile aligns with each of the atlas's frozen stage-effect axes (`rna_beta_E`, `rna_beta_M`, `rna_beta_L`). **This is an alignment readout, not a stage prediction.** No supervised model has been fit — the atlas is frozen and the function is deterministic given (input, atlas). Misalignment with all three stage axes is itself a finding (it does not imply Normal). ```{r align_synthetic, eval = TRUE} ref <- pdactrace:::.get_reference() high_audit <- ref[!is.na(rna_beta_L) & audit_score >= 0.3][seq_len(300)] # Synthetic patient that resembles a Late-stage profile set.seed(7) patient_rna <- setNames( high_audit$rna_beta_L + rnorm(nrow(high_audit), sd = 0.3), high_audit$gene_symbol) aln <- align_patient_profile(patient_rna, weight_by = "audit_score", min_audit_score = 0.3, top_n_genes = 200) print(aln) aln$rna ``` The output reports two complementary signals: the primary audit-weighted Pearson rho between the patient and each stage axis (with Fisher-z 95% CI), and a secondary audit-weighted vote share showing what fraction of genes are closest to each stage's β. If the two signals disagree, both are surfaced — they are not collapsed into a single decision. # Scope boundary The PDAC reference atlas remains the default reference for `compute_audit_score(evidence = NULL)` and `align_patient_profile()`; the v0.4 release does not provide a `build_user_atlas()` function or a layer plugin architecture. The frozen v0.3.0 scoring rule, the audit labels, and the post-freeze evaluation discipline apply identically whether the input is the bundled atlas or a user-supplied evidence table — the framework is the same; only the data source changes. `align_patient_profile()` deliberately does **not** support supervised stage prediction. A Random-Forest / GLM patient-level classifier trained on TCGA / CPTAC stage labels is conceivable but would break the package's central "no supervised fitting under limited ground truth" discipline; it is therefore out of scope. Limitations: rho is unstable below ~50 overlapping genes; the protein layer reports only categorical concordance because the atlas does not carry per-stage numeric protein effect sizes; the audit-score floor restricts the gene dictionary to anchors with multi-layer evidence — varying it changes which dictionary you align against, not the patient. Disease-agnostic atlas building (Level 2) is planned for v0.5. # Appendix A — Stage harmonization (TNM / AJCC -> Normal / Early / Mid / Late) The tissue layer requires `stage` labels in `c("Normal", "Early", "Mid", "Late")`. Most public PDAC cohorts report TNM (T-stage) or AJCC (I / II / III / IV); the canonical mapping pdactrace uses internally for the bundled atlas is: | Source label | Standardised | Notes | |--------------------------------------|--------------|----------------------------------| | Adjacent normal / matched control | Normal | "non-tumor", "normal pancreas" | | AJCC IA / IB; T1; T1N0M0 | Early | localised, resectable | | AJCC IIA / IIB; T2; T3 N0; T1-2 N1 | Mid | locally advanced, resectable / borderline | | AJCC III; T3 N+; T4 any N | Mid | locally advanced unresectable | | AJCC IV; M1; any T any N M1 | Late | metastatic | If a cohort reports only "stage" without TNM detail, use the AJCC roman-numeral row. If you only have **N vs T** at the tissue layer (no clinical staging), pdactrace's tissue layer does **not** apply -- the 12-template trajectory framework needs at least 2 of the 4 stage levels with samples to discriminate trajectory shape. Three options: 1. Look for clinical staging in a separate metadata file (most TCGA / GEO / CPTAC submissions include it; map via the table above). 2. Use a different tool for tissue analysis and bring the **serum** side through `project_user_serum_cohort()` (Appendix B), which needs only PDAC vs HC. 3. Re-cast all tumor samples to one stage (e.g. all -> "Early"). Technically runs but loses framework value -- the trajectory classification degenerates to "up vs down". Not recommended. # Appendix B -- Serum projection (binary group contrast, no stage) The serum side of the framework is binary: the bundled atlas computes `serum_log2fc_PDAC_vs_HC` and `serum_log2fc_Pan_vs_HC` from group means, not from stage-aware regression. So a user with their own serum DIA / TMT cohort can project it through the framework with only PDAC / HC ± Pancreatitis group labels: ```{r serum_projection_demo, eval = FALSE} # Synthetic example: 50 genes x 24 samples (8 PDAC, 8 HC, 8 Pancreatitis) g <- 50; n <- 24 intensity <- matrix(rnorm(g * n, mean = 5), nrow = g, ncol = n, dimnames = list(paste0("GENE", seq_len(g)), paste0("S", seq_len(n)))) intensity[1:10, 1:8] <- intensity[1:10, 1:8] + 2 # PDAC-up signal coldata <- data.frame( sample = colnames(intensity), group = rep(c("PDAC", "HC", "Pancreatitis"), each = 8)) serum_summary <- project_user_serum_cohort( intensity, coldata, group_col = "group", pdac_label = "PDAC", hc_label = "HC", pan_label = "Pancreatitis", # optional test = "limma_eBayes", # or "wilcox" / "t_test" padj_cutoff = 0.05) head(serum_summary) ``` The output is a `data.table` with `gene_symbol`, `serum_log2fc_PDAC_vs_HC`, `serum_padj_PDAC_vs_HC`, `serum_log2fc_Pan_vs_HC`, `serum_padj_Pan_vs_HC`, `translation_class` (joined with bundled atlas RNA direction: A = same-direction tissue<->serum, B = opposite, C = decoupled), and `serum_detected`. Plug it into the assembler: ```{r serum_assembly_demo, eval = FALSE} ev <- assemble_user_evidence(serum_summary = serum_summary) audit <- compute_audit_score(evidence = ev) ``` If you also have a tissue cohort projected through `project_user_cohort()`, pass both into `assemble_user_evidence()` to get the full multi-layer audit score. ```{r session_info} sessionInfo() ```