生信方法入门：把 RNA-seq 学习项目整理成可以复查的目录

生信分析中的困扰不总是“软件不会用”。几周后重新打开项目，可能已经记不清哪张表是输入、哪次运行修改了参数、某个样本为什么被排除。目录结构不能替代方法学，但能为查清这些问题提供基础。

这篇只搭建一个教学项目骨架，用六个虚构样本和三行人工计数练习元数据检查与日志记录。不下载真实 FASTQ，不调用比对或差异分析软件，不生成任何可用于解释真实生物学问题的结论。

1. 将来源、决定和结果分开

建议把原始输入、样本信息、配置、脚本、结果与日志分别保存。目录名字本身没有统一标准，关键是每一层职责清楚，其他人可以看懂。

1
2
3
4
5
6
7
8
9
project/
  inputs/raw/        真实输入的存放位置；本练习留空
  references/        参考文件与版本说明；本练习留空
  metadata/          样本表
  config/            参数与处理约定
  workflow/          正式脚本及执行入口
  results/           可重新生成的产物
  logs/              运行与检查记录
  env/               运行环境信息

不要把未处理文件和经过过滤的文件都叫作 final.csv。也不要一边改输入，一边重跑同一个结果目录，却没有记录输入何时发生变化。

元数据是输入的一部分，不是可有可无的备注。至少需要稳定的样本 ID、实验条件、个体或生物学重复关系；视实验设计还需要批次、组织、采样时间、建库协议及其他协变量。真实患者信息不应直接写进公开仓库，匿名化标识也需要按适用的数据管理要求维护映射。

2. 完整示例：创建全新的教学项目

代码只使用 Python 标准库。保存为 create_learning_project.py，执行 python create_learning_project.py。程序使用 tempfile 创建新目录，所有写入都位于该目录；结束后保留文件供检查，不自动删除或访问现有项目。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
import csv
import hashlib
import json
import platform
import sys
import tempfile
from datetime import datetime, timezone
from pathlib import Path

root = Path(tempfile.mkdtemp(prefix="perry-rnaseq-"))
for relative in ["inputs/raw", "references", "metadata", "config", "workflow",
                 "results", "logs", "env"]:
    (root / relative).mkdir(parents=True, exist_ok=True)

samples = [
    {"sample_id": "C1", "subject_id": "demo_01", "condition": "control", "batch": "B1"},
    {"sample_id": "C2", "subject_id": "demo_02", "condition": "control", "batch": "B1"},
    {"sample_id": "C3", "subject_id": "demo_03", "condition": "control", "batch": "B1"},
    {"sample_id": "T1", "subject_id": "demo_04", "condition": "treated", "batch": "B1"},
    {"sample_id": "T2", "subject_id": "demo_05", "condition": "treated", "batch": "B1"},
    {"sample_id": "T3", "subject_id": "demo_06", "condition": "treated", "batch": "B1"},
]
sample_ids = [record["sample_id"] for record in samples]
metadata_path = root / "metadata/samples.tsv"
with metadata_path.open("w", encoding="utf-8", newline="") as handle:
    writer = csv.DictWriter(handle, fieldnames=list(samples[0]), delimiter="\t")
    writer.writeheader()
    writer.writerows(samples)

counts_path = root / "inputs/toy_counts.tsv"
rows = [
    ["teaching_gene_A", 10, 11, 12, 20, 21, 22],
    ["teaching_gene_B", 5, 6, 4, 7, 8, 6],
    ["teaching_gene_C", 1, 2, 3, 1, 2, 3],
]
with counts_path.open("w", encoding="utf-8", newline="") as handle:
    writer = csv.writer(handle, delimiter="\t")
    writer.writerow(["gene_id", *sample_ids])
    writer.writerows(rows)

