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.
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.
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)
#> [1] 50 24
#>
#> CohortA CohortB
#> Early 3 3
#> Late 3 3
#> Mid 3 3
#> Normal 3 3rna_fit <- fit_stage_de(counts,
stage = sample_meta$stage,
cohort = sample_meta$cohort)
#> converting counts to integer mode
#> -- note: fitType='parametric', but the dispersion trend was not well captured by the
#> function: y = a/x + b, and a local regression fit was automatically substituted.
#> specify fitType='local' or 'mean' to avoid this message next time.
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)])
#> gene_symbol rna_pattern rna_pattern_rho lrt_padj
#> <char> <char> <num> <num>
#> 1: TGENE01 Early_Burst_Up 0.9943036 1.409783e-06
#> 2: TGENE02 Early_Burst_Up 0.9864846 2.597462e-09
#> 3: TGENE03 Early_Burst_Up 0.9732333 8.884354e-12
#> 4: TGENE04 Early_Burst_Up 0.9768965 1.368786e-05
#> 5: TGENE05 Early_Burst_Up 0.9561062 5.841372e-05
#> 6: TGENE06 Early_Burst_Up 0.9850767 1.285955e-10prot_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.")
}
#> gene_symbol prot_pattern prot_pattern_rho
#> <char> <char> <num>
#> 1: TPROT01 Mid_Plateau_Up 0.9744093
#> 2: TPROT02 Late_Burst_Up 0.9709155
#> 3: TPROT03 Early_Peak 0.9768480
#> 4: TPROT04 Mid_Peak 0.9771442
#> 5: TPROT05 Early_Burst_Up 0.8545095
#> 6: TPROT07 Late_Loss_Down 0.9275559user_ev <- assemble_user_evidence(
rna_fit = rna_pat,
prot_fit = prot_pat,
signal_peptide = up_genes)
str(user_ev, max.level = 1)
#> Classes 'data.table' and 'data.frame': 100 obs. of 13 variables:
#> $ gene_symbol : chr "TGENE01" "TGENE02" "TGENE03" "TGENE04" ...
#> $ rna_pattern : chr "Early_Burst_Up" "Early_Burst_Up" "Early_Burst_Up" "Early_Burst_Up" ...
#> $ rna_pattern_rho : num 0.994 0.986 0.973 0.977 0.956 ...
#> $ rna_lrt_padj : num 1.41e-06 2.60e-09 8.88e-12 1.37e-05 5.84e-05 ...
#> $ prot_pattern : chr NA NA NA NA ...
#> $ cell_origin_top : chr NA NA NA NA ...
#> $ serum_detected : logi NA NA NA NA NA NA ...
#> $ serum_log2fc_PDAC_vs_HC: num NA NA NA NA NA NA NA NA NA NA ...
#> $ flt_signal_peptide : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
#> $ flt_direction_match : logi NA NA NA NA NA NA ...
#> $ rna_cohort_agreement : num 1 1 1 1 1 1 1 1 1 1 ...
#> $ max_I2_meta : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ max_abs_beta_meta : num NA NA NA NA NA NA NA NA NA NA ...
#> - attr(*, ".internal.selfref")=<pointer: 0x559bebb46bd0>
#> - attr(*, "sorted")= chr "gene_symbol"
#> - attr(*, "source")= chr "user-supplied evidence (assemble_user_evidence)"user_audit <- compute_audit_score(evidence = user_ev)
table(user_audit$audit_class)
#>
#> high_confidence low
#> 11 89
head(user_audit[order(-audit_score)], 10)
#> gene_symbol evidence_strength biological_coherence translational_relevance
#> <char> <num> <num> <num>
#> 1: TGENE01 0.25 0.7471518 0.3
#> 2: TGENE10 0.25 0.7455863 0.3
#> 3: TGENE02 0.25 0.7432423 0.3
#> 4: TGENE06 0.25 0.7425384 0.3
#> 5: TGENE04 0.25 0.7384482 0.3
#> 6: TGENE03 0.25 0.7366166 0.3
#> 7: TGENE08 0.25 0.7332210 0.3
#> 8: TGENE05 0.25 0.7280531 0.3
#> 9: TGENE07 0.25 0.6296328 0.3
#> 10: TGENE09 0.25 0.4416615 0.3
#> leakage_gate heterogeneity_gate positive_score audit_score audit_class
#> <num> <num> <num> <num> <char>
#> 1: 1 1 0.4365031 1.0000000 high_confidence
#> 2: 1 1 0.4359552 0.9987447 high_confidence
#> 3: 1 1 0.4351348 0.9968653 high_confidence
#> 4: 1 1 0.4348884 0.9963008 high_confidence
#> 5: 1 1 0.4334569 0.9930213 high_confidence
#> 6: 1 1 0.4328158 0.9915526 high_confidence
#> 7: 1 1 0.4316273 0.9888299 high_confidence
#> 8: 1 1 0.4298186 0.9846861 high_confidence
#> 9: 1 1 0.3953715 0.9057701 high_confidence
#> 10: 1 1 0.3295815 0.7550496 high_confidenceThe same evaluate_anchor_enrichment() interface accepts
an optional anchors table so a user can evaluate their own
anchor list:
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)
#> tier tier_key top_n n_anchor_declared n_anchor_in_atlas hits
#> <char> <char> <int> <int> <int> <int>
#> 1: SECONDARY (T1+T2) secondary 20 5 5 5
#> expect fold pval
#> <num> <num> <num>
#> 1: 1 5 0.0002059305SummarizedExperiment (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:
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")
#> converting counts to integer mode
#> -- note: fitType='parametric', but the dispersion trend was not well captured by the
#> function: y = a/x + b, and a local regression fit was automatically substituted.
#> specify fitType='local' or 'mean' to avoid this message next time.
identical(dim(rna_fit_se), dim(rna_fit))
#> [1] TRUEThe 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").
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).
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)
#> Best-match by rho: Late (rho=0.98, p=7e-207, 200 genes); voting agrees: Late (vote_share=0.46).
#>
#> stage cor_to_stage_axis cor_pval cor_lo95 cor_hi95 vote_share
#> <char> <num> <num> <num> <num> <num>
#> Normal NA NA NA NA 0.1101468
#> Early 0.9557575 8.976376e-156 0.9419171 0.9663571 0.2036359
#> Mid 0.9598262 7.110716e-164 0.9472241 0.9694662 0.2259604
#> Late 0.9751032 6.960526e-207 0.9672132 0.9811127 0.4602569
#> weighted_dist n_genes_used
#> <num> <int>
#> 1.0627099 200
#> 0.3116161 200
#> 0.2939051 200
#> 0.2362621 200
#>
#> [atlas v0.4.0; 200/300 input genes used]
aln$rna
#> stage cor_to_stage_axis cor_pval cor_lo95 cor_hi95 vote_share
#> <char> <num> <num> <num> <num> <num>
#> 1: Normal NA NA NA NA 0.1101468
#> 2: Early 0.9557575 8.976376e-156 0.9419171 0.9663571 0.2036359
#> 3: Mid 0.9598262 7.110716e-164 0.9472241 0.9694662 0.2259604
#> 4: Late 0.9751032 6.960526e-207 0.9672132 0.9811127 0.4602569
#> weighted_dist n_genes_used
#> <num> <int>
#> 1: 1.0627099 200
#> 2: 0.3116161 200
#> 3: 0.2939051 200
#> 4: 0.2362621 200The 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.
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.
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:
project_user_serum_cohort() (Appendix B), which needs only
PDAC vs HC.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:
# 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:
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.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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.32.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] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [3] GenomicRanges_1.65.4 Seqinfo_1.3.2
#> [5] IRanges_2.47.5 S4Vectors_0.51.10
#> [7] BiocGenerics_0.59.12 generics_0.1.4
#> [9] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [11] data.table_1.18.6.1 pdactrace_0.99.25
#> [13] rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 ellipse_0.5.0 xfun_0.61
#> [4] bslib_0.12.0 ggplot2_4.0.3 ggrepel_0.9.8
#> [7] lattice_0.23-1 vctrs_0.7.3 tools_4.6.1
#> [10] parallel_4.6.1 tibble_3.3.1 rARPACK_0.11-0
#> [13] pkgconfig_2.0.3 Matrix_1.7-6 RColorBrewer_1.1-3
#> [16] S7_0.2.2 mixOmics_6.37.0 lifecycle_1.0.5
#> [19] stringr_1.6.0 compiler_4.6.1 farver_2.1.2
#> [22] statmod_1.5.2 DESeq2_1.53.4 codetools_0.2-20
#> [25] htmltools_0.5.9 sys_3.4.3 buildtools_1.0.0
#> [28] sass_0.4.10 yaml_2.3.12 tidyr_1.3.2
#> [31] pillar_1.11.1 jquerylib_0.1.4 MASS_7.3-66
#> [34] BiocParallel_1.47.0 limma_3.99.0 DelayedArray_0.39.7
#> [37] cachem_1.1.0 abind_1.4-8 RSpectra_0.16-2
#> [40] locfit_1.5-9.12 tidyselect_1.2.1 digest_0.6.39
#> [43] stringi_1.8.9 purrr_1.2.2 reshape2_1.4.5
#> [46] dplyr_1.2.1 maketools_1.3.2 labeling_0.4.3
#> [49] fastmap_1.2.0 grid_4.6.1 cli_3.6.6
#> [52] SparseArray_1.13.3 magrittr_2.0.5 S4Arrays_1.13.1
#> [55] corpcor_1.6.10 withr_3.0.3 scales_1.4.0
#> [58] XVector_0.53.0 igraph_2.3.3 otel_0.2.0
#> [61] gridExtra_2.3.1 evaluate_1.0.5 knitr_1.52
#> [64] rlang_1.3.0 Rcpp_1.1.2 glue_1.8.1
#> [67] jsonlite_2.0.0 plyr_1.8.9 R6_2.6.1SummarizedExperiment
(v0.99.2)