倍体>2时bcftools merge等位基因拆分及合并报错求助

核心报错原因 #

1. 旧索引文件不匹配：脚本直接复制原始VCF的.tbi索引，但修改后的VCF内容已变更（移除INFO字段），旧索引无法匹配新文件，导致bcftools merge解析索引失败。
2. INFO字段未完全移除：报错行仍存在AC=1;AF=0.333等INFO字段，说明bcftools annotate -x INFO命令可能未正确执行，或文件路径存在问题导致未处理目标文件。
3. 多倍体FORMAT字段残留：即使移除INFO，FORMAT中的GT等字段仍保留多倍体信息（如0/0/1），bcftools merge仍可能因倍体不兼容报错。

修正步骤 #

1. 修复脚本：重新生成索引而非复制旧索引 #

修改脚本中处理索引的部分，删除复制旧tbi的代码，改用bcftools index为修改后的VCF重新建立索引：

#!/bin/bash

cd /nesi/nobackup/vuw03939/combined_vcf/
pwd

vcf_dir="../nextflow/liti_3iploid/sarek_anno_3ploid_snpeff/annotation"
output_dir="./combined_3ploid"
mkdir -p "$output_dir"
file_list="3ploid.list"

echo "Files to process:"
cat "$file_list"

while IFS= read -r vcf_file || [ -n "$vcf_file" ]; do
    filename=$(basename "$vcf_file")
    output_file="${filename%.vcf.gz}_removed.vcf.gz"
    # 移除INFO和FORMAT中与多倍体相关的字段（保留GT用于后续处理）
    bcftools annotate -x INFO,FORMAT/PL,FORMAT/AD -O z -o "${output_dir}/${output_file}" "${vcf_dir}/${vcf_file}"
    # 为新生成的VCF建立索引
    bcftools index "${output_dir}/${output_file}"
done < "$file_list"

• 额外移除FORMAT/PL和FORMAT/AD：这些字段与基因型剂量相关，多倍体的格式会干扰bcftools merge
• 用bcftools index生成匹配的新索引，确保文件与索引一致性

2. 验证INFO字段是否完全移除 #

处理完成后，随机抽取一个修改后的VCF文件验证：

# 查看前10行数据行（跳过表头）
bcftools view -H "${output_dir}/your_sample_removed.vcf.gz" | head -10

确认数据行的第5列（INFO列）为.，无残留的AC=...等字段。

3. 调整合并命令：指定参考基因组与倍体兼容参数 #

执行合并时，添加--force-samples和--ploidy参数（即使移除字段，明确指定倍体可避免潜在检测逻辑报错）：

bcftools merge --force-samples --ploidy 3 combined_3ploid/*.gz -Ov -o combined_3ploid.vcf

• --force-samples：强制合并样本名重复的文件（如果存在）
• --ploidy 3：明确告知工具处理3倍体数据（若有4/5倍体，可统一设为最高倍体，或分批次处理后再合并）

4. 替代方案：转换为单倍型拆分后合并 #

如果上述方法仍失败，可将多倍体VCF拆分为多个单倍型VCF后再合并：

# 拆分单个3倍体VCF为3个单倍型文件
bcftools norm -m - "${vcf_dir}/sample.vcf.gz" -O z -o sample_hap1.vcf.gz
# 对所有样本执行拆分后，再用bcftools merge合并所有单倍型文件

• -m -：将多倍体基因型拆分为单倍型（每个单倍型作为独立样本）
• 合并后可在VEP中查看所有单倍型的SNP信息

额外注意事项 #

• 确保file_list中的文件路径正确：若列表内是相对路径，需与vcf_dir拼接后指向真实文件，避免处理空文件或错误文件
• 检查VCF表头：部分注释工具会添加##contig等信息，合并前可用bcftools norm -f ref.fasta统一参考基因组序列信息

内容的提问来源于stack exchange，提问作者Isaac