08 可复现分析报告与结果交付

分析做完不是终点，能被别人（包括 6 个月后的自己）重新跑一遍才是。这一章讨论 bulk RNA-seq 项目的收尾环节：怎么组织代码和输出、怎么写分析报告、怎么交付给合作者。没有什么"最佳实践"一定对，但有一些经验能让你后续工作顺得多。

可复现的三个层次

从最基本到最严格：

层次	含义	实现方式
相同数据 + 相同代码 → 相同结果	基础要求	代码版本控制 + 随机种子
相同数据 + 相同环境 → 相同结果	跨机器复现	conda / renv / docker 固定依赖
相同实验 + 相同分析流程 → 相似结论	最高要求	清晰的样本 metadata + 分析日志

前两条是技术问题；第三条是组织问题。

目录结构回顾

01 章里给过一份推荐的目录结构。在收尾阶段加一件事：把最终的可交付产物单独放进 deliverables/。

rnaseq_project/
├── README.md
├── metadata/samples.tsv
├── data/
│   ├── raw/               
│   ├── qc/                
│   └── counts/            
├── analysis/              
│   ├── 01-qc.R
│   ├── 02-de.R
│   ├── 03-enrichment.R
│   └── 04-figures.R
├── results/               # 中间结果
│   ├── tables/
│   └── figures/
├── deliverables/          # 交付给合作者的最终产物
│   ├── report.html
│   ├── report.pdf
│   ├── figures_final/     # 发表级图
│   └── tables_final/      # 带 README 的表格
└── envs/
    └── env.yaml / renv.lock

results/ 是 "所有代码跑出来的东西"；deliverables/ 是 "我要给别人看的那一小部分"。两者要分开，因为 results/ 里通常有很多中间文件、测试用的图、探索性结果，不适合直接给合作者。

R Markdown / Quarto 做分析报告

bulk RNA-seq 最常见的交付形式是一份 HTML 报告，包含样本表 + QC 图 + 差异基因结果 + 富集 + 讨论。R Markdown（传统）或 Quarto（新）是目前最成熟的方案。

一个最小可用的 R Markdown 模板：

---
title: "RNA-seq analysis of dex-treated airway epithelium"
author: "Your Name"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: show
    theme: flatly
---

# Background

(简要描述实验目的、样本来源、分析目标)

# Sample overview

```{r samples, message=FALSE}
library(DT)
samples <- read.delim("metadata/samples.tsv")
datatable(samples, rownames = FALSE)

QC

knitr::include_graphics("results/figures/01-library-size.png")

... PCA、距离热图等

Differential expression

library(DESeq2); library(dplyr); library(DT)
# ...

Enrichment

...

Session info

sessionInfo()

几个实用技巧：

- `code_folding: show`：读者能默认看到代码，但可以随时折叠
- `toc: true` + `toc_float: true`：左侧浮动目录，长报告必备
- `DT::datatable()` 嵌入交互式表格：合作者可以自己搜索和排序
- 最后一定要有 `sessionInfo()`：记录所有包版本

**Quarto** 是 R Markdown 的继任者，支持多语言（R、Python、Julia），渲染更快。生态还在迁移期，但新项目推荐直接用 Quarto（`.qmd` 文件），语法几乎完全兼容 `.Rmd`。

## 包版本管理

"我半年前跑过的代码，今天重跑结果变了"是一个很常见的困扰。R 这边最主流的两个方案：

### renv

`renv` 会在项目里记录一份 lockfile，锁定所有包的版本：

```r
install.packages("renv")
renv::init()        # 初始化，生成 renv.lock
renv::snapshot()    # 任何包版本变了，重新 snapshot
renv::restore()     # 别人拿到项目，装回原来的所有包版本

对于发表前要归档的项目，强烈建议用 renv。

conda env.yaml

如果项目里混了 bash 工具（Salmon、STAR、FastQC）和 R 包，用 conda 统一管理更方便：

# envs/env.yaml
name: bulkrnaseq
channels:
  - conda-forge
  - bioconda
dependencies:
  - r-base=4.3.2
  - bioconductor-deseq2=1.42.0
  - bioconductor-edger=4.0.16
  - bioconductor-limma=3.58.1
  - bioconductor-clusterprofiler=4.10.0
  - salmon=1.10.2
  - fastqc=0.12.1
  - multiqc=1.21
  - r-tidyverse

用 conda env create -f envs/env.yaml 或者更快的 mamba env create -f envs/env.yaml，就能在任何机器上装回一样的环境。

交付给合作者的结果表

给合作者的 CSV / Excel 表不能只是 write.csv(res, ...)。至少要做三件事：

