基于大体积GVCF文件开展群体遗传学分析的技术咨询

第一步：预处理GVCF，生成精简的变异位点VCF #

GVCF包含大量非变异位点，直接转换会导致文件臃肿，R工具无法处理。必须先用NGS工具完成这一步：

• 用GATK合并样本GVCF并调用变异，生成只含变异位点的VCF：

  gatk GenotypeGVCFs -R reference.fasta -V gendb://gvcf_db -O raw_variants.vcf.gz
  

• 对原始变异做硬过滤，保留高质量SNP：

  gatk VariantFiltration -V raw_variants.vcf.gz -O filtered_variants.vcf.gz \
    --filter-expression "QUAL < 30.0 || QD < 2.0 || FS > 60.0 || SOR > 3.0 || MQ < 40.0" \
    --filter-name "LowQualSNP"
  

• 提取纯SNP位点（排除Indel）：

  gatk SelectVariants -V filtered_variants.vcf.gz -O snp_only.vcf.gz --select-type-to-include SNP
  

第二步：将过滤后的VCF转换为adegenet兼容格式 #

处理后的VCF体积大幅缩小，可选择以下两种稳定路径：

路径1：用PLINK转Genepop格式 #

PLINK处理大VCF效率远高于R，适合批量转换：

• 先转PLINK二进制格式：

  plink --vcf snp_only.vcf.gz --make-bed --out plink_data --allow-extra-chr
  

• 再转Genepop格式：

  plink --bfile plink_data --recode genepop --out genepop_data --allow-extra-chr
  

  生成的genepop_data.gen可直接用adegenet读取：

  library(adegenet)
  genepop_obj <- read.genepop("genepop_data.gen", ncode=3)
  

路径2：用bcftools+PGDSpider转多格式 #

如果需要GENETIX/STRUCTURE/FSTAT格式，用bcftools导出纯文本VCF后，用PGDSpider做格式转换：

• 导出纯文本VCF：

  bcftools view snp_only.vcf.gz -O v -o snp_only_plain.vcf
  

• 打开PGDSpider，选择输入格式为VCF，输出格式选对应目标格式（GENETIX/STRUCTURE/FSTAT），按向导完成转换，生成的文件用adegenet对应函数读取：
  • GENETIX：read.genetix()
  • STRUCTURE：read.structure()
  • FSTAT：read.fstat()

第三步：R中高效处理大群体遗传数据 #

如果转换后的数据仍较大，可在adegenet中做以下优化：

• 转换为genpop对象减少内存占用：

  genpop_obj <- genind2genpop(genepop_obj)
  

• 计算统计量时优先用向量化操作，必要时用parallel包并行计算。

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