如何从新BAM文件获取指定SNV的基因型数据（REF/ALT覆盖度）

前提准备 #

1. 把你手里的SNV位点整理成标准格式文件，推荐用VCF格式（命名为target_sites.vcf），内容如下：

##fileformat=VCFv4.2
#CHROM  POS ID REF ALT QUAL FILTER INFO
chr1   1000  .   A   C   .   .   .
chr1   2000  .   T   A   .   .   .

或者用BED格式（注意BED是0起始坐标，POS要减1，命名为target_sites.bed）：

chr1    999 1000    A/C
chr1    1999 2000    T/A

2. 确保你有对应参考基因组序列文件（reference.fasta），且BAM文件已建立索引（tumor1.bam.bai、tumor2.bam.bai）。

方法1：用bcftools（快速高效，推荐） #

bcftools是处理VCF/BAM的轻量级工具，能直接生成符合需求的格式：

1. 用mpileup扫描目标位点，结合call生成原始变异文件，再添加DR/DV字段：

bcftools mpileup -f reference.fasta -R target_sites.vcf tumor1.bam tumor2.bam | \
bcftools call -mv -Ov | \
bcftools +fill-tags -Ov -- -t DR,DV > final_output.vcf

• -R target_sites.vcf：指定只分析目标位点
• +fill-tags -t DR,DV：自动添加参考等位基因深度（DR）和变异等位基因深度（DV）字段
• 输出的final_output.vcf会包含你需要的GT:DR:DV格式信息。

2. 若需要更定制化的输出格式，可使用bcftools query提取字段：

bcftools mpileup -f reference.fasta -R target_sites.vcf tumor1.bam tumor2.bam | \
bcftools call -mv -Ov | \
bcftools query -f '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="Number of reference reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="Number of variant reads">
#CHROM  POS ID REF ALT FORMAT tumor1.bam tumor2.bam
%CHROM\t%POS\t%ID\t%REF\t%ALT\tGT:DR:DV\t[%GT:%DR:%DV]\t[%GT:%DR:%DV]\n' > custom_output.vcf

方法2：用GATK（严谨流程，适合临床/科研场景） #

GATK是基因组分析的标准工具，基因型分型更准确：

1. 准备目标区间文件（target.intervals）：

chr1:1000-1000
chr1:2000-2000

2. 生成单个样本的gVCF：

gatk HaplotypeCaller -R reference.fasta -I tumor1.bam -L target.intervals -ERC GVCF -O tumor1.g.vcf.gz
gatk HaplotypeCaller -R reference.fasta -I tumor2.bam -L target.intervals -ERC GVCF -O tumor2.g.vcf.gz

3. 合并gVCF并进行基因型分型：

gatk GenomicsDBImport -V tumor1.g.vcf.gz -V tumor2.g.vcf.gz -L target.intervals --genomicsdb-workspace-path gdb_workspace
gatk GenotypeGVCFs -R reference.fasta -V gendb://gdb_workspace -O raw_variants.vcf.gz

4. 添加等位基因深度注释（DR/DV对应DepthPerAlleleBySample的结果）：

gatk VariantAnnotator -R reference.fasta -V raw_variants.vcf.gz -I tumor1.bam -I tumor2.bam -O annotated_variants.vcf.gz -A DepthPerAlleleBySample

之后可通过bcftools提取需要的FORMAT字段，转换成GT:DR:DV格式。

方法3：用samtools+脚本（灵活定制） #

如果需要完全自定义统计逻辑，可先用samtools输出原始比对数据，再用脚本处理：

1. 生成pileup文件：

samtools mpileup -f reference.fasta -l target_sites.bed tumor1.bam tumor2.bam > pileup.txt

2. 用Python脚本统计REF/ALT覆盖度并生成目标格式：

import pandas as pd

# 读取pileup数据
pileup_df = pd.read_csv("pileup.txt", sep="\t", header=None, 
                        names=["CHROM", "POS", "REF", "t1_bases", "t1_quals", "t2_bases", "t2_quals"])
# 读取目标位点信息
target_df = pd.read_csv("target_sites.vcf", sep="\t", comment="#", 
                        names=["CHROM", "POS", "ID", "REF", "ALT"])

# 定义统计函数
def count_alleles(bases_str, ref_base, alt_base):
    ref_count = bases_str.count(ref_base.upper()) + bases_str.count(ref_base.lower())
    alt_count = bases_str.count(alt_base.upper()) + bases_str.count(alt_base.lower())
    return ref_count, alt_count

# 合并数据并统计
merged_df = pd.merge(pileup_df, target_df, on=["CHROM", "POS", "REF"])
merged_df[["t1_dr", "t1_dv"]] = merged_df.apply(
    lambda x: count_alleles(x["t1_bases"], x["REF"], x["ALT"]), axis=1, result_type="expand"
)
merged_df[["t2_dr", "t2_dv"]] = merged_df.apply(
    lambda x: count_alleles(x["t2_bases"], x["REF"], x["ALT"]), axis=1, result_type="expand"
)

# 生成目标VCF格式
output_df = merged_df[["CHROM", "POS", "ID", "REF", "ALT"]].copy()
output_df["FORMAT"] = "GT:DR:DV"
output_df["tumor1.bam"] = "0/1:" + merged_df["t1_dr"].astype(str) + ":" + merged_df["t1_dv"].astype(str)
output_df["tumor2.bam"] = "0/1:" + merged_df["t2_dr"].astype(str) + ":" + merged_df["t2_dv"].astype(str)

# 写入文件
with open("custom_result.vcf", "w") as f:
    f.write("##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")
    f.write("##FORMAT=<ID=DR,Number=1,Type=Integer,Description=\"Number of reference reads\">\n")
    f.write("##FORMAT=<ID=DV,Number=1,Type=Integer,Description=\"Number of variant reads\">\n")
    f.write("#CHROM\tPOS\tID\tREF\tALT\tFORMAT\ttumor1.bam\ttumor2.bam\n")
    output_df.to_csv(f, sep="\t", index=False, header=False)

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