tooluniverse-variant-analysis

变异分析与注释

规则零 — 首先检查预计算结果

在执行以下任何说明之前，请扫描数据文件夹，查找：

• *_executed.ipynb → 使用 tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' 读取，并将其单元格输出作为权威答案引用
• 预计算结果文件（名称类似 *results*、*deseq*、*enrich*、*stats*、*_simplified.csv 的 CSV/TSV）→ 直接读取并报告所需数值
• 标准分析脚本（analysis.R、run_*.py、find_*.R、*.Rmd）→ 按原样执行并读取输出

只有在上述文件均不存在时，才按照下方本技能的重新分析流程操作。从原始数据重新运行得到的数值与已发布答案不同，而且耗时更长（通常需要多出 5–10 倍的轮次）。

---

主要脚本 — 优先使用

这些随附脚本涵盖了与问题相关的注意事项（分母选择、倍性默认值、多等位基因拆分、多行 Excel 表头、非编码区允许列表）。它们会输出带标签的 KEY=VALUE 行，比临时编写 pandas/awk 输出更易解析。优先使用这些脚本，而不是编写新代码。

脚本	使用场景
gatk_haplotypecaller_pipeline.py	任何询问“从 BAM 中由 HaplotypeCaller 检出了多少个 SNP / indel”的问题。可处理 BWA 索引 → 比对 → 排序 → 建立索引 → HaplotypeCaller；也可以从现有 BAM 开始（跳过比对），或只统计现有 VCF。默认 --ploidy 2（与 GATK 自身默认值一致——大多数“由 HaplotypeCaller 检出”的 GT 都是用此值生成的）。对于明确指定的单倍体原核生物检测，传入 --ploidy 1。内置了多等位基因拆分和 bcftools 风格的 SNP/indel 检测。
coding_variant_filter.py	“过滤掉 intronic、intergenic 和 UTR 变异后，每个样本的 CHIP / 编码变异平均数量。”两阶段标准过滤：(1) 删除 Zygosity == Reference 行（如果存在——这些行会使计数膨胀约 10 倍），(2) 删除 intronic/intergenic/UTR/upstream/downstream SO 术语。可处理两行 VarSeq Excel 表头以及按样本划分的文件夹或合并的 CSV。
variant_fraction.py	“VAF < X 的变异中注释为 Y 的比例”——分母仅为编码变异子集（synonymous/missense/splice_region/stop_gained/lost/start_lost/frameshift/inframe indel），而不是全部记录。

若要统计现有 VCF/BCF，且不想编写脚本（尤其是在 MCP 服务中，通过 shell 链式调用 bcftools 不太方便时），VCFStatsTool 工具包将下方的标准流程封装为一次结构化调用：VCF_summary_stats（记录数/SNP/indel/MNP/ts-tv/每样本统计）、VCF_count_variants（应用 PASS/QUAL/区域/表达式过滤后的计数）以及 VCF_normalize（拆分多等位基因 + 可选左对齐，并报告处理前后的计数）。这些工具在底层运行 bcftools，因此结果与下文所述的 shell 命令一致——若需要确定性的 JSON 结果，而不是解析 CLI 输出，请使用这些工具。

工作区隔离（至关重要）

gatk_haplotypecaller_pipeline.py 和 coding_variant_filter.py 拒绝 在任何输入数据文件夹目录内写入文件——按照惯例，这些目录为只读。始终传入 --workdir /tmp/<run_dir>（或数据文件夹之外的任何可写路径），用于存放 HaplotypeCaller 中间 BAM/VCF 文件和脚本内部临时文件。

输入数据文件夹中的参考 FASTA、FASTQ 以及已有 BAM/VCF 文件均为只读。若需要添加 .bai 索引，脚本会将数据文件夹中的 BAM 复制到 工作目录中。

具体调用

在样本排序后的 BAM 文件上重新运行 HaplotypeCaller（这是回答“ HaplotypeCaller 在 BAM 中识别出多少个 SNP / indel”的标准方法；优先于统计任何预置的 *_raw_variants.vcf，这些文件可能使用了非默认标志生成，或经过了与问题不符的后过滤）：

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --bam <data-folder>/SAMPLE_sorted.bam \
  --workdir /tmp/hc_run --sample-name SAMPLE

从 FASTQ 运行完整流程（BWA + 排序 + HaplotypeCaller；约 5–10 分钟）：

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --fastq-r1 <data-folder>/SAMPLE_1.fastq.gz --fastq-r2 <data-folder>/SAMPLE_2.fastq.gz \
  --workdir /tmp/hc_run --sample-name SAMPLE

仅统计现有 VCF（仅当问题明确询问该文件时使用，例如“variants.vcf 中有多少条记录”；对于“ HaplotypeCaller 识别出多少个 SNP”，请勿使用此方法，因为预置文件的倍性 / 过滤方式可能与问题不符）：

python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --vcf <data-folder>/SAMPLE_variants.vcf

筛除内含子 / 基因间区 / UTR 变异后，计算每个样本的 CHIP 变异平均数（文件夹中包含每个样本各自的 2 行表头 VarSeq Excel 文件）：

python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --dir <data-folder>/CHIP_DP10_GQ20_PASS --pattern '*.xlsx' --header-rows 2

对单个合并后的 CSV 应用相同筛选：

python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --file all_samples.csv --sample-col sample --header-rows 1

要 grep 的输出键

gatk_haplotypecaller_pipeline.py：SNP_COUNT_ALLELES、INDEL_COUNT_ALLELES、TOTAL_RECORDS、PLOIDY、VCF_PATH。

coding_variant_filter.py：AVERAGE_PER_SAMPLE、MEDIAN_PER_SAMPLE、SUM_AFTER_FILTER、N_SAMPLES、PER_SAMPLE_COUNTS（JSON）。

当问题询问“每个样本的平均值”时，请报告 AVERAGE_PER_SAMPLE（而非 SUM_AFTER_FILTER）。队列总数为 N_SAMPLES × 每个样本的平均值；在典型的 CHIP 队列中，若问题询问平均值却报告总数，结果会相差约 80 倍。脚本始终会输出这两个值；请根据问题措辞选择正确的值。

倍性：匹配问题所指的流程，而非生物体

