尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

变异表示与归一化:用 left-alignment 与 parsimony 规则消除 VCF 变异记录的歧义(scientific-agent-skills 实战指南)

变异表示与归一化:用 left-alignment 与 parsimony 规则消除 VCF 变异记录的歧义(scientific-agent-skills 实战指南) 变异表示与归一化用 left-alignment 与 parsimony 规则消除 VCF 变异记录的歧义scientific-agent-skills 实战指南【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills同一处基因组改变可以写成无数条字段各不相同的变异记录两条记录可能没有一个字段值相同却描述同一个变异两条POS完全相同的记录也可能是两个不同的变异。任何在归一化之前进行的比较、合并、去重或注释查询都会静默丢失真实的匹配——不会报错交集只是比应有的更小。本指南基于scientific-agent-skills仓库的genomic-coordinates技能完整讲解 VCF 变异的归一化parsimony 修剪 left-alignment 左对齐原理、normalize_variant.py脚本的全部用法以及多等位基因拆分、参考基因组校验、HGVS 反向规则和结构变异处理读完即可在任何跨数据集的变异比较前正确执行归一化。为什么一个变异有无数种拼写考虑这样一条参考序列position 1 2 3 4 5 6 7 8 9 10 base G G C A C A C A C T从CACACAC重复区段中删除一个AC无论删的是哪一对相邻的AC得到的序列都是GGCACACT。下面四条记录全部是同一个变异POS7 REFCAC ALTC POS5 REFCAC ALTC POS3 REFCAC ALTC POS2 REFGCA ALTG任何测序流程都可能输出其中任意一条。而重复区段正是 indel 高度聚集的地方也是最需要归一化、歧义最严重的区域。仓库中的测试夹具 tests/genomic-coordinates/fixtures/tiny.fa 特意构造了一条长度为 10 的参考序列GGCACACACT第 3–9 位是CACACAC重复目的就是在离线环境里让左对齐的用例真正有趣见 tests/genomic-coordinates/test_scripts.py。冗余侧翼碱基构成第二个歧义维度。POS3 REFCA ALTCT与POS4 REFA ALTT是同一个 SNV前者只是带上了一个没有发生变化的碱基。VCF 的POS指向的是REF的第一个碱基区间是POS到POS len(REF) - 1对 indel 而言POS是锚定碱基——事件之前那个本身未发生改变的碱基详见 skills/genomic-coordinates/references/format-conventions.md。归一化规则parsimonious left-aligned一个变异被认为是归一化的当它同时满足parsimonious最简碱基数尽可能少但至少保留一个left-aligned左对齐在不改变其描述的序列的前提下向 contig 起始方向尽可能平移。这是 Tan、Abecasis 与 Kang 于 2015 年发表的Unified representation of genetic variantsBioinformatics 31(13):2202–2204给出的定义也是bcftools norm与vt normalize所实现的规则。skills/genomic-coordinates/scripts/normalize_variant.py 的模块注释与源码正是对这一过程的忠实实现。完整算法分两步当所有等位基因都以相同碱基结尾时若某个等位基因已只剩一个碱基则向左取参考序列的一个碱基扩展每个等位基因并将POS减一然后删除每个等位基因的最后一个碱基。当每个等位基因都至少有两个碱基、且都以相同碱基开头时删除每个等位基因的第一个碱基并将POS加一。第 1 步让变异在重复区段中向左行走第 2 步剥离冗余的填充碱基。两步都会终止。注意POS的语义第 1 步中POS递减向左移动锚点第 2 步中POS递增修剪后锚点右移到更短、等价的记录上。源码实现normalize_variant.py 核心逻辑normalize_variant.py对单条双等位基因记录的归一化实现在normalize()函数中关键流程对应 skills/genomic-coordinates/scripts/normalize_variant.py前置校验POS 1或REF为空直接抛出VariantErrorREF ALT断言没有发生改变同样拒绝。REF 对照用reference.fetch(contig, pos - 1, pos - 1 len(ref))从 FASTA 取实际序列与REF比对不一致则返回ref_check MISMATCH并停止归一化。右端修剪与左向扩展循环判断ref[-1] alt[-1]若某个等位基因只剩 1 个碱基则向左取一个参考碱基前置并pos - 1随后统一丢弃尾部碱基--window默认 1000 bp限制向左扩展的距离触底时在detail中提示重复区段可能延伸更远。左端修剪当len(ref) 1且len(alt) 1且首碱基相同时删除首碱基并pos 1为 indel 保留一个锚定碱基。运行方式需 Python 3.11仅标准库无第三方依赖、无网络访问python3 normalize_variant.py --fasta ref.fa chr1 7 CAC C # chr1:7:CAC:C - chr1:2:GCA:G pos_shift 5pos_shift源码字段shifted为正表示左对齐让锚点穿过了重复区段向左移动为负表示修剪冗余侧翼碱基后锚点向右落到了更短、等价的记录上。例如3 CACT实际是 SNV4 ATpos_shift -1。变异类型由classify()判定snv长度均为 1、mnv长度相等且大于 1、insertionlen(ref)1且alt.startswith(ref)、deletionlen(alt)1且ref.startswith(alt)、否则为complex。脚本支持的完整命令行参数参数说明--fasta必填参考 FASTA 路径存在.fai索引时使用索引做随机读取否则整条 contig 读入内存--input输入 VCF5 列以上按CHROM POS ID REF ALT读取或裸四列 TSVCHROM POS REF ALT跳过#、track、browser、开头行--compare归一化两条或多条CONTIG:POS:REF:ALT记录并报告它们是否为同一变异--split将逗号分隔的多等位基因 ALT 拆成每条记录后再归一化--window左对齐窗口默认 1000 bp用于限制向左扩展距离--format输出格式tsv默认或json-o, --output输出文件路径缺省写 stdout退出码0表示所有记录都通过了参考序列校验1表示至少一条记录的REF不匹配或记录无效2表示用法或参考序列错误例如 FASTA 不存在。等价性检查--compare归一化之后比较两个变异只需比较归一化后的四个字段。把两条记录都归一化再比较即可python3 normalize_variant.py --fasta ref.fa \ --compare chr1:7:CAC:C chr1:3:CAC:C chr1:2:GCA:G # verdict: identical -- all 3 records normalise to chr1:2:GCA:G判定逻辑见源码main()所有记录的ref_check都是ok时才可比归一化后的 key 集合大小为 1 输出identical否则输出distinct。verdict 输出到 stderrstdout 上的逐记录表格保持可解析。归一化是幂等的——对已归一化的记录再次归一化结果不变这一性质由 tests/genomic-coordinates/test_scripts.py 中的test_normalisation_is_idempotent覆盖验证。归一化必须使用正确的参考序列左对齐需要读取参考碱基。如果传入错误的组装版本脚本会自信地给出一个错误的答案。因此每条记录的REF会先与 FASTA 核对不匹配就停止该记录ref_check MISMATCH REF says A but the reference has C at chr1:3REF不匹配是成本最低的组装错配探测器。如果失败记录超过寥寥数条说明变异与 FASTA 属于不同构建——此时应当运行check_contigs.py而不是去调整任何坐标python3 check_contigs.py --identify unknown.fa.fai python3 check_contigs.py variants.vcf annotation.gtf --genome GRCh38.fa.faicheck_contigs.py读取.fai、.chrom.sizes、VCF 头部、SAM 头部、FASTA、BED、GTF/GFF按主染色体长度识别组装并报告两文件 join 可能出错的每一种原因命名不一致、长度冲突、坐标越过 contig 末端、只出现在单侧文件中的 contig任何不兼容时退出码为 1。组装签名表与 GRCh37 vs hg19二者仅在线粒体上不同16,569 bp 的 rCRS 与 16,571 bp等细节见 skills/genomic-coordinates/references/reference-builds.md。测试 tests/genomic-coordinates/test_scripts.py 用GRCh37.fa.fai、hg19.chrom.sizes与grch38.vcf等夹具验证了按线粒体长度区分 GRCh37/hg19、GRCh38 VCF 对上 hg19 基因组报different assemblies、BED 坐标越过 contig 末端被捕获等行为。多等位基因记录必须先拆分再归一化ALTG,GG是共享一行的两个变异。它们必须在归一化之前拆分因为让它们能同处一行的共享REF对其中任何一个都不是最简REFpython3 normalize_variant.py --fasta ref.fa --split --input cohort.vcf先归一化再拆分、或者把多等位基因记录当作整体归一化得到的每条记录都是单独错误的。bcftools norm -m -any -f ref.fa按正确顺序同时完成这两件事。注意拆分会重写基因型genotype与INFO字段——NumberA的逐等位基因INFO条目随拆分一并切分其余条目则复制到两条记录上。--split同样作用于命令行直接给出的变异--split chr1 7 CAC C,CACAC会先展开成两条记录再分别归一化。另一个方向HGVS 向右移动VCF 左对齐HGVS 恰恰相反在存在歧义的情况下参考序列最 3 端的位置被任意指定为发生改变的位置。两个标准刻意方向相反且这种差异是真实的——同一个删除在 VCF 与临床报告HGVS中有不同的坐标。更麻烦的是HGVS 的3是相对于被描述的参考序列而言的描述向何处移动在正义链基因上在负义链基因上VCFPOScontig 起点基因组最左基因组最左HGVSg.contig 终点基因组最右基因组最右HGVSc./n./p.转录本 3 端基因组最右基因组最左因此对于负义链基因HGVSc.描述可能与左对齐的 VCF 记录恰好一致对正义链基因则系统性地不一致。永远不要通过调整坐标在两者之间转换要通过掌握转录本模型的工具往返bcftools csq、VEP、Mutalyzer、Python 的hgvs包。转录本、CDS、蛋白坐标之间的映射规则c.1是起始 ATG 的 A、无c.0、5 UTR 为负、3 UTR 用*、GFF phase 是需移除的碱基数而非start % 3见 skills/genomic-coordinates/references/transcript-coordinates.md。符号等位基因与结构变异DEL、DUP、INV、CNV、INS以及断点BND记录不携带字面序列REF只是POS处单个锚定碱基范围存在于INFO/END与INFO/SVLEN中。它们无法归一化normalize_variant.py对alt.startswith((, *, .))的记录直接以ref_check skipped原样放行并在detail中注明symbolic, missing, or spanning-deletion ALT: left as written而不是假装可以处理。*作为 ALT 等位基因表示该样本的等位基因被此处重叠的另一条记录删除了。它不是一个变异把*等位基因计作备选观测会虚增等位基因频率归一化与统计时应显式排除。比较两组变异集之前应该做的事在CHROM:POS:REF:ALT上 join 之前按以下顺序执行# 1. 相同组装、相同 contig 命名 python3 check_contigs.py setA.vcf setB.vcf --genome ref.fa.fai # 2. 结构约定完好 python3 audit_intervals.py setA.vcf --genome ref.fa.fai # 3. 拆分、校验 REF、修剪、左对齐——两组都用同一参考序列 python3 normalize_variant.py --fasta ref.fa --split --input setA.vcf -o A.norm.tsv python3 normalize_variant.py --fasta ref.fa --split --input setB.vcf -o B.norm.tsv只有完成第 3 步之后才能 join。在第 3 步之前计算的交集是未知规模的低估而且是有偏的它会系统性低估重复区段中的 indel——恰恰是大多数有价值的 indel 所在的位置。audit_intervals.py的完整发现规则start_below_one证明 0-based 数据混进了 1-based 文件、many_zero_length证明 1-based 单碱基特征被写进了 0-based 文件、past_contig_end证明组装错误或 contig 边缘 off-by-one、mixed_contig_naming、first_block_offset证明 BED12 的blockStarts被写成了绝对坐标、not_parsimonious警告未修剪等位基因、bad_alt_allele捕获 Ensembl/VEP 的无锚点-记号等可参看 skills/genomic-coordinates/SKILL.md其中not_parsimonious为警告级其余致命项会让退出码变为 1可作为数据目录的 CI 门禁。总结与验证归一化的核心要点可归纳为四条比较之前先归一化任何 join、去重、注释查找如果发生在归一化之前都会静默丢失真实匹配且丢失集中在重复区段里的 indel。归一化 parsimony 修剪 参考序列上的 left-alignment这是 Tan et al. (2015) 的定义bcftools norm与vt normalize遵循同一规则。始终校验 REF归一化基于参考序列REF与 FASTA 不匹配意味着组装版本不同先运行check_contigs.py识别组装不要手工调整坐标。分清边界情况多等位基因先--split再归一化符号等位基因与*不做归一化处理VCF 与 HGVS 方向相反必须通过掌握转录本模型的工具往返转换。整套脚本与夹具在仓库中可离线复现运行uv run --with pytest python -m pytest tests/genomic-coordinates -q即可验证归一化左对齐、幂等性、--compare判定、MISMATCH退出码、--split展开以及组装识别等全部行为参考序列即 tests/genomic-coordinates/fixtures/tiny.fa。上游技能总入口 skills/genomic-coordinates/SKILL.md 与本篇所在的 skills/genomic-coordinates/references/variant-representation.md 相互呼应后者是该技能的等位基因表示、归一化算法、等价性检查与多等位基因拆分的完整参考。【免费下载链接】scientific-agent-skillsTurn any AI agent into an AI Scientist. The #1 Agent Skills library for science, used by 190,000 scientists worldwide. 165 ready-to-use validated skills plus 100 scientific databases covering biology, chemistry, medicine, and drug discovery. Compatible with Cursor, Claude Code, Codex, Pi, Antigravity, and the open Agent Skills standard.项目地址: https://gitcode.com/GitHub_Trending/cl/scientific-agent-skills创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
返回列表