从vcftools上的Weir-fst推断差异成分

Question

vcftools --vcf ALL.chr1.phase3_shapeit2_mvncall_integrated_v5.20130502.genotypes.vcf --weir-fst-pop POP1.txt --weir-fst-pop POP2.txt --out fst.POP1.POP2

上面的脚本使用Weir和Cokerham的1984年公式计算1000个基因组总体数据的FST距离。该公式使用3个方差成分，即A，B，C（种群之间；在人群中的个体之间；人群中的个体之间的配子之间）。

输出直接提供了公式的结果，而不是该程序计算得出的组件以达到最终结果。如何要求VCFTools输出A，B，C？

Answer 1

如果您可以将数据输入HIERFSTAT的格式，则可以从varcomp.glob获取差异组件。我通常要做的是：

1. 使用--012使用vcftools获取基因型
2. 将0/1/2/-1转换为hierfstat格式（例如，11/12/22/na）
3. 将数据加载到hierfstat并计算（见下文）

r示例：

library(hierfstat)
data = read.table("hierfstat.txt", header=T, sep="t")
levels = data.frame(data$popid)
loci = data[,2:ncol(data)]
res = varcomp.glob(levels=levels, loci=loci, diploid=T)
print(res$loc)
print(res$F)

Fst对于每个基因座（行）是（无层次设计），来自res$loc：res$loc[1]/sum(res$loc)。如果您的采样更复杂，则需要以不同的方式解释方差组件。

- 根据您的评论更新 -

我在熊猫中这样做，但是任何语言都会做。这是文本替代练习。只需将您的.012文件放入数据框架中，然后转换如下即可。我逐行阅读numpy b/c我有大量的SNP，但是Read_CSV也可以。

import pandas as pd
import numpy as np
z12_data = []
for i, line in enumerate(open(z12_file)):
    line = line.strip()
    line = [int(x) for x in line.split("t")]
    z12_data.append(np.array(line))
    if i % 10 == 0:
        print i
z12_data = np.array(z12_data)
z12_df = pd.DataFrame(z12_data)
z12_df = z12_df.drop(0, axis=1)
z12_df.columns = pd.Series(z12_df.columns)-1
hierf_trans = {0:11, 1:12, 2:22, -1:'NA'}
def apply_hierf_trans(series):
    return [hierf_trans[x] if x in hierf_trans else x for x in series]
hierf = df.apply(apply_hierf_trans)
hierf.to_csv("hierfstat.txt", header=True, index=False, sep="t")

然后，您会在R中读取该文件hierfstat.txt，这些是您的基因座。您需要在抽样设计（例如人口）中指定您的水平。然后致电varcomp.glob（）获取方差组件。如果您想使用它，我在这里有一个并行版本。

请注意，在这种情况下，您将0指定为参考等位基因。可能是您想要的，也许不是。我经常计算次要等位基因频率并使2个小等位基因，但这取决于您的学习目标。