从 VCF 到完整突变转录本与蛋白：mutatedSequence 2.0

过去一年里，我把 2023 年的 mutatedSequence 流程整个重写了一遍。旧的 bcftools/seqkit 版本能跑，但它对负链、跨外显子变异、partial model 和 SV 的处理都不够可靠。新版本改成 genome-first block 模型：先把变异应用到基因组坐标，再按 exon/CDS 切取，最后翻译成蛋白。

这篇文章记录三个最容易踩的坑、工具的设计、输出，以及哪些情况它明确处理不了。代码和完整文档都在 GitHub: mutatedSequence。

7 分 14 秒 · 中文女声旁白 · 1080p · Remotion

为什么重写

旧流程把每个 CDS exon 单独写成 FASTA，再用 bcftools consensus 逐个应用变异。这个设计有几个根本问题：

负链：CDS 已经反向互补，但 VCF 的 REF/ALT 是 forward strand，REF 自然对不上；
跨 exon junction：VCF 的 REF 可能包含 intron，而 per-exon FASTA 没有这段，变异无法应用；
多等位/phasing：一个 consensus 只能选一个 ALT，bcftools 默认也不区分 haplotype；
符号型 SV：bcftools consensus 只支持 <DEL>、<*>、<NON_REF>，不支持 <INS>/<DUP>/<INV>；
3' 延伸：旧流程在 CDS 后接 5 kb 基因组序列，会把 intron 也翻译进去。

新版本先构建 genomic blocks，再按 exon/CDS 取序列。这样 deletion 跨 intron、整个 exon 被删除、负链、5' N pad 和 split stop 都能在同一个模型里处理。

三个常见问题

1. REF 错位，参考序列前面多 1–3 个 X/N？

旧代码里有一行 mode = b.replace(a, 'X')：a 是从基因组切出的 CDS，b 是 Ensembl CDS FASTA。这里的 X 只是占位符，不是碱基，也不是 alignment。

Ensembl GTF 的 CDS 行不包含终止密码子，终止密码子单独写成 stop_codon；而 cds.all.fa 包含 stop。所以 GTF CDS 拼接结果通常比 Ensembl CDS 短 3 bp。
5' 不完整模型（mRNA_start_NF/cds_start_NF）在 cds.all.fa 里会用 1–2 个 N 补齐第一个密码子，而 GTF CDS 从第一个已知碱基开始。
因此 b = "NN" + a + "TAG" 时，mode = "NNXTAG"：2 个 N 是 5' padding，X 是匹配部分，TAG 是 3' stop。看起来像“前面多了一个 X”，实际上多出的 3 个碱基在 3' 端。

真实 chr21 例子：

Transcript	Raw GTF CDS	Ensembl CDS	差值	解释
ENST00000432178 (CCT8)	319	324	+5	5' 2 个 N + 3' stop
ENST00000460679 (JAM2)	782	786	+4	5' 1 个 N + 3' stop
ENST00000284885 (TMPRSS15)	3057	3060	+3	无 N；仅 3' stop
ENST00000381284 (ITSN1)	949	951	+2	5' 2 个 N；无 stop

chr21 全部 1,050 条 CDS 的 raw length difference：

+3: 831   0: 136   +5: 40   +4: 29   +1: 8   +2: 6
leading N count: 0: 967, 1: 37, 2: 46
terminal codon: stop: 900, no_stop: 150

如果你在 VCF 层面看到 REF 对不上，还要检查：负链是否反向互补、indel 的 POS 是 anchor（REF=AT, ALT=A 删除的是 POS+1）、VCF 是否 left-align 过、大小写（dna_sm 是小写）、assembly/contig 是否一致。新工具的 diagnose 会输出这些分布；--strict 可以让 REF mismatch 直接报错。

2. 为什么很多转录本没有 stop codon？

Ensembl 109 的 121,766 条 CDS FASTA 中，99,812 条（81.97%）以 TAA/TAG/TGA 结尾，21,954 条（18.03%）不是。这不是 bug：

