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

DNA 序列数据分析与处理逐步操作手册你要完成的工作是把一个序列文件变成来源清楚、坐标正确、标签可信、能交给模型的数据集。阅读顺序按数据流组织。前九章用 CPU 即可完成。本项目有两条实践线。大肠杆菌基因组用于学习“序列如何与注释、方向、逐碱基标签对应”人类启动子数据用于学习“短序列分类数据如何检查、去重、编码和评估”。两套任务的标签分别构建。0. 运行位置与环境打开 PowerShell复制以下命令。Set-Location-LiteralPathG:\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-envGPU 模型章另装 PyTorch数据处理阶段无需等待模型依赖。本次环境torch固定为Windows/Python3.12的CUDA12.6构建迁移到其他系统时应按PyTorch官方兼容表选择对应构建数据处理依赖可独立安装。1. 先学会辨认文件1.1 生成可反复练习的小文件$py-m dna_learning examplesGet-Content-LiteralPathdata\examples\toy.fastaGet-Content-LiteralPathdata\examples\toy.fastqFASTA 示例为toy_chr 人工教学序列 acgtACGTNNRYACGTACGT后第一项是序列 ID其余内容是描述。换行只用于排版不能把每行当作独立序列。一个文件可包含多条以开头的记录。FASTQ 示例为toy_read_1 artificial_Phred33 ACGTNACG IIII!5?IFASTQ格式是一种基于ASCII编码的文本格式主要用于存储生物序列及其测序质量评价数据为高通量测序领域的标准格式并被NCBI Short Read Archive等数据库采用。该格式包含四行文本首行以开头记录序列标识符及描述信息第二行为碱基序列第三行以起始可附带注释第四行通过特定ASCII字符表示每个碱基的测序质量评分。 [1]FASTQ的质量评分采用Phred算法将错误概率转换为ASCII字符如Q30对应字符’?。Illumina测序仪生成的FASTQ文件包含flowcell、lane、tile等设备位置参数。该格式可通过Biopython、EMBOSS等工具与FASTA等格式互相转换适用于Hiseq系列测序平台的数据存储与分析。这里采用 Phred33。I对应 Q40!对应 Q05对应 Q20?对应 Q30。质量描述的是测序读段中每个碱基的可信度序列与质量字符数必须相同。一般 FASTQ 可有换行正式读取使用解析器不能一律按四行手工分组。在二代/高通量测序FASTQ 格式中Q40、Q30、Q20 这里的 Q 指的是 Phred 质量得分Phred Quality Score它是用来量化测序仪对某个碱基识别的“准确率”或“测错概率”的指标。$py-m dna_learning inspectdata/examples/toy.fastq预期format为fastqrecord_count为 1长度为 8平均 Phred 为 31.25。这是格式教学样例尚未执行测序读段去接头或比对。1.2 同时认识其他文件的职责文件保存什么首先检查什么FASTA /.fa/.fna序列 ID 与碱基序列ID、记录数、长度、字符FASTQ /.fq读段序列与质量分数格式、长度配对、质量编码GenBank /.gbff序列、元信息、基因等特征版本、feature、位置与方向GFF3注释特征与坐标seqid、type、start/end、strand、phase、属性BED基因组区间染色体、起止位置、方向VCF相对参考序列的变异参考版本、POS、REF、ALT、样本字段分类 CSV已提取的序列与任务标签标签含义、来源、划分和重复$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-LiteralPathdata\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。这些数值是本次真实文件的检查结果。逐项理解以下指标。指标计算或含义解释长度记录中碱基字符数检查是否读取完整、是否误读了 CDS 文件GC 比例(GC)/(ACGT)本教程的分母排除 N 和其他模糊碱基N 比例N/序列长度描述未知碱基的占比模糊碱基比例非 A/C/G/T 字符占比包含 N、R、Y 等小写数量原始小写字符数可能承载软屏蔽信息应在大写化前记录重复相同序列或反向互补等价序列用于后续划分审计核心代码在dna_learning/io.py的normalize与base_stats。示例$py-cfrom dna_learning.io import normalize,base_stats; s,changesnormalize(aCGtNRY); print(s); print(changes); print(base_stats(s))预期规范化序列为ACGTNNN长度仍为 7R/Y 的位置为[5, 6]。原始文件不改动。若直接把 N 删除后续所有注释坐标都会发生偏移。本项目把合法模糊碱基映射为 N-、?等非法字符会报错提示检查数据来源或格式。GC 异常用于提示进一步调查单凭这一项不能确定污染。4. 理解坐标最容易造成静默错误的一步本项目内部统一采用 Python 风格的0 基、左闭右开[start0,end0)。GFF3 的原始起止坐标采用 1 基闭区间转换为start0start1-1、end0end1。序列 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-csACGTAACTGGCA; 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↔GN 保持。$py-cfrom dna_learning.io import read_fasta,reverse_complement; pnext(iter(read_fasta(reports/first_300_plus.fasta).values())); mnext(iter(read_fasta(reports/first_300_minus.fasta).values())); print(mreverse_complement(p)); print(reverse_complement(m)p)两个结果都应为True。5. 用注释提取 CDS并反查生物学含义CDS 是基因组中真正能够被细胞“翻译”成蛋白质氨基酸的那一段连续的核苷酸序列。$py-m dna_learning prepare--only genomes产物位于data/processed/ecoli文件用途features.tsv/features.json每条 CDS 的来源、分段坐标、方向、蛋白 ID 与例外cds.fasta按注释方向提取的 CDS 序列translation_checks.tsvCDS 翻译与原注释的核对结果coordinate_examples.json真实 GFF3 坐标到 Python 坐标的示例本次提取4,318 条 CDS。4,297 条完整且无相关例外标记的 CDS 翻译匹配21 条带例外标记而跳过直接核对。跳过原因可从表中追溯。Get-Content-LiteralPathdata\processed\ecoli\features.tsv-TotalCount 4Get-Content-LiteralPathdata\processed\ecoli\translation_checks.tsv-TotalCount 8提取核心使用feature.extract(record.seq)让 Biopython 处理负链和复合位置。翻译根据注释的transl_table与codon_start完整 CDS 使用cdsTrue从而正确处理细菌的替代起始密码子。这里区分 gene、CDS 和启动子。gene 注释范围与 CDS 范围可能不同启动子依赖转录起始位点和任务定义。细菌 CDS 区间与非 CDS 区间没有被拿来充当 GUE 启动子正负样本。6. 先划分来源再切成样本若对同一染色体滑窗后随机打散相邻样本可能共享大量碱基测试表现会失真。这里先在原始坐标上划分再在每个集合内切窗。大肠杆菌采用连续区段约 80%/10%/10%边界两侧各留 512 bp 隔离。酵母把核染色体 I 留给验证、II 留给测试其他核染色体用于训练线粒体不进入该教学训练。默认窗口长度和步长都为 512。不能凑足一个窗口的末尾短片段丢弃并在区间与输出数量中核对PAD 用于独立的批次演示。N 比例大于 10% 时过滤整个窗口序列内部不删字符。10% 是本教程可修改的配置选择。Get-Content-LiteralPathdata\processed\ecoli\regions_512.tsvGet-Content-LiteralPathdata\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的主要数组数组本次形状含义input_ids[9060,512]每个碱基的 tokenlabels[9060,512,2]最后两通道依次为正链和负链 CDSsupervision_mask[9060,512,2]该位置、该方向是否用于可靠监督attention_mask[9060,512]输入是否为真实序列位置sample_id/split[9060]与窗口来源表关联$py-cimport numpy as np; dnp.load(data/processed/ecoli/cds_windows_512.npz,allow_pickleFalse); 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-LiteralPathdata\processed\gue\audit.json原始训练/验证/测试集为4,904 / 613 / 613条每条 300 bp。标签 1 表示该基准的启动子正例0 表示按基准规则构建的负例。本次清理规则预先固定保留原始文件与 dev/test。训练副本隔离与留出集 exact/RC 重复的记录匹配 dev 的 14 条、匹配 test 的 6 条。仅根据训练集内部的标签发现 10 个冲突组共 20 条训练记录全部隔离不替它们猜测正确标签。最终训练集为4,864 条验证、测试仍各613 条。跨全部集合一共发现 15 个 canonical 序列标签冲突组这一全局统计只用于数据审计训练冲突判断只读取训练标签。Get-Content-LiteralPathdata\processed\gue\removed_train.tsv-TotalCount 5Get-Content-LiteralPathdata\processed\gue\distributions.tsvremoved_train.tsv逐条解释删除原因。缺少原始坐标的 GUE 样本把坐标字段留空并明确来源而不猜测它来自哪个基因。继续检查近重复$py-m dna_learning similarityGet-Content-LiteralPathreports\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。输入编码固定为PAD0 MASK1 A2 C3 G4 T5 N6MLM 的输出类别是A0,C1,G2,T3,N4。输入 token 与预测类别是两套编号-100是“不参与损失”的位置。预览批次形状为[2,300]有效长度为[300,292]。第二条仅在此演示里缩短 8 bp右侧补 8 个 PAD对应 attention mask 为 false。分类标签形状是[2]MLM 标签形状是[2,300]。$py-cimport numpy as np; from dna_learning.encoding import batch_encode; ids,valid,lengthsbatch_encode([ACGTN,AG]); onehotnp.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 数归一化。模型只在清理后的训练集拟合参数固定为C1.0、种子 42。结果在reports/kmer_baseline/metrics.json逐样本预测与系数另存。系数能帮助理解哪些特征与本数据标签关联这种关联本身不构成功能机制证明。这个基线提供一个很快能跑通的完整闭环也让你在深度学习之前观察数据里已有的分类信号。