DNA 序列数据分析与处理：逐步操作手册

你要完成的工作，是把一个序列文件变成来源清楚、坐标正确、标签可信、能交给模型的数据集。阅读顺序按数据流组织。前九章用 CPU 即可完成。

本项目有两条实践线。大肠杆菌基因组用于学习“序列如何与注释、方向、逐碱基标签对应”；人类启动子数据用于学习“短序列分类数据如何检查、去重、编码和评估”。两套任务的标签分别构建。

0. 运行位置与环境

打开 PowerShell，复制以下命令。

Set-Location-LiteralPath'G:\000_school\001_科研\1_全基因组选择测序\0_资料\20261013TrinityDNA\TrinityDNA_learning'$py='.\.venv\Scripts\python.exe'&$py-m dna_learning--help

应看到examples、inspect、download、prepare、extract-region、preview、report、baseline、check-env、train等命令。

$py是本教程使用的 Python。直接指定它可避免调用到 ArcGIS 的 Python 2.7 或其他环境。每次新开 PowerShell，重新运行前两行即可。

项目内的.venv已在当前电脑建立。若将项目拷贝到另一台电脑，应重新创建环境，不直接复制虚拟环境。

&'G:\Python\python.exe'-m venv.venv &$py-m pip install-r requirements-data.txt &$py-m dna_learning check-env

GPU 模型章另装 PyTorch；数据处理阶段无需等待模型依赖。

本次环境torch固定为Windows/Python3.12的CUDA12.6构建，迁移到其他系统时应按PyTorch官方兼容表选择对应构建，数据处理依赖可独立安装。

1. 先学会辨认文件

1.1 生成可反复练习的小文件

&$py-m dna_learning examplesGet-Content-LiteralPath'data\examples\toy.fasta'Get-Content-LiteralPath'data\examples\toy.fastq'

FASTA 示例为：

>toy_chr 人工教学序列 acgtACGTNNRYACGTACGT

>后第一项是序列 ID，其余内容是描述。换行只用于排版，不能把每行当作独立序列。一个文件可包含多条以>开头的记录。

FASTQ 示例为：

@toy_read_1 artificial_Phred33 ACGTNACG + IIII!5?I

FASTQ格式是一种基于ASCII编码的文本格式，主要用于存储生物序列及其测序质量评价数据，为高通量测序领域的标准格式，并被NCBI Short Read Archive等数据库采用。
该格式包含四行文本：首行以@开头记录序列标识符及描述信息，第二行为碱基序列，第三行以+起始可附带注释，第四行通过特定ASCII字符表示每个碱基的测序质量评分。 [1]FASTQ的质量评分采用Phred算法，将错误概率转换为ASCII字符（如Q30对应字符’?'）。Illumina测序仪生成的FASTQ文件包含flowcell、lane、tile等设备位置参数。
该格式可通过Biopython、EMBOSS等工具与FASTA等格式互相转换，适用于Hiseq系列测序平台的数据存储与分析。

这里采用 Phred+33。I对应 Q40，!对应 Q0，5对应 Q20，?对应 Q30。质量描述的是测序读段中每个碱基的可信度；序列与质量字符数必须相同。一般 FASTQ 可有换行，正式读取使用解析器，不能一律按四行手工分组。

在二代/高通量测序（FASTQ 格式）中，Q40、Q30、Q20 这里的 Q 指的是 Phred 质量得分（Phred Quality Score），它是用来量化测序仪对某个碱基识别的“准确率”或“测错概率”的指标。

&$py-m dna_learning inspectdata/examples/toy.fastq

预期format为fastq，record_count为 1，长度为 8，平均 Phred 为 31.25。这是格式教学样例，尚未执行测序读段去接头或比对。

1.2 同时认识其他文件的职责

&$py-m dna_learning inspectdata/examples/toy.gff3 &$py-m dna_learning inspectdata/examples/toy.bed &$py-m dna_learning inspectdata/examples/toy.vcf

这些小文件互相独立，用于辨认格式。压缩后缀.gz说明文件被 gzip 压缩，不能据此判断内部是序列还是注释。

2. 下载真实数据并建立来源清单

&$py-m dna_learning download

该命令下载大肠杆菌 FASTA、GenBank、GFF3，酵母 FASTA，以及启动子数据的三个 CSV。成功时逐项打印“已校验”。再次运行时，已通过校验的数据直接复用。

Get-Content-LiteralPath'data\raw\manifest.json'

重点看assembly、url、bytes和校验值。大肠杆菌版本为GCF_000005845.2，酵母为GCF_000146045.2。

3. 检查序列的基本质量

&$py-m dna_learning inspectdata/raw/ecoli/GCF_000005845.2_ASM584v2_genomic.fna.gz--out reports/ecoli_inspect.json

预期包含 1 条记录，总长度4,641,652 bp(1 bp 代表 “一对配对的碱基”,1 个字母 = 1 bp)，序列 ID 为NC_000913.3。这些数值是本次真实文件的检查结果。

逐项理解以下指标。

核心代码在dna_learning/io.py的normalize与base_stats。示例：