config = {
    "schema_version": 1,
    "data_kind": "artificial_teaching_counts",
    "metadata": "metadata/samples.tsv",
    "counts": "inputs/toy_counts.tsv",
    "reference_release": None,
    "strandedness": "not_applicable_to_this_toy_table",
    "analysis": "column_totals_only; no differential testing",
}
(root / "config/project.json").write_text(
    json.dumps(config, indent=2, ensure_ascii=False) + "\n", encoding="utf-8"
)

# Read the files back: validation must check written input, not only Python variables.
with metadata_path.open(encoding="utf-8", newline="") as handle:
    records = list(csv.DictReader(handle, delimiter="\t"))
with counts_path.open(encoding="utf-8", newline="") as handle:
    table = list(csv.reader(handle, delimiter="\t"))

metadata_ids = [record["sample_id"] for record in records]
count_ids = table[0][1:]
assert len(set(metadata_ids)) == len(metadata_ids), "Duplicate sample ID"
assert metadata_ids == count_ids, "Metadata and count columns differ in order"
assert all(record["condition"] in {"control", "treated"} for record in records)
gene_ids = [row[0] for row in table[1:]]
assert len(set(gene_ids)) == len(gene_ids), "Duplicate gene ID"
totals = [0] * len(count_ids)
for row in table[1:]:
    assert len(row) == len(count_ids) + 1, "Wrong number of columns"
    assert all(value.isdecimal() for value in row[1:]), "Toy counts must be non-negative integers"
    for index, value in enumerate(row[1:]):
        totals[index] += int(value)

with (root / "results/toy_column_totals.tsv").open("w", encoding="utf-8", newline="") as handle:
    writer = csv.writer(handle, delimiter="\t")
    writer.writerow(["sample_id", "toy_column_total"])
    writer.writerows(zip(count_ids, totals))

manifest = {}
for relative in ["metadata/samples.tsv", "inputs/toy_counts.tsv", "config/project.json"]:
    manifest[relative] = hashlib.sha256((root / relative).read_bytes()).hexdigest()
(root / "config/input_sha256.json").write_text(
    json.dumps(manifest, indent=2) + "\n", encoding="utf-8"
)
environment = {"python": sys.version, "platform": platform.platform()}
(root / "env/python_environment.json").write_text(
    json.dumps(environment, indent=2) + "\n", encoding="utf-8"
)
log = {
    "recorded_at_utc": datetime.now(timezone.utc).isoformat(),
    "status": "PASS",
    "samples": len(count_ids),
    "genes": len(gene_ids),
    "analysis": "teaching input validation and column totals only",
}
(root / "logs/check.json").write_text(json.dumps(log, indent=2) + "\n", encoding="utf-8")
print("Project:", root)
print("PASS: samples=6 genes=3")
print("Toy column totals:", dict(zip(count_ids, totals)))

assert 在这里用于小练习的内部检查。运行时不要使用 Python 的 -O 优化选项，否则这些检查会被移除；正式流水线可以改成显式条件与异常，让校验成为不可跳过的执行步骤。

3. 核对几个真正有意义的结果

程序应给出六个样本、三个教学基因，列和依次为 C1=16、C2=19、C3=19、T1=28、T2=31、T3=31。它们只是人工整数的算术结果，不是测序深度质量判断，也不是差异表达证据。

打开 metadata/samples.tsv 与 inputs/toy_counts.tsv，确认样本顺序一致。真实工作中，一张表按照字母排序，另一张保持原始顺序，如果直接按位置拼接，可能悄悄把标签分配给错误样本。应按样本 ID 对齐并检查，不要只比较列数。

还要打开 logs/check.json，检查状态是否为 PASS；打开 config/project.json，确认它明确说明“教学计数、只求列和”。这能避免后来的人将演示产物误认为完成了真实分析。

SHA-256 用于检查记录的文件内容是否发生变化。它不证明数据来源可信，也不证明分析正确；如果文件变化，就需要说明为什么变化，并更新对应运行记录，而不是悄悄重写旧清单。算法接口见 Python 的 hashlib 官方文档。

4. 目录准备好了，还缺哪些真实分析信息

