已有全基因组数据，还能回答什么新问题？一次线粒体数据复用研究的启示

一批全基因组测序数据已经完成 SNP 检测或群体分析，原始 FASTQ 随后进入归档。这时，如果研究问题变成“样本之间的母系谱系有什么差异”，或者“现有标记是否足以区分材料”，下一步未必需要立即重新送样。

更值得先做的是数据可行性评估：旧数据中有没有足够的目标信号，能否形成可信的新结果，以及继续分析的成本是否合理。 本文以一项蛔虫线粒体研究为起点，讨论可复用的方法，给出原创代码示例，并说明哪些结论需要额外证据。

这项研究提供了什么证据

Woolfe 等在 Scientific Reports 发表的研究，对 149 份蛔虫 WGS 数据比较四种处理路线，再进行线粒体组装。按文中“单条 contig 含全部 12 个蛋白编码基因”的口径，结果如下：研究原文。

这些百分比由原文计数换算，不是通用成功率。研究识别出新的线粒体分支 D；但并非 149 份数据都自动得到无缺口结果，其中一份经过人工处理仍有缺口。原始 reads 回贴检查覆盖的是其中 30 份。论文也没有给出适用于其他项目的固定节省比例。方法与结果说明

从“再出一张图”转向一个明确的新问题

数据复用的价值取决于它能补上哪一段证据。针对已经持有动物或寄生虫 WGS 的团队，可以先把问题拆成以下几类；这是项目设计建议，不代表原论文逐一验证了这些应用。

这些分析不能直接替代物种界定、真实传播链追踪或临床诊断。一个线粒体聚类首先说明序列关系；把它解释为宿主间传播、杂交或独立物种，还需要核基因组、流行病学或其他独立证据。

可复用的技术亮点：把流程变成可以比较的方案

亮点一：先核查输入，避免对着结果表寻找已经丢失的序列。 MitoZ 的方法研究表明，动物全基因组 shotgun reads 可以用于线粒体组装，但公共数据也可能事先移除了线粒体 reads。因此，“做过 WGS”不等于“现在交到手上的文件仍保留目标信息”。只有 SNP 表或过滤后的 VCF，通常不能替代原始 reads 完成这里讨论的组装任务。MitoZ 方法论文

亮点二：保留直接分析的基线，检查额外步骤是否真正有用。 去宿主、目标富集和参数调整都增加了处理环节。可以对同一组代表性数据比较候选路线，而不是先假定复杂方案一定更优。作者公开仓库中的预处理、去宿主、组装与绘图代码可以作为流程核对入口；具体参数仍需与正式论文及输入读长对应。作者代码仓库

亮点三：把“程序完成”与“结果合格”分开。 我们建议至少保留以下四层验收记录：

1. 数据层： 读长、质量、有效数据量及样本标识是否可靠。
2. 序列层： 候选长度、未知碱基、异常重复和预期基因集合是否合理。
3. 支持层： reads 是否覆盖候选序列，低覆盖区域、接头附近和冲突位点是否需要检查。
4. 解释层： 不同分析口径是否给出一致结论；元数据是否足以支持生物学解释。

这套验收表比“某个软件运行成功”更适合变成可交付的服务。它也方便把不同路线的结果放在同一尺度上比较。

代码示例：检查一个已经组装出的候选序列

下面给出一段独立编写的 Linux/Bash 示例，使用双端短读长 FASTQ 和已经得到的候选线粒体 FASTA。它完成质控、回贴和覆盖度汇总，不包含组装步骤，也不是原论文的完整复现。

先将 mt_depth_summary.py 放在工作目录，准备 sample.R1.fq.gz、sample.R2.fq.gz 和 candidate.fa。需要已安装 fastp、Bowtie2、Python 3，以及 SAMtools 1.13 或更高版本；实际项目应固定并记录具体版本。为了便于阅读，下方展示主要命令；完整代码包中的 check_candidate.sh 另有输入和版本检查。

set -euo pipefail

# 使用一个不存在的新目录，避免覆盖之前的结果。
out=mt_pilot
mkdir "$out"

fastp -i sample.R1.fq.gz -I sample.R2.fq.gz \
  -o "$out/clean.R1.fq.gz" -O "$out/clean.R2.fq.gz" \
  --detect_adapter_for_pe --thread 4 \
  --html "$out/fastp.html" --json "$out/fastp.json"