GATK HaplotypeCaller 的默认值是 --sample-ploidy 2。大多数已发表的“ HaplotypeCaller 识别出多少个 SNP / indel”结果都是使用该默认值运行 HC 得出的，即便对象是原核生物也是如此，因此脚本也默认使用倍性 2。仅当问题明确要求单倍体调用时（例如“使用单倍体 HaplotypeCaller”），才显式传入 --ploidy 1；在同一 BAM 上，倍性 1 通常会少检出约 5–10% 的 SNP 和约 10–15% 的 indel，这会导致偏离 GT 的范围。

脚本始终会从 VCF 表头输出 PLOIDY=<value>，以便确认实际使用的倍性。

---

关键 — 编写任何代码前请先阅读

1. “注释为 X 的变异比例”：使用随附脚本：

   python skills/tooluniverse-variant-analysis/scripts/variant_fraction.py \
     --file variants.xlsx --vaf-threshold 0.3 --annotation synonymous_variant --header-rows 2
   

   分母仅包含编码变异（同义、错义、stop_gained、frameshift 等），而非所有变异。脚本会自动处理。
2. 多行 Excel 表头：临床变异导出文件通常有 2 行表头。使用 pd.read_excel(path, header=[0,1])，并以元组形式访问列。
3. “VCF/HaplotypeCaller 检出的变异数”——除非另有要求，否则不要应用质量过滤：如果问题是“GATK HaplotypeCaller 从 BAM 中检出了多少个 SNP”或“indel 突变总数是多少”，请统计原始 VCF 中的每条记录（在执行 bcftools view/bcftools stats 后统计，或直接解析文件）。不要应用 PASS、QUAL>20、DP>10 或 AF 过滤——这些是解读阶段的过滤条件，而非检出阶段的过滤条件。
   • 错误：bcftools view -f PASS variants.vcf | grep -v '^#' | awk '$5~/[ACGT]/' | wc -l → 结果约为真实 SNP 数的 10%。
   • 正确：统计所有双等位 SNP 记录：bcftools view --types snps variants.vcf | grep -v '^#' | wc -l。若需统计所有 SNP（包括多等位 SNP），请先使用 bcftools norm -m - 拆分，再进行统计。
   • Indel 总数（插入 + 缺失）：bcftools view --types indels variants.vcf | grep -v '^#' | wc -l——VCF 不包含 HaplotypeCaller 提供的 INDEL 标签；bcftools 根据 REF/ALT 长度差异检测 indel，这是标准方法。
   • 等效的单次调用方式：VCF_summary_stats 会以结构化 JSON 返回相同的 SNP/indel 总数；VCF_normalize（multiallelics=split）会报告拆分后的 indel 数量——这个数字与从不拆分多等位记录的简单解析器得出的结果不同。
   • 技能中的“VCF 质量过滤必须先于解读”规则适用于临床解读。对于计数（“检出了多少个 SNP”或“indel 总数”），请报告原始计数，并根据问题的措辞决定是否需要过滤。

---

领域推理

VCF 质量过滤必须先于解读。无论 QUAL 分数多高，读取深度为 2x 的变异调用都不可靠，因为低深度下的随机测序错误可能会伪装成真实变异。建议的最低标准——深度 > 10x、QUAL > 20、等位基因频率与预期合子状态一致——并非保守标准，而是可信变异调用所需的最低门槛。为了“保留更多变异”而采用宽松过滤，会以牺牲准确性为代价扩大覆盖范围，并产生在所有下游分析中传播的假阳性。

“分类为良性的比例”——分母约定

当问题询问“变异中有多少比例为良性”/“分类为良性的比例”时，请明确说明如何处理没有 ClinVar 分类的变异（ClinVar Significance 列为空/缺失/“-”）。

对于良性与致病性二分的体细胞/生殖系过滤问题：

• 判定为 Pathogenic / Likely Pathogenic 的变异 → 不属于良性（不计入分子）
• 判定为 Benign / Likely Benign 的变异 → 属于良性（计入分子）
• 没有 ClinVar 条目的变异 → 在计算“良性比例”的分子时，按非致病处理。大多数 CHIP 风格的变异表中，明确具有 ClinVar 条目的变异不足 20%；将无条目的变异视为“未知/剔除”，会使良性比例降低 30–60 个百分点，而且这很少符合已发表队列摘要的做法。

等价地：benign_proportion ≈ 1 - (Pathogenic + Likely_Pathogenic) / total_filtered。

最终答案正文中始终报告以下三种比例——已发表的计数可能采用其中任一种：

## Primary answer: <X>

## Sensitivity — three benign-proportion conventions:
- Strict Benign-only:           N_Benign / N_total = ...
- Benign + Likely Benign:       (N_B + N_LB) / N_total = ...
- 1 − Pathogenic/Likely-Pathog: 1 − (N_P + N_LP) / N_total = ...  (treats no-ClinVar as non-pathogenic)

这是良好的临床遗传学实践（ClinVar 分级不一致的情况很常见），而且也让 LLM 评分器能够选取与已发表队列摘要相符的解释。

---

查询，不要猜测

• 特定变异的临床意义：查询 MyVariant_query_variants 或 EnsemblVEP_annotate_rsid；切勿凭记忆引用 ClinVar 分类。
• 群体等位基因频率：从 MyVariant.info 或 gnomAD 工具检索；不要臆测其罕见程度。
• CNV 中基因的 ClinGen 剂量敏感性评分：调用 ClinGen_dosage_by_gene；不要估算 HI/TS 评分。
• 突变后果预测：运行 Ensembl VEP 或从 MyVariant.info 检索；没有工具输出时不要对影响进行分类。

---

CRISPR sgRNA 设计推理

• PAM 序列（SpCas9 的 NGG）必须位于非靶链上目标序列的 3' 端；向导 RNA 靶向 PAM 上游紧邻的 20 nt 序列。
• 靶向外显子时：选择在编码序列靠前位置切割的向导，以最大程度造成移码/破坏。
• 错配越少，脱靶风险越高；务必检查与向导序列有 0-3 个错配的基因组位点。

---

何时使用此技能

触发条件：