正式 RNA-seq 流程通常还需要数据来源与授权、读长与单端/双端信息、链特异性、质控规则、参考基因组及注释的精确版本、软件版本、线程和参数，以及计数或定量的具体方法。把这些记录在配置和环境文件中，比把软件名字写进 README 更具体。

原始整数计数与 TPM 不是可以任意互换的输入。本篇没有进行差异检验；如继续学习 DESeq2，需要按它要求的计数输入和实验设计构造对象，而不是把 TPM、VST 或 rlog 值当作原始计数交进去。基于转录本定量结果的导入也有专门流程，不能靠四舍五入 TPM 代替。详见 DESeq2 官方教程。

本例的六个 subject_id 都是不同的虚构个体。如果真实数据来自配对设计、重复时间点或复杂批次，样本表与统计设计都必须体现这些关系。复制一个模板不会自动得到正确的设计矩阵。

5. 环境文件和日志应该怎样补全

示例只记录了 Python 与操作系统，不包含真实生信软件清单。以后引入依赖时，可以在独立环境中安装并记录实际版本；Python 的环境隔离可参考 venv 官方文档。虚拟环境目录本身通常不是需要直接复制给其他人的交付物，版本、安装说明和复建方法更重要。

如果使用 R，应保留实际会话的 sessionInfo()；如果使用命令行工具，要记录完整命令、工作目录、起止时间、退出码，以及标准输出和错误日志。容器或锁定依赖有助于减少环境差异，但它们不能代替参考数据版本和分析参数。

将本次运行的脚本副本保存进 workflow/，并在自己的说明中记录执行入口。上面的校验清单只包含三个教学输入文件，没有自动跟踪后来添加的脚本；添加脚本或参考文件后，应扩展清单，而不是误以为它已经覆盖整个项目。

6. 常见误区与下一步练习

常见误区包括：结果表只有截图没有原始表；同一个样本在不同文件里采用不同名字；把技术重复当成独立生物学重复；只记“用了最新参考基因组”；把实际失败的运行与成功结果混进同一个日志目录。

可复查还需要解释性记录：排除某个样本的原因是什么？为什么改参数？新的运行相较旧运行改变了哪些输入？应写出判断依据，不只写“重新运行后正常”。

先做一个低风险练习：只在新教学目录里，把计数表的 T1 与 T2 两列交换，再执行下面的只读校验，确认样本顺序检查能够发现变化。再把一个样本 ID 改成重复值，观察唯一性检查是否触发。不要在真实输入上做这样的破坏性测试。

将检查保存为 check_saved_project.py，以创建脚本打印的完整目录为唯一参数运行。不要为了检查改过的文件而重跑创建脚本：它会生成另一个新目录，而不是读取你编辑的那一个。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import csv
import sys
from pathlib import Path

if len(sys.argv) != 2:
    raise SystemExit("Usage: python check_saved_project.py PROJECT_DIRECTORY")
root = Path(sys.argv[1]).resolve(strict=True)
with (root / "metadata/samples.tsv").open(encoding="utf-8", newline="") as handle:
    samples = list(csv.DictReader(handle, delimiter="\t"))
with (root / "inputs/toy_counts.tsv").open(encoding="utf-8", newline="") as handle:
    counts = list(csv.reader(handle, delimiter="\t"))
ids = [sample["sample_id"] for sample in samples]
assert len(ids) == len(set(ids)), "Duplicate sample ID"
assert ids == counts[0][1:], "Sample columns do not match metadata order"
print("PASS: sample IDs are unique and ordered correctly")

最后给项目写一份很短的 README：这是教学数据还是正式数据、如何重新生成结果、成功时应检查哪些文件、哪些步骤尚未运行。这四句话能够阻止很多误解，也为之后接入真实 RNA-seq 工作流保留清晰边界。

官方参考

继续阅读与公开资源

NCBI FASTQ 格式说明

官方格式参考

打开公开资源 ↗

FastQC 项目文档

官方质量控制资源

打开公开资源 ↗