cp candidate.fa "$out/candidate.fa"
samtools faidx "$out/candidate.fa"
bowtie2-build "$out/candidate.fa" "$out/mt_index"

bowtie2 --very-sensitive -p 4 -x "$out/mt_index" \
  -1 "$out/clean.R1.fq.gz" -2 "$out/clean.R2.fq.gz" \
  2>"$out/map.log" \
  | samtools view -u -F 2308 - \
  | samtools sort -@ 2 -o "$out/mapped.bam" -

samtools index "$out/mapped.bam"
samtools depth -aa -s -q 20 -Q 20 \
  "$out/mapped.bam" >"$out/depth.tsv"

python3 mt_depth_summary.py \
  "$out/candidate.fa.fai" "$out/depth.tsv" \
  >"$out/coverage_summary.tsv"

fastp 的双端输入、接头检测及 HTML/JSON 报告参数参见官方说明；Bowtie2 的 --very-sensitive 是比对敏感性预设，并不意味着结果一定真实。Bowtie2 手册

这里有几个影响结果口径的细节：

• -F 2308 排除未比对、次级和补充比对记录；覆盖统计仍会受到比对软件如何选择主比对的影响。
• depth -aa 保留所有参考位置，包含完全没有 reads 的序列。如果只统计有覆盖的位置，平均深度和覆盖广度可能被高估。
• -s 减少同一对 reads 重叠区的重复计数；它不是 PCR 去重。
• 在 samtools depth 中，-q 是最低碱基质量，-Q 是最低比对质量。示例中的 20 是演示参数，需要按数据调整；不要与 mpileup 的同名字母参数混用。SAMtools depth 手册

用 Python 检查分母，而不只是求平均数

下面是代码包中的完整覆盖统计脚本。它要求单个 BAM 的三列 depth -aa 输出，并用 FASTA 索引核对每个 contig 的位置数。缺少零覆盖位置、重复位置或意外输入多个 BAM 时会报错，不会悄悄给出偏高结果。

#!/usr/bin/env python3
"""Summarize ONE BAM's samtools depth -aa output; not an assembly validator."""
import argparse
import csv
import sys


def summarize(fai_path, depth_path):
    lengths = {}
    with open(fai_path, encoding="utf-8") as handle:
        for line in handle:
            fields = line.rstrip("\n").split("\t")
            if len(fields) < 2:
                raise ValueError("Invalid FASTA index row")
            name, length = fields[0], int(fields[1])
            if not name or name in lengths or length <= 0:
                raise ValueError("Empty/duplicate reference name or invalid length")
            lengths[name] = length
    if not lengths:
        raise ValueError("Empty FASTA index")
    # count, depth sum, bases >=1x, bases >=10x
    stats = {name: [0, 0, 0, 0] for name in lengths}
    with open(depth_path, encoding="utf-8") as handle:
        for line in handle:
            fields = line.rstrip("\n").split("\t")
            if len(fields) != 3:
                raise ValueError("Expected CHROM, POS, DEPTH for exactly one BAM")
            name, pos, depth = fields[0], int(fields[1]), int(fields[2])
            if name not in stats:
                raise ValueError("Reference absent from FASTA index: " + name)
            row = stats[name]
            if pos != row[0] + 1 or pos > lengths[name] or depth < 0:
                raise ValueError("Missing, repeated, unordered, or invalid position")
            row[0] += 1
            row[1] += depth
            row[2] += depth >= 1
            row[3] += depth >= 10
    result = []
    for name, length in lengths.items():
        count, total, covered, covered10 = stats[name]
        if count != length:
            raise ValueError("Incomplete positions: use depth -aa without a region filter")
        result.append((name, length, total / length,
                       100 * covered / length, 100 * covered10 / length))
    return result


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("fai", help="Index of exactly the candidate reference used for mapping")
    parser.add_argument("depth", help="Three-column depth -aa output for one BAM")
    args = parser.parse_args()
    try:
        rows = summarize(args.fai, args.depth)
    except (OSError, ValueError) as exc:
        parser.error(str(exc))
    writer = csv.writer(sys.stdout, delimiter="\t", lineterminator="\n")
    writer.writerow(["contig", "length_bp", "mean_depth", "breadth_1x_pct", "breadth_10x_pct"])
    for name, length, mean, one, ten in rows:
        writer.writerow([name, length, f"{mean:.3f}", f"{one:.2f}", f"{ten:.2f}"])