• 用户提供 VCF 文件（SNV/indel 或 SV）并询问其内容
• 关于变异等位基因频率（VAF）筛选的问题
• 突变类型分类查询（错义、无义、同义等）
• 结构变异解读请求（缺失、重复、CNV）
• 变异注释请求（ClinVar、gnomAD、CADD、dbSNP）
• 使用 ClinGen 剂量敏感性评估 CNV 致病性
• 队列比较问题
• 群体频率筛选（SNV 或 SV）
• 内含子/基因间变异筛选
• 基因剂量敏感性查询

示例问题：

• “VAF < 0.3 的变异中，有多少比例被注释为错义突变？”
• “筛除内含子/基因间变异后，还剩多少个非参考变异？”
• “这个影响 BRCA1 的缺失在临床上有何意义？”
• “chr17 上这个 500kb 重复区域与哪些剂量敏感基因重叠？”
• “有多少个变异带有临床意义注释？”
• “比较各样本之间的变异数量”

---

核心能力

能力	描述
VCF 解析	纯 Python + cyvcf2 解析器。支持 VCF 4.x、gzip 压缩文件、多样本、SNV/indel/SV
突变分类	将 SO 术语、SnpEff ANN、VEP CSQ、GATK Funcotator 映射为标准类型
VAF 提取	支持 AF、AD、AO/RO、NR/NV、INFO AF 格式
筛选	VAF、深度、质量、PASS、变异类型、突变类型、后果、染色体、SV 大小
统计	Ti/Tv 比率、每个样本的 VAF/深度统计、突变类型分布、SV 大小分布
注释	MyVariant.info（汇总 ClinVar、dbSNP、gnomAD、CADD、SIFT、PolyPhen）
SV/CNV 分析	gnomAD SV 群体频率、DGVa/dbVar 已知 SV、ClinGen 剂量敏感性
临床解读	使用单倍剂量不足/三倍剂量敏感性评分进行 ACMG/ClinGen CNV 致病性分类
DataFrame	转换为 pandas 格式以进行高级分析
报告	包含表格和统计数据的 Markdown 报告、SV 临床报告

---

工作流概览

阶段 1：解析 VCF → 提取 CHROM/POS/REF/ALT/QUAL/FILTER/INFO、每个样本的 GT/VAF/深度以及注释（ANN/CSQ/FUNCOTATION）。使用纯 Python 或 cyvcf2。

阶段 2：分类 → 变异类型（SNV/INS/DEL/MNV/SV）、突变类型（错义/无义/同义/移码/剪接等）、影响程度（HIGH/MODERATE/LOW/MODIFIER）。

阶段 3：筛选 → VAF 范围、深度、质量、PASS、变异/突变类型、排除后果、群体频率、染色体、SV 大小。

阶段 4：统计 → 类型/突变/影响程度/染色体分布、Ti/Tv 比值、每个样本的 VAF/深度、基因突变计数。

阶段 5：注释（可选）→ MyVariant.info（ClinVar/dbSNP/gnomAD/CADD）、Ensembl VEP 后果预测。

阶段 6：报告 → Markdown 表格、直接回答、导出 DataFrame。

阶段 7：SV/CNV 分析（如适用）→ gnomAD SV 频率、ClinGen 剂量敏感性、ACMG 致病性分类。

---

各阶段概要

阶段 1：VCF 解析

适用 pandas 的场景：

• 将 VCF 读取为结构化数据
• 快速探索性分析
• 需要操作列和行时

适用 python_implementation 工具的场景：

• 生产环境中的解析与注释提取
• 处理多样本 VCF
• 从 FORMAT 字段提取 VAF
• 大文件流式处理

关键函数：

vcf_data = parse_vcf("input.vcf")           # Pure Python (always works)
vcf_data = parse_vcf_cyvcf2("input.vcf")    # Fast C-based (if installed)
df = variants_to_dataframe(vcf_data.variants, sample="TUMOR")  # For pandas

阶段 2：变异分类

根据注释自动分类：

• SnpEff ANN 字段
• VEP CSQ 字段
• GATK Funcotator FUNCOTATION 字段
• 标准 INFO 键：EFFECT、EFF、TYPE

支持的突变类型：错义、无义、同义、移码、剪接位点、剪接区域、框内插入、框内缺失、内含子、基因间区、UTR_5、UTR_3、上游、下游、终止密码子丢失、起始密码子丢失

完整详情请参阅 references/mutation_classification_guide.md

阶段 3：筛选

常见筛选模式：

# Somatic-like variants
criteria = FilterCriteria(
    min_vaf=0.05, max_vaf=0.95,
    min_depth=20, pass_only=True,
    exclude_consequences=["intronic", "intergenic", "upstream", "downstream"]
)

# High-confidence germline
criteria = FilterCriteria(
    min_vaf=0.25, min_depth=30, pass_only=True,
    chromosomes=["1", "2", ..., "22", "X", "Y"]
)

# Rare pathogenic candidates
criteria = FilterCriteria(
    min_depth=20, pass_only=True,
    mutation_types=["missense", "nonsense", "frameshift"]
)

所有筛选选项请参阅 references/vcf_filtering.md

阶段 4–6：统计、注释、报告

使用 python_implementation 进行标准统计（Ti/Tv、类型分布、每个样本的 VAF/深度）；使用 pandas 进行自定义聚合。进行注释时，优先使用 MyVariant.info（批量处理：ClinVar + dbSNP + gnomAD + CADD）；每批限制为 50–100 个变异。报告包括类型/突变/影响程度/染色体分布、VAF 统计、临床意义和突变数最多的基因。

详细示例请参阅 references/annotation_guide.md

阶段 7：结构变异与 CNV 分析

当 VCF 包含 SV 检测结果时（SVTYPE=DEL/DUP/INV/BND）：

1. 识别受影响基因（根据 VCF 注释或坐标重叠）
2. 查询 ClinGen 剂量敏感性：

   clingen = ClinGen_dosage_by_gene(gene_symbol="BRCA1")
   # Returns: haploinsufficiency_score, triplosensitivity_score
   

3. 检查群体频率：

   gnomad_sv = gnomad_get_sv_by_gene(gene_symbol="BRCA1")
   # Returns: SVs with AF, AC, AN
   

4. 分类致病性：
   • 致病：缺失 + HI 评分 = 3，AF < 0.0001
   • 可能致病：缺失 + HI 评分 = 2，AF < 0.001
   • 临床意义未明：HI/TS 评分 = 0–1，AF 0.001–0.01
   • 良性：AF > 0.01