18,559 条带 mRNA_end_NF：3' 不完整，stop 在已知转录本之外；
2,614 条带 mRNA_start_NF：5' 不完整（可与 end_NF 重叠）；
约 574 个 IG/TR V/D/J 基因片段，本来就不是完整蛋白；
114 个 non_stop_decay，无 stop，被 NSD 降解；
还有 758 条以 stop triplet 结尾，但 CDS 长度不是 3 的倍数，即不在阅读框内。

严格无框内 stop 的集合是 22,712（21,954 + 758）。旧流程的 isStopWithCodon 强制 len(seq) % 3 == 0，会排除 758 个伪 stop，但把 partial model 都标成 non-stop；tailCdsSeq 还使用了错误的负链坐标。新工具给每条转录本输出 reference_class（complete、partial_3、partial_5_complete_3 等）和 reference_stop_reason，不强行要求 stop，并同时输出 annotated protein 和 3' 延伸蛋白。

3. 等位基因和 SV 要注意什么？

不要盲目应用所有 ALT；按 GT/phase 选。支持 alt、genotype、hom、hap1、hap2。
未 phasing 的 0/1 用 hap1/hap2 只是按 GT 顺序约定，结果是 mosaic，不是真实单倍型。
1/2 不能在一个线性 consensus 里同时表示两个 ALT；分别跑两次。
重叠变异默认跳过后来的，TSV 标记 overlap_skipped；--overlap error 直接报错。
符号型 <DEL> 用 POS anchor + END；<INS> 用 LEFT/RIGHT_SVINSSEQ（建议只取一个）；<DUP> 按 tandem 处理；<INV> 反转 POS+1..END。
<TRA>/<BND>/<CNV> 无法表示成一条线性 transcript，会标记 unresolved。

dysgu issue #65 的核心就是这件事：bcftools consensus 不支持符号型 SV，SVLEN/END/SVINSSEQ 又经常不一致，phasing 还会改变结果。新工具会写一个 normalized.vcf（单 ALT、显式 REF/ALT），方便交叉检查，但主输出不依赖 bcftools。

设计：genome-first blocks

genome FASTA + GTF/GFF3 + VCF
        ↓
解析 transcript / exon / CDS / stop_codon
        ↓
REF 与 forward genome 校验（--strict 可报错）
        ↓
选择 allele → 构建 Edit → 处理 unresolved/overlap
        ↓
按 chromosome 生成 reference/ALT blocks
        ↓
逐个 transcript 切 exon / CDS
        ├─ mutated cDNA
        ├─ mutated CDS（CDS + stop_codon + 5' N pad）
        ├─ annotated protein（停在注释 stop）
        └─ extended protein（CDS + 剪接 3' UTR）
        ↓
per-variant + aggregate TSV / FASTA / normalized VCF

核心是 block 列表：reference block、alternate block、insertion block。exon/CDS 和 block 求交，因此 deletion 跨 intron 时只删除 exonic 部分；整个 exon 被删除时，剩余 exon 正确拼接；负链先取 plus-strand block，再对 exon/CDS 整体反向互补。

输出

每次 run 的主要文件：

mutation_report.tsv：每个 variant × transcript 一行，含原始 VCF REF/ALT、选中 allele、REF 校验、codon/AA 变化，以及完整 before/after cDNA/CDS/protein；
transcripts.tsv：每个 transcript 一行的聚合结果；
variants.tsv：每个 VCF 记录/选中 ALT 一行；
normalized.vcf：单 ALT、显式 REF/ALT；
sequences/：transcript/CDS/protein/基因组区域 FASTA。

示例（测试数据）：

variant	transcript	consequence	codon	AA	protein diff
A_missense	TX_A	missense	CTG → ATG	L → M	p.Leu11Met
A_nonsense	TX_A	nonsense	AAA → TAA	K → *	p.Lys30_Asn233del
A_frameshift_ins	TX_A	frameshift	TAT → TCA	Y → S	p.Tyr41_Asn233delins3aa
B_missense	TX_B	missense	CTG → ATG	L → M	p.Leu11Met

B_missense 是负链：VCF ALT 是 forward 的 T，转录本上反向互补成 ATG。

明确不能处理的情况

这部分很重要。遇到以下情况，不要假设工具“默默做对了”：