if __name__ == "__main__":
    main()

例如，一条长度为 4 bp 的演示序列，逐位点深度是 0、10、20、0，平均深度应为 7.5×，至少 1× 和至少 10× 的覆盖广度都为 50%。如果把两个零值丢掉，就会得到误导性的 15× 均深和 100% 覆盖。这里的演示序列是合成数据，不是研究样本。

脚本只输出统计量，不给出“组装通过”的自动判定。10× 是一个汇总观察点，不是本项目或所有物种的合格阈值。 随包测试覆盖零值分母、完全无覆盖的 contig、缺失位置、重复位置和多 BAM 误输入。上述五项合成数据测试已通过；本文没有对原论文 149 份样本做独立复现，也没有宣称已完成 fastp—Bowtie2 全流程的真实样本验证。

不足与边界：哪些地方仍需要专业判断

候选序列能被 reads 覆盖，不等于它一定来自真实线粒体。 核基因组中的线粒体来源片段（NUMTs）和污染序列都可能干扰判断。MitoZ 的方法设计也包含对这类假阳性的处理。上面的示例只向候选 FASTA 回贴，不能独立排除 NUMTs；必要时需要联合核参考、候选结构、基因完整性和其他证据复核。MitoZ 方法论文

一个高均深不能遮住局部缺口。 应查看逐位点覆盖、低质量区和冲突位点。线性表示的候选序列两端还可能受比对边界影响；要判断环化，需专门检查跨接头支持，不能仅因为 FASTA 只有一条序列就称其为闭合环状基因组。

目标富集可能改变你看见的多样性。 如果过滤库离真实谱系较远，严格匹配可能丢失有用 reads。这是选择富集策略时需要检查的潜在偏差，可以通过保留未富集基线和复核异常样本来评估，不应在未比较前认定某条路线更好。

动物项目的参数不能直接搬到植物或其他数据类型。 建库方式、基因组结构、注释数据库和遗传密码表都应与研究对象匹配。MitoZ 面向动物线粒体，迁移时应核对具体版本的支持范围。MitoZ 项目说明

复用数据仍然有成本。 算力、人工核查、失败重试和补充实验都应计入。对于路线选择，一个实用指标是“新增总成本 ÷ 合格结果数”，并确保所有路线使用相同验收标准。没有合格结果时，应报告试做未通过，而不是给出虚假的单位成本优势。

如何把它做成一个可交付的分析项目

建议先从小范围评估开始，再决定是否批量推进。以下是本类需求可以讨论的交付范围，不是已经完成的客户案例：

如果你已经有 WGS 数据，希望判断是否适合继续挖掘线粒体或其他生物学信息，可以通过云生信定制分析服务提交需求。首次沟通先准备四项信息：研究对象与样本数、测序平台和建库方式、现有文件类型、最希望解决的一个问题。先提交数据概况即可，原始数据传输方式和分析范围在沟通后确认。

云生信提供需求评估、分析方案与交付范围沟通；具体可行性、周期和费用以项目评估为准。能否重建目标序列、得到什么结论，需要由你的数据决定。服务说明

原文、工具与代码

• Woolfe L, et al. Scalable assembly of Ascaris mitogenomes from whole-genome data reveals a novel clade. Scientific Reports, 2026. DOI: 10.1038/s41598-026-65562-w。
• 研究作者的流程代码。
• Meng G, et al. MitoZ: a toolkit for animal mitochondrial genome assembly, annotation and visualization. Nucleic Acids Research, 2019. DOI: 10.1093/nar/gkz173。
• fastp 官方文档、Bowtie2 官方手册、SAMtools depth 官方手册。
• 下载本文原创代码与合成数据测试。代码用于教学和方案试做，不替代针对项目的数据验证。

本文为文献技术解读与工程示例。论文结果、本文建议和示例测试已分别标明；文中未使用原论文图片，也未将研究团队结果表述为云生信已完成的客户成果。