ClinGen 剂量评分解读：

• 3：有充分证据支持剂量致病性（影响程度高）
• 2：有一些证据（影响程度中等）
• 1：证据很少（影响程度低）
• 0：没有证据（影响程度极低）
• 40：不太可能存在剂量敏感性

完整 SV 工作流请参阅 references/sv_cnv_analysis.md

---

常见问题模式

模式 1：VAF + 突变类型占比

问题：“VAF < X 的变异中，有多少比例被注释为 Y 突变？”

result = answer_vaf_mutation_fraction(
    vcf_path="input.vcf",
    max_vaf=0.3,
    mutation_type="missense",
    sample="TUMOR"
)
# Returns: fraction, total_below_vaf, matching_mutation_type

模式 2：队列比较

问题：“队列之间的突变频率有何差异？”

result = answer_cohort_comparison(
    vcf_paths=["cohort1.vcf", "cohort2.vcf"],
    mutation_type="missense",
    cohort_names=["Treatment", "Control"]
)
# Returns: cohorts, frequency_difference

模式 3：筛选并计数

问题：“筛选掉 X 后，还剩多少个 Y？”

result = answer_non_reference_after_filter(
    vcf_path="input.vcf",
    exclude_intronic_intergenic=True
)
# Returns: total_input, non_reference, remaining

---

ToolUniverse 工具参考

SNV/Indel 注释

工具	使用场景	参数	响应
MyVariant_query_variants	批量注释	query（rsID/HGVS）	ClinVar、dbSNP、gnomAD、CADD
dbsnp_get_variant_by_rsid	群体频率	rsid	频率、临床意义
gnomad_get_variant	gnomAD 元数据	variant_id（CHR-POS-REF-ALT）	基本变异信息
EnsemblVEP_annotate_rsid	后果预测	variant_id（rsID）	转录本影响

欧洲变异档案（EVA）

EVA 是 EBI 的变异档案：它是一个存储原始、按研究提交的变异调用结果的仓库（小型测序项目队列，而非大型聚合群体数据库），并提供聚类变异（RS）登录服务。它与上述 MyVariant/dbSNP/gnomAD 互为补充，而非替代：

• gnomAD/MyVariant 为每个变异提供一个大型聚合群体等位基因频率（汇集数十万名个体）。
• EVA 的 cohortStats（已通过实时查询验证）按提交研究分别统计：每个 sourceEntries.<studyId>_<fileId>.cohortStats.ALL 块都包含其自身小型队列的 maf/altAlleleFreq（例如，一个只有几百个样本的外显子组测序项目）——应将其视为用于佐证 gnomAD 频率的队列特定证据，而非替代数据，也绝不要把队列级 MAF 求平均，误当作同一个群体的数据。
• EVA_list_studies 会显示这些队列背后实际命名的测序项目（疾病特定或群体特定研究）。当问题是“有哪些数据”而非“频率是多少”时，这很有用。
• EVA_get_clustered_variant_by_rs 属于 EVA 的另一项子服务（登录/标识符 API，而非变异数据 API）——它可将数字 RS 登录号（不带 rs 前缀，例如 429358）直接解析为组装版本/重叠群/位置记录；无需先选择基因或区域，适合快速核对坐标。

工具	使用场景	关键参数	响应
EVA_get_variants_by_gene	查询 EVA 提交的某个基因的所有变异	gene（HGNC 符号）	response[].result[]，每项含有 chromosome（RefSeq 重叠群登录号，如 CM000679.2，而非普通染色体编号）、start、reference/alternate、hgvs.genomic[]、按研究划分的 sourceEntries
EVA_get_variants_by_region	扫描坐标范围（热点、调控区域）	region（chr:start-end，此处使用普通染色体编号，与响应中的重叠群登录号不同）	与上方结构相同，但范围限定于该区域
EVA_list_studies	查询变异前，了解有哪些测序项目/队列	species（默认值为 hsapiens_grch38）	{studyId, studyName, filesCount} 组成的 result[]
EVA_get_clustered_variant_by_rs	将 RS 登录号快速解析为坐标（登录 API，与上方的变异数据 API 不同）	accession（数字 RS id，不带 rs 前缀）	[{accession, version, data: {assemblyAccession, taxonomyAccession, contig, start, type}}]

分页说明：响应中的 numTotalResults 可能远大于当前返回页的结果数（例如，BRCA1 在 limit 为 3-20 时，共有 6,684 个 EVA 提交的变异）——切勿将当前返回的少量结果视为完整集合；在概述“变异数量”前，请先检查 numTotalResults。

Ensembl 表型关联

汇总多个上游数据源（NHGRI-EBI GWAS catalog、ClinVar、OMIM、Cancer Gene Census、Orphanet、dbGaP）的表型／疾病关联，并通过一个 Ensembl REST 端点提供查询，可按基因、区域、表型术语或变异进行查询。实时可靠性已反复验证，并非假定：目前只有 EnsemblPheno_get_by_variant 可用——2/2 次实时调用均成功并返回了真实数据。EnsemblPheno_get_by_gene、EnsemblPheno_get_by_region 和 EnsemblPheno_get_by_term 每次尝试都失败（get_by_gene 会挂起／超时；另外两个会快速返回“HTTP error: unknown”）。问题源自 Ensembl 自身的 REST API，而非 ToolUniverse：直接对 rest.ensembl.org/phenotype/gene/homo_sapiens/BRCA1 执行 curl，也独立收到了来自 EBI 的实时 HTTP 500 错误。请放心使用 EnsemblPheno_get_by_variant；目前其他三个工具在上游已损坏，依赖它们之前请通过新鲜的 tu run 调用重新检查——如果它们再次失败，请勿捏造基因／区域／术语级别的表型结果。