&$py-c"from dna_learning.io import normalize,base_stats; s,changes=normalize('aCGtNRY'); print(s); print(changes); print(base_stats(s))"

预期规范化序列为ACGTNNN，长度仍为 7，R/Y 的位置为[5, 6]。原始文件不改动。若直接把 N 删除，后续所有注释坐标都会发生偏移。

本项目把合法模糊碱基映射为 N；-、?等非法字符会报错，提示检查数据来源或格式。GC 异常用于提示进一步调查，单凭这一项不能确定污染。

4. 理解坐标：最容易造成静默错误的一步

本项目内部统一采用 Python 风格的0 基、左闭右开[start0,end0)。GFF3 的原始起止坐标采用 1 基闭区间；转换为start0=start1-1、end0=end1。

序列 A C G T A A C T G G C A 1基位置 1 2 3 4 5 6 7 8 9 10 11 12 0基索引 0 1 2 3 4 5 6 7 8 9 10 11

原始注释3..9对应 Python[2:9]，提取GTAACTG，长度为 7。

&$py-c"s='ACGTAACTGGCA'; print(s[2:9]); print(len(s[2:9]))"

Biopython 读取 GenBank 后已经转换了位置对象，不能再对feature.location.start减一次 1。features.tsv同时提供location_biopython0、0 基分段坐标和换算后的 1 基坐标，原始注释文件另行保留。

&$py-m dna_learning extract-region--fastadata/raw/ecoli/GCF_000005845.2_ASM584v2_genomic.fna.gz--seqid NC_000913.3--start0 0--end0 300--strand+--out reports/first_300_plus.fasta &$py-m dna_learning extract-region--fastadata/raw/ecoli/GCF_000005845.2_ASM584v2_genomic.fna.gz--seqid NC_000913.3--start0 0--end0 300--strand---out reports/first_300_minus.fasta

这两组命令的核心作用是：从大肠杆菌参考基因组中分别提取同一区间的「正链」和「负链」序列，并用代码验证它们之间严格的“反向互补”关系。

两条结果长度都为 300；负链结果是同一参考区间的反向互补。DNA 互补关系为 A↔T、C↔G，N 保持。

&$py-c"from dna_learning.io import read_fasta,reverse_complement; p=next(iter(read_fasta('reports/first_300_plus.fasta').values())); m=next(iter(read_fasta('reports/first_300_minus.fasta').values())); print(m==reverse_complement(p)); print(reverse_complement(m)==p)"

两个结果都应为True。

5. 用注释提取 CDS，并反查生物学含义

CDS 是基因组中真正能够被细胞“翻译”成蛋白质氨基酸的那一段连续的核苷酸序列。

&$py-m dna_learning prepare--only genomes

产物位于data/processed/ecoli：

本次提取4,318 条 CDS。4,297 条完整且无相关例外标记的 CDS 翻译匹配，21 条带例外标记而跳过直接核对。跳过原因可从表中追溯。

Get-Content-LiteralPath'data\processed\ecoli\features.tsv'-TotalCount 4Get-Content-LiteralPath'data\processed\ecoli\translation_checks.tsv'-TotalCount 8

提取核心使用feature.extract(record.seq)，让 Biopython 处理负链和复合位置。翻译根据注释的transl_table与codon_start；完整 CDS 使用cds=True，从而正确处理细菌的替代起始密码子。

这里区分 gene、CDS 和启动子。gene 注释范围与 CDS 范围可能不同；启动子依赖转录起始位点和任务定义。细菌 CDS 区间与非 CDS 区间没有被拿来充当 GUE 启动子正负样本。

6. 先划分来源，再切成样本

若对同一染色体滑窗后随机打散，相邻样本可能共享大量碱基，测试表现会失真。这里先在原始坐标上划分，再在每个集合内切窗。

大肠杆菌采用连续区段约 80%/10%/10%，边界两侧各留 512 bp 隔离。酵母把核染色体 I 留给验证、II 留给测试，其他核染色体用于训练，线粒体不进入该教学训练。

默认窗口长度和步长都为 512。不能凑足一个窗口的末尾短片段丢弃并在区间与输出数量中核对；PAD 用于独立的批次演示。N 比例大于 10% 时过滤整个窗口，序列内部不删字符。10% 是本教程可修改的配置选择。

Get-Content-LiteralPath'data\processed\ecoli\regions_512.tsv'Get-Content-LiteralPath'data\processed\ecoli\windows_512.tsv'-TotalCount 3

本次 512 bp 数据集为训练7,251、验证904、测试905条，总计9,060条；256 bp 数据另存，供原核阶段 MLM 使用。

• 512 bp 是单条数据（单张“试卷/片段”）的长度，也就是模型的窗口大小（Window Size）。
• 9,060 条 是样本的数量（总行数）。
• MLM 的全称是 Masked Language Modeling（掩码语言模型）。

每条样本包含来源版本、seqid、起止位置、方向与 split，可以用reference[start0:end0]原路核对。脚本同时清除跨集合 exact/RC 等价窗口。

7. 构建逐碱基标签

标签文件cds_windows_512.npz的主要数组：

