> Status: `draft`
>
> Template class: SOURCE-BACKED WORKFLOW

## Purpose

Fit negative-binomial tradeSeq GAMs to raw counts using an explicit
trajectory/pseudotime contract. This notebook consumes pseudotime generated by
a trajectory method; it does not infer a trajectory and does not bin
pseudotime.

## Inputs and parameters

The Seurat object must contain raw counts and metadata columns containing
pseudotime and cell weights.

```{r}
library(Seurat)
library(tradeSeq)
library(BiocParallel)
library(Matrix)
library(qs2)
library(readr)

input_object <- "input/seurat_with_pseudotime.qs2"
output_fit <- "output/tradeseq_fit.qs2"
output_results <- "output/tradeseq_start_end.tsv"

assay_name <- "RNA"
counts_layer <- "counts"
pseudotime_column <- "pseudotime_1"
cell_weight_column <- "weight_1"
min_cells <- 20
feature_set <- NULL
n_knots <- 5
family <- "nb"
n_workers <- 2
test_method <- "start_vs_end"
```

## Analysis

```{r}
object <- qs2::qs_read(input_object)
counts <- GetAssayData(object, assay = assay_name, layer = counts_layer)
metadata <- object[[]]
stopifnot(all(c(pseudotime_column, cell_weight_column) %in% colnames(metadata)))
stopifnot(identical(colnames(counts), rownames(metadata)))

pt <- metadata[[pseudotime_column]]
weights <- metadata[[cell_weight_column]]
stopifnot(all(is.finite(pt)), all(is.finite(weights)), all(weights >= 0))

keep_expressed <- Matrix::rowSums(counts > 0) > min_cells
genes_after_filter <- rownames(counts)[keep_expressed]
genes_use <- if (is.null(feature_set)) {
  genes_after_filter
} else {
  intersect(feature_set, genes_after_filter)
}
stopifnot(length(genes_use) > 0)

pseudotime_matrix <- matrix(
  pt,
  ncol = 1,
  dimnames = list(colnames(object), "lineage_1")
)
weights_matrix <- matrix(
  weights,
  ncol = 1,
  dimnames = list(colnames(object), "lineage_1")
)

BPPARAM <- MulticoreParam(workers = n_workers)
set.seed(1234)
gam_fit <- fitGAM(
  counts = counts[genes_use, , drop = FALSE],
  pseudotime = pseudotime_matrix,
  cellWeights = weights_matrix,
  nknots = n_knots,
  family = family,
  parallel = TRUE,
  BPPARAM = BPPARAM,
  verbose = TRUE
)

stopifnot(test_method %in% c("association", "start_vs_end"))
test_results <- if (test_method == "association") {
  associationTest(gam_fit)
} else {
  startVsEndTest(gam_fit)
}
test_results <- data.frame(
  gene = rownames(test_results),
  test_results,
  row.names = NULL
)
test_results$FDR <- p.adjust(test_results$pvalue, method = "BH")
```

## Diagnostics

```{r}
print(gam_fit)
print(head(test_results[order(test_results$FDR), ], 20))
print(summary(test_results$FDR))
```

## Outputs

```{r}
dir.create(dirname(output_fit), recursive = TRUE, showWarnings = FALSE)
dir.create(dirname(output_results), recursive = TRUE, showWarnings = FALSE)
qs2::qs_save(gam_fit, output_fit)
write_tsv(test_results, output_results)
```

## Method notes

The source notebook uses a one-lineage score as pseudotime, filters genes
expressed in more than 20 cells, uses five knots, and runs start-versus-end
testing. With `feature_set <- NULL`, this template fits all genes passing the
declared expression filter; an explicit feature vector can be supplied when a
smaller hypothesis-driven set is intended. `associationTest()` asks whether a
gene changes along pseudotime, whereas `startVsEndTest()` asks about the
trajectory endpoints. They are not interchangeable. Lineage inference belongs
in a separate trajectory workflow.