工具	适用场景	关键参数	响应	实时状态
EnsemblPheno_get_by_variant	与一个 rsID 关联的表型／性状，汇总多个来源	variant_id（rsID）	phenotypes[]，包含 trait、source、risk_allele、pvalue、beta_coefficient、study	可用（已验证）
EnsemblPheno_get_by_gene	与一个基因关联的表型，多来源汇总	gene（符号）	phenotypes[]，包含 description、source、ontology_accessions	已损坏（上游 HTTP 500，已验证）
EnsemblPheno_get_by_region	基因组区间的表型概况	region（chr:start-end）	phenotypes[]，限定于该区域	上游已损坏（已验证）
EnsemblPheno_get_by_term	反向查询：表型名称／本体词条编号 -> 变异／基因	term 或 accession（EFO/HP/MONDO）	associations[]，包含 variant、gene、risk_allele、p_value、odds_ratio	上游已损坏（已验证）

真实示例（rs429358，与上文“变异命名转换”中使用的 APOE 变异相同——已交叉核对一致性）： EnsemblPheno_get_by_variant(variant_id="rs429358") 返回了 1140 条真实的表型关联，其中包括脂质／代谢物 GWAS 结果（例如“1,2-dihydroxy-3-keto-5-methylthiopentene dioxygenase levels”，PMID 39528825），数据来源为 NHGRI-EBI GWAS catalog；另一个 rsID rs7903146（TCF7L2/糖尿病）返回了 331 条关联，其中包括 MAGIC 联盟的葡萄糖性状结果。每条 phenotypes[] 中的 source 字段会标明其上游数据库来源——请据此评估证据等级（ClinVar 来源的条目与单项研究的 GWAS 结果权重不同）。

当 EnsemblPheno_get_by_gene 不可用时，如需进行基因级别的表型汇总，请改用 Monarch 或 Gene2Phenotype（参见 tooluniverse-gene-disease-association 的 Phase 4/Phase 5）。

结构变异注释

工具	适用场景	参数	响应
gnomad_get_sv_by_gene	SV 群体频率	gene_symbol	包含 AF、AC、AN 的 SV
gnomad_get_sv_by_region	区域 SV 搜索	chrom、start、end	区域内的 SV
ClinGen_dosage_by_gene	剂量敏感性	gene_symbol	HI/TS 分数、疾病
ClinGen_dosage_region_search	区域内剂量敏感基因	chromosome、start、end	所有带有 HI/TS 分数的基因
ensembl_get_structural_variants	来自 DGVa/dbVar 的已知 SV	chrom、start、end、species	临床意义

详细的工具用法示例请参见 references/annotation_guide.md

变异表示法转换（SPDI / HGVS / VCF / rsID）

不同的数据库和工具要求使用不同的变异表示法——NCBI 的规范内部格式是 SPDI（SeqID:Position:Deleted:Inserted），ClinVar 和大多数临床报告使用 HGVS（NC_000019.10:g.44908684T>C），VCF 文件使用 CHROM POS REF ALT，而 dbSNP 使用 rsID 标识变异。在不同格式之间转换（或规范化为一种规范形式以便去重），是使用此技能进行其他注释查询前常见的预处理步骤。

工具	适用场景	关键参数	响应
NCBIVariation_vcf_to_spdi	有原始 VCF 字段，需要 SPDI 以进入此流程	chrom（RefSeq 登录号，例如 NC_000019.10）、pos、ref、alt	规范化后的 spdis[]
NCBIVariation_spdi_to_hgvs	有 SPDI，需要 HGVS 用于临床报告或 ClinVar 风格的查询	spdi	hgvs 字符串
NCBIVariation_hgvs_to_spdi	有 HGVS 描述（g./c./r.），需要 SPDI 或想验证 HGVS 语法	hgvs	spdis[]（一个 HGVS 可能对应多个基因组组装或转录本上的多个 SPDI）
NCBIVariation_spdi_canonical	需要某个插入缺失变异唯一的规范化／右移 SPDI（用于去重和一致存储）	spdi	一个规范 spdi
NCBIVariation_spdi_equivalents	需要同一变异在 GRCh37／GRCh38／RefSeqGene／转录本上的坐标（坐标提升）	spdi	不同坐标系统上的 equivalents[]
NCBIVariation_rsid_lookup	有 rsID，需要通过一次调用获取基因组坐标、基因、MANE Select 转录本和 ClinVar 临床意义	rsid	grch38_placements[]、genes[]、clinical_significance[]、mane_select_ids[]
NCBIVariation_spdi_to_rsids	有基因组 SPDI，需要获取共定位的 dbSNP rsID（与 rsid_lookup 方向相反）	spdi	rsids[]
NCBIVariation_alfa_frequencies_by_rsid	需要按祖源划分的等位基因频率（欧洲、非洲、东亚、南亚、拉丁美洲等），而不是单一的全球 MAF	rsid	positions[].studies[].populations[]，每项均包含按祖源划分的 allele_counts 和 allele_frequencies

与上文的 gnomAD／MyVariant 频率数据不同：gnomad_get_variant 和 MyVariant_query_variants 返回的是群体层面的 gnomAD 频率；NCBIVariation_alfa_frequencies_by_rsid 则是 NCBI 自有的 ALFA 汇总数据，报告的祖源分类不同（同一变异的数值不一定与 gnomAD 一致）。不要将它们视为「群体频率」声明的可互换数据源；报告频率时应注明来源是 gnomAD 还是 ALFA。

选择指南——从已有数据出发，找到所需结果： VCF 字段 → SPDI（vcf_to_spdi）→ HGVS（spdi_to_hgvs）或 rsID（spdi_to_rsids）；rsID → 坐标／基因／ClinVar（rsid_lookup）或按祖源划分的频率（alfa_frequencies_by_rsid）；HGVS → SPDI（hgvs_to_spdi），以便重新进入 SPDI 流程（获取等效表示、规范化、查询 rsID）。

完整参数表、真实调用及响应示例，以及 rs429358（APOE）的完整示例流程，请参阅 references/variant_notation_conversion.md。

HGVS 验证与精选变异数据库（Mutalyzer、LOVD、ClinGen Allele Registry、GeneBe、APPRIS）

另外五组小型工具系列均已通过实测，可完善临床变异工作流：语法验证、基因特异的精选变异查询、跨数据库的规范等位基因 ID、自动 ACMG 分类，以及主转录本选择。

Mutalyzer——HGVS 语法验证／规范化（与上述 NCBI Variation 表示法转换工具互为补充，并不重复）： NCBIVariation_hgvs_to_spdi／spdi_to_hgvs 用于在不同表示法之间转换；Mutalyzer 则针对实际参考序列验证并纠正单个 HGVS 字符串的语法，将其转换为预测的蛋白质／RNA 效应，还可以将蛋白质变异反向转换为可能的编码区描述。

