最近的工作涉及到处理大量种群数据，其中有个分析需要统计每个个体携带的突变数量，在网上搜寻了诸多方法后确定了一些比较高效的策略，固有此文记录。

本文的部分流程参考源：

BioStars "how to count variants par sample per chromosome in a vcf file?"

提取每个个体携带的突变数量

该流程主要通过 plink2 实现：

1
2
plink2 --vcf [vcf file] --out [plink output1]
plink2 --pfile [plink output1] --sample-counts --out [plink output2]
  • [vcf file]：此处填入 vcf 文件的路径。
  • [plink output1]：此处填入 plink2 输出的文件前缀。
  • [plink output2]：此处填入 plink2 输出的文件前缀。

运行结束后，可以得到以 .scount 结尾的文件，其中各列的含义如下：

  1. #IID：标识每个样本的编号。
  2. HOM_REF_CT：纯合参考基因型计数，个体纯合参考基因型的位点数量。
  3. HOM_ALT_SNP_CT：纯合变异基因型计数，个体纯合变异基因型的单核苷酸多态性（SNP）数量。
  4. HET_SNP_CT：杂合基因型计数，个体在该位点上是杂合基因型的SNP数量（一条染色体是参考基因型，另一条是变异基因型）。
  5. DIPLOID_TRANSITION_CT：二倍体转换变异的数量（转换：嘌呤与嘌呤或嘧啶与嘧啶之间的替换）。
  6. DIPLOID_TRANSVERSION_CT：二倍体颠换变异的数量（颠换：嘌呤与嘧啶之间的替换）。
  7. DIPLOID_NONSNP_NONSYMBOLIC_CT：二倍体非SNP非符号变异计数。包括一些插入、缺失（INDELs）等。
  8. DIPLOID_SINGLETON_CT：二倍体单例变异计数，表示在整个样本集中只出现在该个体中的变异数量。
  9. HAP_REF_INCL_FEMALE_Y_CT：单倍体参考基因型计数，即在单倍体区域（Y 染色体）的参考等位基因数量。
  10. HAP_ALT_INCL_FEMALE_Y_CT：单倍体变异基因型计数，在单倍体区域中出现的变异等位基因的数量。
  11. MISSING_INCL_FEMALE_Y_CT：缺失基因型计数，在样本中缺失的基因型数据计数。

更多可输出列的详细信息可见：https://www.cog-genomics.org/plink/2.0/formats#scount

具体参数的用法可见：https://www.cog-genomics.org/plink/2.0/basic_stats#sample_counts

通过 HOM_ALT_SNP_CT HET_SNP_CT DIPLOID_NONSNP_NONSYMBOLIC_CT 计算每个个体携带的突变数量。

提取每个突变对应的携带者信息

1
bcftools norm -m -any [vcf file] | bcftools query -i 'GT="alt"' -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE:%GT]\n' - > [output]
  • [vcf file]：此处填入 vcf 文件的路径。
  • [output]：此处填入输出文件的路径。

此处通过 bcftools norm -m -any 对多等位基因进行了拆分，可视自身需求去除该部分，此处 GT="alt" 将所有携带替代等位基因的个体提取出并打印。

输出结果示例：

1
chr1    228424930       G       A       4922792:0/1     4884261:0/1     3988378:0/1     2479308:0/1     4414951:0/1     3862053:0/1

前四列用于表示特定突变（分别指示染色体编号、位置、ref 和 alt），后续以制表符分隔记录该突变的携带者及基因型。