情况	行为	建议
未 phasing 的 compound/overlapping variants	生成 mosaic；后一个 edit 被 skip	先 phasing；用 Haplosaurus/VCF2Prot 做 phased haplotype
splice-site 导致 exon skipping、cryptic splice、intron retention	只注释 splice region，exon 结构固定	VEP + SpliceAI/MMSplice，或 alternative transcript
<TRA>/<BND>/<CNV>、复杂重排	unresolved + skip	assembly/pangenome、AnnotSV、Varcode、vcf2diploid
没有 SVINSSEQ 的 <INS>、非 tandem <DUP>	无法重建/按 tandem 近似	assembly-based insertion，或重新表示 VCF
线粒体密码、硒蛋白、alternative start、programmed frameshift	使用标准核密码；CTG 翻译成 L，不强制 M	专用遗传密码/注释工具
assembly/release 不一致、transcript 坐标 VCF	大量 REF mismatch；warning/skip 或 --strict 报错	统一 reference；liftover；bcftools norm
BCF/MAF/gVCF、BAM/RNA、HLA	不接受/不使用	转 VCF；Isovar；pVACseq
exon 边界 insertion、partial-overlap replacement block	边界包含是约定；ALT 线性切片近似	检查边界；拆成更小的显式 edit

另外，未 phasing 的多样本、polyploid、somatic VAF/subclone、全基因组 mutated FASTA、超大 TSV、临床/群体/regulatory 注释、HGVS、neoantigen 都不是它的强项。

和现有工具的关系

它不是 VEP 的替代品，也不是“全面更好”。更准确地说：

需求	更合适的工具
临床/群体/regulatory/domain 注释、SIFT/PolyPhen、HGVS	VEP、SnpEff、Mutalyzer
phased CDS/protein haplotype、compound variants	Haplosaurus、VCF2Prot
极快的 VCF consequence 注释	bcftools csq
新抗原/HLA	pVACseq
RNA reads 证据	Isovar
完整 cDNA/CDS/annotated+extended protein、SV/partial model/审计 TSV	mutatedSequence

我们的优势主要在 sequence-level：完整 cDNA/CDS、5' N pad、split stop、负链、跨 exon junction、常见 SV、逐 variant + aggregate 审计表。劣势在注释广度、真实 phasing、splicing 建模和超大规模。完整的对比（VEP/Haplosaurus、bcftools csq、SnpEff、VCF2Prot、ORFan、Varcode、Isovar、Mutalyzer/TransVar、pVACseq、vcf2diploid）写在 README 里。

怎么运行

# 直接用 Python（无第三方依赖）
python -m mutatedseq run \
  -g genome.fa -a annotation.gtf -v sample.vcf -o out \
  --cds cds.all.fa.gz --cdna cdna.all.fa.gz --pep pep.all.fa.gz \
  --sample S1

# 诊断 REF 偏移 / N pad / 非终止转录本
python -m mutatedseq diagnose \
  -g genome.fa -a annotation.gtf -v sample.vcf -o diag \
  --cds cds.all.fa.gz --pep pep.all.fa.gz

# REF 不匹配直接报错
python -m mutatedseq run ... --strict

镜像已发布到 Docker Hub：wenyuhaobuemoji/mutatedsequence。镜像只包含环境，代码通过 -v 挂载：

docker pull wenyuhaobuemoji/mutatedsequence:2.0.0
docker run -it --rm \
  -v "${PWD}:/workspace" -w /workspace \
  wenyuhaobuemoji/mutatedsequence:2.0.0

# 容器内
python -m mutatedseq --version
python -m mutatedseq run ...

验证

合成测试数据覆盖正负链、跨外显子 deletion、partial model、split stop、MNV、frameshift/in-frame、multiallelic/phased、overlap、符号型 SV 和 REF mismatch；10 个端到端测试全部通过。真实 Ensembl 109 chr21：1,050 条 CDS 重建全部一致；952/1,050 个蛋白与 Ensembl pep 完全一致，98 个差异全部是 3' 不完整模型的末端 partial codon。

链接

如果这篇文章只能留下一句话：把 VCF 变成完整蛋白，真正的难点不是“翻译”，而是坐标、等位基因、转录本结构和 SV 的边界。