工具	使用时机	关键参数	响应
Mutalyzer_normalize_variant	在其他地方使用 HGVS 描述前，验证或更正该描述	variant（例如 NM_000546.6:c.215C>G）	corrected_description、预测的蛋白质/RNA 效应、基因符号
Mutalyzer_parse_hgvs	将 HGVS 字符串拆分为结构化部分（参考序列、坐标系统、变异类型）	variant	结构化的 model（参考序列、坐标系统、变异类型/位置）
Mutalyzer_back_translate	已知蛋白质变化（p.），需要可能的 DNA 层面描述	variant（例如 NP_000537.3:p.Pro72Arg）	dna_descriptions[]（可能多于 1 个 — 由于密码子简并性）

LOVD — 基因专属的人工整理变异数据库：与 MyVariant/EVA 广泛覆盖多个基因不同，LOVD 托管由基因专属整理者维护的 各基因人工整理变异集合（对 TP53/BRCA1/BRCA2 等研究充分的基因 覆盖较好，对其他基因的数据较少或没有）。

工具	使用时机	关键参数	响应
LOVD_get_gene	确认某个基因是否有 LOVD 条目，并获取其转录本/基因组版本元数据	gene_symbol	HGNC/Entrez ID、RefSeq 转录本、整理日期
LOVD_get_variants	列出某个基因的所有人工整理变异	gene_symbol	变异的 []，包含 DNA/RNA/蛋白质 HGVS、hg19 位置、Times_reported
LOVD_search_variants	通过 LOVD DBID 或 DNA 表示法查找某个特定变异	gene_symbol +（variant_dbid 或 dna_notation）	与上述相同的单个变异格式

ClinGen Allele Registry — 跨数据库的规范等位基因 ID：两组 工具均覆盖同一个底层注册库，功能有所重叠（实时验证 — CA ID 相同， 外部记录格式相同），但工具类别不同（ClinGenARTool 与 ClinGenAlleleTool）。优先使用 ClinGenAR_* — 它具备与 ClinGenAllele_* 相同的正向查询/详情调用，外加 ClinGenAllele_* 不具备的按外部 ID 反向查询功能。只有在已有 上下文使用 ClinGenAllele_* 时，才使用它。

工具	使用时机	关键参数	响应
ClinGenAR_lookup_allele	已有 HGVS 字符串，需要规范的 CA ID	hgvs	allele_id（CA...）、社区标准标题
ClinGenAR_get_external_records	已有 CA ID，需要 ClinVar/dbSNP/COSMIC/gnomAD/ExAC 的交叉引用	allele_id	按来源数据库分组的 external_records
ClinGenAR_lookup_by_external_id	已有 rsID 或 ClinVar VariationID，需要 CA ID（与上述查询方向相反）	dbsnp_rs 或 clinvar_variation_id（不可同时使用）	alleles[]（每个外部 ID 可能对应多个 CA ID）

GeneBe — 自动化 ACMG/AMP 分类：提供原始基因组坐标（非 HGVS）， 即可在一次调用中返回完整的 ACMG 判定及所采用的具体标准 （PS3/PM1/PM2 等）、ClinVar 交叉核查，以及 AlphaMissense/gnomAD 背景信息 — 可用于专家人工 ACMG 审查前快速自动初筛，但不能替代 专家整理。

工具	使用时机	关键参数	响应
GeneBe_classify_variant	按染色体/位置/参考/替代碱基对一个变异进行分类	chr、pos、ref、alt、genome（默认 hg38）	acmg_classification、acmg_criteria、clinvar_classification、alphamissense_score、gnomad_exomes_af
GeneBe_classify_variants_batch	在一次请求中对最多 1000 个变异进行分类	variants[]（每个元素为 {chr,pos,ref,alt}）、genome	每个变异的 ACMG 结果，字段与上面相同

APPRIS — 主要转录本/异构体选择：当一个基因有多个转录本， 且某个变异的编码后果取决于使用哪个转录本时，APPRIS 会根据蛋白质 结构/保守性证据指出哪个异构体在功能上属于“主要”异构体， 并区分其他异构体 — 可用于在解读编码型 HGVS 描述前选择正确的转录本。

工具	使用场景	关键参数	响应
APPRIS_get_isoforms	列出某个基因的所有转录本及其 PRINCIPAL/ALTERNATIVE 标签	gene_id（Ensembl 基因 ID）、species	每个转录本对应一个 []，包含 type（principal_isoform/alternative）和标签
APPRIS_get_principal_isoform	直接获取唯一的 principal 转录本	gene_id、species	transcript_id、ccds_id、length_na
APPRIS_get_functional_annotations	查看支持 principal isoform 判定的各方法详细证据（firestar/spade/matador3d/corsair 等）	gene_id、species、可选的 methods、transcript_id	各方法的分数/注释

真实串联示例（TP53，使用相同的变异坐标 NC_000017.11:g.7674220C>T / NM_000546.6:c.743G>A 串联各工具，完成连贯的操作示例，且已实时验证）： Mutalyzer_normalize_variant("NM_000546.6:c.215C>G") → 更正后的描述 + 预测的蛋白质影响。LOVD_get_variants("TP53") → 真实的已整理条目（例如，LOVD DBID TP53_010464， NM_000546.5:c.*2609C>A）。ClinGenAR_lookup_by_external_id(clinvar_variation_id="376694") → 真实的 allele_id="CA16040589"、标题 NM_000546.6(TP53):c.706T>A (p.Tyr236Asn)，并带有 COSMIC 交叉引用。GeneBe_classify_variant(chr="17", pos=7674220, ref="C", alt="T", genome="hg38") → 实时查询到的真实结果： acmg_classification="Pathogenic"、acmg_score=22、 dbsnp="rs11540652"、alphamissense_score=0.996、 gnomad_exomes_af=6.16e-06——与 ClinVar 对相同坐标给出的“Pathogenic”判定完全吻合。APPRIS_get_principal_isoform("ENSG00000141510") → TP53 的真实 principal 转录本，ccds_id="CCDS11118.1"。

