---
title: "5. Reproducibility from raw data"
author: "Jibeom Ko"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{5. Reproducibility from raw data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# Goal

`pdactrace` ships a frozen 10,113-gene atlas plus the audit, query,
report, and projection API that consumes it. This vignette
documents **four reproducibility layers**, ordered from "no
external data needed" to "the entire FASTQ-to-atlas pipeline":

1. **Layer 1 — offline, bundled.** Use the prebuilt `data/*.rda`
   directly. No network, no rebuild.
2. **Layer 2 — user cohort.** Project a user-supplied count or
   intensity matrix through the same audit framework via
   `project_user_cohort()`.
3. **Layer 3 — bundled phase chain.** Re-derive the bundled atlas
   from a git checkout of the source repository, which carries
   every processed input the build reads.
4. **Layer 4 — full pipeline (FASTQ → counts → fits).**
   Out of scope for the package; this layer is the
   manuscript-monorepo's job. Pointer at the bottom.

# Layer 1 — Offline, bundled

```{r layer1, eval = FALSE}
library(pdactrace)

# All bundled atlas objects load lazily from `data/`:
data(pdactrace_reference)
data(pdactrace_protein_betas)
data(default_templates)
data(pdactrace_external_anchors)
data(meta_analysis)
data(atlas_metadata)
data(pdactrace_data_sources)

query_gene("LTBP1")             # uses pdactrace_reference internally
explain_score("LTBP1")          # same
list_data_sources(layer = "RNA")
atlas_provenance()
```

This is the canonical use case — most users never need anything
beyond Layer 1.

# Layer 2 — User cohort

If you have your own staged RNA-seq counts (and optionally a
matched tissue-protein intensity matrix), feed it into the same
framework:

```{r layer2, eval = FALSE}
data(toy_counts); data(toy_coldata); data(toy_protein)

res <- project_user_cohort(
  rna        = toy_counts,
  coldata    = toy_coldata,
  stage_col  = "stage",
  cohort_col = "cohort",
  protein    = toy_protein)

res$summary
head(res$audit[order(-audit_score)])
```

`project_user_cohort()` chains four steps internally:
`fit_stage_de()` → `classify_trajectory()` → `assemble_user_evidence()`
→ `compute_audit_score()`. Each can be called individually for
custom workflows; see `vignette("user_cohort_extension")`.

`SummarizedExperiment` input works identically — pass an `SE` whose
`assay()` is the count matrix and whose `colData()` carries the
stage / cohort columns (`fit_stage_de(se, stage_col, cohort_col)`).

# Layer 3 — Re-derive the bundled atlas

`data-raw/build_reference.R` rebuilds the atlas from processed,
gene-level inputs. Seven of them are bundled in the package as
`inst/extdata/*.csv.xz` (`multi_cohort_consistency` plus the six
downstream tables `phase{2c,29,42,60,77,80}`). The two larger upstream
fits — `phase33` (RNA) and `phase34` (protein), ~3.5 MB together — are
kept out of the package tarball and live in `data-raw/` of the source
repository. `data-raw/` is not installed with the package, so the
rebuild runs from a git checkout:

```{r layer3, eval = FALSE}
# From a shell:
#   git clone https://github.com/jibeomko/pdactrace
#   cd pdactrace
#   Rscript data-raw/build_reference.R
```

Each input is looked up in `inst/extdata/` and then in `data-raw/`, so
no environment variable or second repository is needed. Run this way,
`build_reference.R` reproduces the 64 base columns of the bundled
`pdactrace_reference` exactly, for the same 10,113 genes in the same
order. The derived audit, Pareto and TRACE-D columns are added by the
`data-raw/attach_*.R` scripts.

If you only want the two upstream fits, without cloning,
`download_phase_csvs()` fetches them into your `BiocFileCache`:

```{r layer3-download, eval = FALSE}
paths <- download_phase_csvs("both")
p33 <- data.table::fread(paths[["phase33"]])
```

# Layer 4 — From FASTQ

This layer — raw FASTQ → STAR/salmon counts → DESeq2 LRT → phase33
CSV — is **out of scope for `pdactrace`**. The cohort-specific
pipelines, ~110 numbered `phase*` scripts, FragPipe configs, scVI /
liana notebooks and ~100 GB of intermediate artefacts live in a
separate manuscript working repository that is not public. Every input
dataset is public, and `list_data_sources()` lists the accessions so
the raw data can be retrieved independently.

# Provenance summary

```{r prov, eval = FALSE}
ap <- atlas_provenance()
ap$package_doi
ap$manuscript_doi
ap$data_layers
list_data_sources()
```

Every gene's `provenance` field is renderable as a plain-English
evidence summary via [`format_provenance()`].

# Session information

```{r session-info}
sessionInfo()
```