如何对比hg38与灵长类共同祖先fasta并生成差异VCF文件？

针对你用hg38和灵长类共同祖先基因组生成人类特异性SNP的VCF需求，有两个成熟的工具链可以实现：

1. BWA + Samtools + BCFtools 标准流程 #

这是基因组变异检测的通用方案，适合高精度的两两基因组对比：

• 步骤1：预处理基因组文件
  先解压祖先基因组压缩包，再为hg38构建BWA索引：

  tar -xzf homo_sapiens_ancestor_GRCh38.tar.gz
  bwa index hg38.fa
  

• 步骤2：基因组比对
  将祖先基因组序列比对到hg38参考序列上：

  bwa mem -t 8 hg38.fa homo_sapiens_ancestor_GRCh38.fa > ancestor_vs_hg38.sam
  

  -t 8代表使用8个线程，可根据硬件资源调整数值。

• 步骤3：SAM文件处理
  把SAM转为BAM，排序并去除PCR重复：

  samtools view -Sb ancestor_vs_hg38.sam > ancestor_vs_hg38.bam
  samtools sort ancestor_vs_hg38.bam -o ancestor_vs_hg38_sorted.bam
  samtools markdup -r ancestor_vs_hg38_sorted.bam ancestor_vs_hg38_sorted_dedup.bam
  samtools index ancestor_vs_hg38_sorted_dedup.bam
  

• 步骤4：调用SNP生成VCF
  用bcftools检测变异并输出VCF：

  bcftools mpileup -f hg38.fa ancestor_vs_hg38_sorted_dedup.bam | bcftools call -mv -Ov -o ancestor_vs_hg38_snps.vcf
  

  参数说明：

  • -mv：仅调用SNP（排除插入缺失）并启用多等位基因检测
  • -Ov：输出文本格式的VCF文件

2. MUMmer 快速比对方案 #

MUMmer是专为近缘基因组比对设计的工具，速度快，适合快速筛选差异位点：

• 步骤1：执行基因组比对

  nucmer hg38.fa homo_sapiens_ancestor_GRCh38.fa
  

• 步骤2：提取SNP并转VCF
  先提取差异位点到文本文件，再转换为标准VCF：

  show-snps -Clr out.delta > snps.txt
  

  snps.txt包含了所有差异位点的详细信息，你可以用Python/R脚本自行转换为VCF格式，也可以用mummer2vcf工具完成格式转换。

关键注意事项 #

• 确保两个基因组的染色体命名规则一致（比如hg38的chr1和祖先基因组的1），不一致时需先统一命名，否则比对结果会出现大量错位。
• 要筛选人类特异性SNP，可以用bcftools过滤掉祖先基因组与人类基因组一致的位点：

  bcftools filter -i 'ALT != REF' ancestor_vs_hg38_snps.vcf -o human_specific_snps.vcf
  

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