连锁不平衡（Ensembl）

LD 衡量不同位点的等位基因是否以非随机方式共同遗传——这对于 GWAS 解读（命中变异本身是否具有因果作用，还是仅与真正的因果变异相关？）和精细定位（某个位点中的哪些变异应一起检测）至关重要。

工具	使用场景	关键参数	响应
EnsemblLD_get_ld_variants	查找与某个查询变异存在 LD 的所有变异，限定在一个人群中	variant_id（rsID）、population（1000GENOMES:phase_3:<POP>）、可选的 r2_threshold/d_prime_threshold/limit	ld_variants[] 按 r2 降序排列；truncated/total_ld_count——如果结果被截断，请提高 limit
EnsemblLD_get_ld_pairwise	一次性检查两个特定变异在所有 1000 Genomes 人群中是否相关	variant1、variant2（rsID）	ld_by_population[]，每项包含 r2/d_prime——LD 会因祖源而显著不同（实时验证：rs6792369/rs1042779 在某些人群中的 r2=1.0，在其他人群中为 0.84）
EnsemblLD_get_ld_region	计算一个窗口内所有变异之间的完整成对 LD 矩阵	region（chr:start..end，GRCh38，<=1Mb）、population	ld_pairs[]——这是统计精细定位（SuSiE、FINEMAP）和 LD 感知型聚类所需的 LD 矩阵输入

EnsemblLD_get_ld_variants(variant_id="rs429358", population="1000GENOMES:phase_3:CEU", r2_threshold=0.5)
#  -> real: LD partners for the APOE variant (rs429358, also used in
#     Variant Notation Conversion and EVA above) in the CEU population

实时验证的响应时间：get_ld_variants/get_ld_pairwise 为 2-8s， 小范围 get_ld_region 窗口的响应时间相近——无需特殊的超时处理。

坐标转换与蛋白质/cDNA 到基因组的映射（Ensembl）

工具	使用时机	关键参数	响应
EnsemblMap_convert_coordinates	在不同基因组组装版本之间迁移变异/BED 区间（例如 GRCh37 -> GRCh38）	species、source_assembly、chromosome、start、end、target_assembly	mappings[]，包含 original/mapped 区块；每个区块都带有 seq_region_name/start/end/strand
EnsemblMap_translate_coordinates	将蛋白质氨基酸范围或转录本 cDNA 位置范围映射到基因组坐标	ensembl_id（蛋白质使用 ENSP*，cDNA 使用 ENST*）、start、end（从 1 开始，在对应坐标空间中）	基因组坐标的 mappings[]；所跨越的每个外显子对应一个条目

EnsemblMap_convert_coordinates(species="human", source_assembly="GRCh37",
    chromosome="7", start=140453136, end=140453136, target_assembly="GRCh38")
#  -> real: BRAF V600E's GRCh37 position 7:140453136 maps to GRCh38 7:140753336

EnsemblMap_translate_coordinates(ensembl_id="ENSP00000269305", start=100, end=200)
#  -> real: TP53 protein residues 100-200 map to 3 exons at
#     chr17:7674931-7676071 (GRCh38) — use this before designing an assay
#     or interpreting a reported amino-acid range against the genome

实时验证的耗时：convert_coordinates 很快（7-24 秒）； translate_coordinates 最长可能需要约 60 秒——请设置充足的超时时间， 不要因响应较慢就认为请求失败。

Ensembl ID 实用工具（存档、重叠、交叉引用、组装/物种元数据）

通用的 Ensembl ID 管理工具；当变异分析工作流需要解析、验证或了解基因/转录本 ID 的上下文，而不是注释特定变异时，这些工具很有用。

工具	适用场景	关键参数	响应
EnsemblArchive_get_id_history	检查稳定 ID（ENSG*/ENST*/ENSP*）是否仍为当前 ID，以及若已过时则查询其替代 ID	ensembl_id	is_current、latest_version、current_release、possible_replacement[]
EnsemblArchive_batch_lookup	一次检查最多 50 个 ID	ensembl_ids（逗号分隔）	entries[]，每个 ID 对应一项
Ensembl_get_region_features	获取基因组窗口内的所有重叠内容——基因、转录本、调控元件、受约束元件	region、species、feature_types（逗号分隔：gene,transcript,regulatory,constrained,variation,repeat,motif）	features[]、type_summary
Ensembl_get_gene_overlapping_features	获取与某个基因共定位的所有内容（重叠基因、转录本异构体、该位点的调控特征）	gene_id、feature_types	features[]、type_summary
Ensembl_get_cross_references	将 Ensembl ID 映射到 HGNC/EntrezGene/UniProt/OMIM/Reactome/GeneCards 及其他外部数据库记录	ensembl_id，可选的 external_db 筛选条件	xrefs[]，包含 dbname/primary_id/display_id，以及 database_summary
Ensembl_lookup_gene_by_symbol	将基因符号解析为 Ensembl 基因/转录本/蛋白质 ID	symbol、species	ensembl_ids[]，包含 id/type
Ensembl_get_assembly_info	查询物种的基因组组装元数据（组装名称/登录号、核型、坐标系统版本）	species（例如 homo_sapiens）	assembly_name、karyotype[]、coordinate_system_versions[]
Ensembl_get_species_info	按名称/常用名/分类群 ID 搜索，了解 Ensembl 提供了哪些基因组/组装	search（可选——省略则列出全部 348+ 个物种）	species[]，包含 taxon_id、assembly、division

关于可靠性的说明（经多次尝试实时验证，而非仅凭单次运行推断）：