&$py-c"import numpy as np; d=np.load('data/processed/ecoli/cds_windows_512.npz',allow_pickle=False); print(d.files); print(d['input_ids'].shape); print(d['labels'].shape); print(d['labels'][0,:30].tolist())"

两个方向可以同时标为 1。本次整个参考序列上有2,994 个碱基被相反方向的 CDS 注释覆盖；单一三分类标签会丢失这种重叠关系，所以使用两个通道。标签 0 表示该方向无 CDS 注释。

进一步逐碱基追踪一个实际样本：

&$pytools/trace_sample.py--index 0

打开 单个样本逐步追踪，逐行对照参考坐标、窗口内坐标、CDS注释与0/1标签。用--index换一个样本即可重复检查。

窗口截断的 CDS 仍可提供窗口内逐碱基标签，不能据此把该片段宣称为完整基因。带 partial、pseudo、exception 等标记的相关位置保留标记，并使用监督掩码区别处理。

&$py-m dna_learning report

打开reports/cds_labels_20kb.png：上方是注释区间，下方是两条逐碱基标签。逐段对照边界，是验证标签是否偏移的直接办法。

反向互补一个窗口时，需要把序列与标签同时变换。长度 L 的区间[s,e)变为[L-e,L-s)，标签倒序并交换正负链。代码在encoding.reverse_cds_labels。

8. 检查真实启动子分类数据

&$py-m dna_learning prepare--only gueGet-Content-LiteralPath'data\processed\gue\audit.json'

原始训练/验证/测试集为4,904 / 613 / 613条，每条 300 bp。标签 1 表示该基准的启动子正例，0 表示按基准规则构建的负例。

本次清理规则预先固定：

1. 保留原始文件与 dev/test。
2. 训练副本隔离与留出集 exact/RC 重复的记录：匹配 dev 的 14 条、匹配 test 的 6 条。
3. 仅根据训练集内部的标签，发现 10 个冲突组，共 20 条训练记录；全部隔离，不替它们猜测正确标签。

最终训练集为4,864 条，验证、测试仍各613 条。跨全部集合一共发现 15 个 canonical 序列标签冲突组；这一全局统计只用于数据审计，训练冲突判断只读取训练标签。

Get-Content-LiteralPath'data\processed\gue\removed_train.tsv'-TotalCount 5Get-Content-LiteralPath'data\processed\gue\distributions.tsv'

removed_train.tsv逐条解释删除原因。缺少原始坐标的 GUE 样本把坐标字段留空，并明确来源，而不猜测它来自哪个基因。

继续检查近重复：

&$py-m dna_learning similarityGet-Content-LiteralPath'reports\similarity\near_pairs.tsv'

这一步使用等长无缺口Hamming距离，检查300bp内最多6个替换，并比较RC方向。本次exact/RC清理后仍找到1对训练—测试近重复，相差4个碱基。它单独记录在报告中，说明exact去重尚不能代表完全独立。本教程保留预先定义的清理规则；若建立自己的严格基准，应先确定相似度分组规则，再做训练，不能根据测试成绩反复调规则。这个检查不覆盖indel或全部同源关系。

查看图片比较三组 GC 分布。长度、GC、N 和类别分布都可能构成模型可利用的信号或偏差。

9. 从字符到模型输入

&$py-m dna_learning preview

输出保存到reports/batch_preview.json。输入编码固定为：

PAD=0 MASK=1 A=2 C=3 G=4 T=5 N=6

MLM 的输出类别是A=0,C=1,G=2,T=3,N=4。输入 token 与预测类别是两套编号；-100是“不参与损失”的位置。

预览批次形状为[2,300]，有效长度为[300,292]。第二条仅在此演示里缩短 8 bp，右侧补 8 个 PAD，对应 attention mask 为 false。分类标签形状是[2]；MLM 标签形状是[2,300]。

&$py-c"import numpy as np; from dna_learning.encoding import batch_encode; ids,valid,lengths=batch_encode(['ACGTN','AG']); onehot=np.eye(5)[ids[0]-2]; print(ids); print(valid); print(onehot)"

one-hot 展示每个碱基的一组互斥指示值；模型实际使用可训练 embedding。处理大数据时按需取批次，不把整个基因组展开为巨大的 one-hot 张量。

MLM 随机选择约 15% 的位置，其中 80% 替换 MASK、10% 替换随机碱基、10% 保留原值。短序列至少选择一个目标，这是教学实现选择。反向互补输入由已遮盖的正链生成，避免另一条链泄漏答案。

10. 先做一个可解释的序列分析基线

&$py-m dna_learning baseline

该实验使用 64 个三碱基组合的频率，再加 GC 和 N 两项特征，用逻辑回归分类，共 66 个特征。含 N 的 k-mer 跳过，频率按有效 k-mer 数归一化。

模型只在清理后的训练集拟合，参数固定为C=1.0、种子 42。结果在reports/kmer_baseline/metrics.json，逐样本预测与系数另存。系数能帮助理解哪些特征与本数据标签关联；这种关联本身不构成功能机制证明。

这个基线提供一个很快能跑通的完整闭环，也让你在深度学习之前观察数据里已有的分类信号。