1. 加一份 README

每张表都配一个 _README.txt 或 _README.md，说明列的含义：

de_genes_airway.tsv
====================

8 个 airway 样本，trt (dex-treated) vs untrt (control)。
用 DESeq2 v1.42.0 做差异分析，design = ~ cell + dex。

列：
  gene_id          — Ensembl gene ID
  gene_symbol      — HGNC symbol（可能为空）
  baseMean         — DESeq2 归一化后的平均表达量
  log2FoldChange   — apeglm shrunken log2FC，正值 = 在 trt 里更高
  lfcSE            — LFC 的标准误
  pvalue           — raw p-value from Wald test
  padj             — Benjamini-Hochberg 校正后的 q-value

分析者：xxx
分析日期：2024-xx-xx

2. 加入 gene_symbol

矩阵里基因 ID 一般是 Ensembl，但合作者关心的是 SYMBOL。一定要把 SYMBOL 合并进去再交付：

library(org.Hs.eg.db); library(AnnotationDbi)
id_map <- AnnotationDbi::select(
  org.Hs.eg.db,
  keys = rownames(res), keytype = "ENSEMBL",
  columns = c("SYMBOL")
)
res_out <- merge(as.data.frame(res), id_map, by.x = 0, by.y = "ENSEMBL")

3. 分层交付

别把所有 60000 行都给合作者。分成三张：

de_all_genes.tsv — 全部基因，用于后续分析
de_significant.tsv — padj < 0.05 的，通常几百到几千行
de_top_genes.xlsx — top 50 padj + 显著通路里的关键基因，带格式化的 Excel

Methods 段怎么写

投稿时 Methods 段的 RNA-seq 部分通常要覆盖：

样本来源：几个对照 / 处理、生物学重复数、测序平台
比对 / 定量工具和版本：STAR v2.7.11、Salmon v1.10.2 等
参考基因组版本：GRCh38, GENCODE v44
定量方式：featureCounts 或 tximport from Salmon
DE 工具 + 设计公式 + 阈值：DESeq2 v1.42.0, design = ~ batch + condition, padj < 0.05 & |LFC| > 1
富集工具 + 数据库版本：clusterProfiler v4.10.0, GO BP, MSigDB Hallmarks v2023.2

两三段话就够。不要模糊处理（别写"我们做了标准的差异分析"），每个工具都要带版本号。审稿人常问的是"你用的 DESeq2 是哪个版本，padj 阈值是多少"。

常见坑

最后几条经验：

随机种子：任何用到随机算法的地方（UMAP、聚类、抽样）都 set.seed(42)，不然每次跑结果都不一样
绝对路径：别在代码里写 /Users/你的名字/xxx，用 here::here() 或 file.path() 相对路径
数据不进 git：raw FASTQ / BAM / 大 counts 文件别 commit，用 .gitignore 屏蔽。如果要共享，用 Zenodo / figshare / GSA / SRA
日志：长分析用 sink() 把 sessionInfo() 和关键输出存一份
测试：bulk07 那个"三工具对比"本身就是一种测试 —— 新分析完成后用不同工具验证一遍，结果差太多时回头查

工具小清单

每个环节推荐的工具：

任务	推荐
报告	`R Markdown` → 过渡到 `Quarto`
表格交互	`DT::datatable`
环境	`renv`（纯 R）或 `conda`（混合工具）
项目模板	`usethis::create_project` + `here`
代码风格	`styler` + `lintr`
git 钩子	`precommit` 自动格式化

下载资源

本章不配脚本。给合作者的交付物相关模板（README / Methods 段 / R Markdown 骨架）可以从 BioF3 GitHub 仓库找到（规划中）。

下一步

bulk RNA-seq 专栏 8 个模块到这里就完成了。后续可能扩展：

转录本级别分析（DTU、DTE）
拼接与可变剪接
小 RNA / miRNA / lncRNA 专题
单细胞 pseudobulk 与 bulk 的互补

也可以回到单细胞实践教程或其他组学专栏。

参考资源

静态文件

离线资料下载

手册 HTML / PDF 已在后台预生成，点击后直接下载网站静态资源。

本篇 HTML 本篇 PDF说明

可复现的三个层次​

目录结构回顾​

R Markdown / Quarto 做分析报告​

QC

Differential expression

Enrichment

Session info

conda env.yaml​

交付给合作者的结果表​

1. 加一份 README​

2. 加入 gene_symbol​

3. 分层交付​

Methods 段怎么写​

常见坑​

工具小清单​

下载资源​

下一步​

参考资源​