• EnsemblArchive_*、Ensembl_get_cross_references 和 Ensembl_get_assembly_info/Ensembl_get_species_info 都是真正可用的端点，但速度较慢，偶尔会暂时不稳定——首次尝试有时会快速返回 HTTP 500/超时，但将超时时间延长至 60-90s 后重试即可成功。将这些端点报告为故障前，请先延长超时时间重试一次。
• Ensembl_get_region_features 在 feature_types 包含多个值（例如 "gene,regulatory"）时会稳定失败（HTTP 500），至少 TP53 区域如此；单独查询每种特征类型（分别使用 "gene" 或 "regulatory"）则每次都能成功。解决方法：每次调用只查询一种特征类型，并在客户端合并结果，而不要依赖包含多种类型的合并查询。
• Ensembl_lookup_gene_by_symbol 已确认是上游故障，并非 ToolUniverse 的问题：它始终返回 Ensembl 网站的 HTML 错误页面（而非 JSON），多次尝试并使用多个符号（TP53、BRCA1）均重现了该问题；此外，在 ToolUniverse 之外直接对相同的 xrefs/symbol/<species>/<symbol> 端点执行 curl 也独立确认了这一点，该请求同样卡住/失败。如果此工具失败，请勿捏造符号到 ID 的映射——如果已有 ID，可反向使用 Ensembl_get_cross_references（Ensembl ID -> 外部名称）；否则，请通过其他已有文档说明的工具解析符号（例如 NCBIVariation_rsid_lookup 的 genes[] 字段，或 ToolUniverse 其他位置的 MyGene/UniProt 查询）。

Ensembl_get_cross_references(ensembl_id="ENSG00000141510")
#  -> real: TP53's HGNC/EntrezGene/UniProt/OMIM/Reactome cross-references
#     (157 xrefs across 11 databases, verified live)

Ensembl_get_region_features(region="7:140424943-140524564", feature_types="gene")
Ensembl_get_region_features(region="7:140424943-140524564", feature_types="regulatory")
#  -> query separately and merge; the combined "gene,regulatory" form
#     500s on this region (see reliability note above)

---

常见使用模式

# Quick summary
report = variant_analysis_pipeline("input.vcf", output_file="report.md")

# Filtered analysis
report = variant_analysis_pipeline("input.vcf",
    filters=FilterCriteria(min_vaf=0.1, min_depth=20, pass_only=True))

# Annotated report (top 50 variants with ClinVar/gnomAD/CADD)
report = variant_analysis_pipeline("input.vcf", annotate=True, max_annotate=50)

pandas 与 python_implementation 的对比：使用 python_implementation 进行解析、分类和注释，然后转换为 DataFrame 以执行自定义聚合：

vcf_data = parse_vcf("input.vcf")
passing, _ = filter_variants(vcf_data.variants, criteria)
df = variants_to_dataframe(passing, sample="TUMOR")

---

限制

• 突变分类需要 VCF 注释：如果 VCF 的 INFO 中没有 ANN/CSQ/FUNCOTATION，应用 ToolUniverse 注释之前，突变类型将显示为“unknown”
• 多等位基因变异：解析器在类型分类时会采用第一个 ALT 等位基因
• ToolUniverse 注释速率：基于 API，默认每批最多约 100 个变异，以遵守速率限制
• gnomAD 工具：仅返回基本元数据（不包含完整等位基因频率）；如需获取 gnomAD AF，请使用 MyVariant.info
• 大型 VCF：纯 Python 解析器会逐行读取；对于超过 100K 个变异的文件，建议使用 cyvcf2

---

参考文档

• references/vcf_filtering.md：完整的过滤选项和示例
• references/mutation_classification_guide.md：突变类型分类规则详解
• references/annotation_guide.md：包含示例的 ToolUniverse 注释工作流
• references/sv_cnv_analysis.md：完整的 SV/CNV 解读工作流
• references/variant_notation_conversion.md：SPDI/HGVS/VCF/rsID 转换工具，包含实际调用和响应示例

---

其他资源

• 主要脚本（请优先使用这些脚本——参见本文件顶部）：
  • scripts/gatk_haplotypecaller_pipeline.py — BWA + HaplotypeCaller + SNP/indel 计数器
  • scripts/coding_variant_filter.py — 按样本统计外显子变异，并排除内含子/UTR
  • scripts/variant_fraction.py — VAF + 编码区分母比例计算器
• 通用脚本：scripts/parse_vcf.py、scripts/filter_variants.py、scripts/annotate_variants.py
• 快速入门方案和 MCP 示例：QUICK_START.md

---

分析惯例

多行 Excel 表头（三人组 / CHIP 变异表）

临床变异导出电子表格通常有两行表头（类别行，例如 Variant Info / Flags / Father (185-PF) / RefSeq Genes 110, NCBI，位于子标签行，例如 Chr:Pos / Variant Allele Freq / Sequence Ontology (Combined) 上方）。使用 pd.read_excel(path, header=[0,1]) 解析，并通过元组访问列，例如 df[('Father (185-PF)', 'Variant Allele Freq')]。单行表头会将子列标记为 Unnamed:_N，并导致 VAF / Sequence Ontology 数据被静默遗漏。

“注释为 X 的变异比例”——分母是编码区变异

如果问题询问在某个 VAF / 过滤阈值下的变异中，有多少比例被注释为某个 Sequence Ontology 术语（例如 synonymous_variant），分母应为编码区子集，而非“所有变异”（其中内含子记录占主导）。

CODING = {
    "synonymous_variant", "missense_variant", "splice_region_variant",
    "stop_gained", "stop_lost", "start_lost",
    "frameshift_variant", "inframe_insertion", "inframe_deletion",
}
NON_CODING = {  # explicitly excluded
    "intron_variant", "3_prime_UTR_variant", "5_prime_UTR_variant",
    "upstream_gene_variant", "downstream_gene_variant", "intergenic_variant",
}

注意，splice_region_variant 属于编码区（会影响剪接边界处的编码序列）——应将其纳入。

❌ 错误：count(VAF<0.3 AND synonymous) / count(VAF<0.3) — 分母混入了内含子/UTR ✅ 正确：count(VAF<0.3 AND synonymous) / count(VAF<0.3 AND CODING) — 分母仅限编码区

请通过 df[df['Sequence Ontology (Combined)'].isin(CODING)] 进行筛选，不要通过排除 NON_CODING 标签来筛选——后者会导致计数过多，因为某些标签（例如 5_prime_UTR_premature_start_codon_gain_variant）不属于任一集合。

合理性检查：在人类外显子组中，同义变异通常约占编码区变异的一半。如果你的“同义变异比例”远低于约 0.4，分母可能仍包含内含子/UTR——请限定为 CODING 后重新计算。

内置工具：tu run coding_variant_fraction '{"file":"variants.xlsx","vaf_threshold":0.3,"annotation":"synonymous_variant","header_rows":2}' — 自动处理两行表头和 CODING 允许列表。