Vcftools上堰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

如果可以将数据转换为等级统计格式，则可以从varcomp.glob中获取方差组件。我通常做的是：

1. 用vcftools结合--012获得基因型
2. 将0/1/2/-1转换为象形文字格式(例如，11/12/22/NA)
3. 将数据加载到递阶统计和计算中(见下文)

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 --我有大量的snps，但是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")

然后，你会把文件hierfstat.txt读入R，这些是你的位点。您需要在抽样设计中指定您的级别(例如，您的人口)。然后调用varcomp.glob()来获取方差组件。如果您想使用这个这里，我有一个并行版本。

请注意，在本例中，您指定0作为参考等位基因。也许是你想要的，也许不是。我经常计算小等位基因的频率，并使2个小等位基因，但这取决于你的学习目标。