2 
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 的值？

4

1 回答 1

4

如果您可以将数据转换为hierfstat的格式，则可以从中获取方差分量varcomp.glob。我通常做的是：

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

示例：

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)。如果您有更复杂的抽样，则需要以不同方式解释方差分量。

--根据您的评论更新--

我在 Pandas 中这样做，但任何语言都可以。这是一个文本替换练习。只需将您的 .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 设为次要等位基因，但这取决于您的研究目标。

于 2015-05-14T20:57:27.213 回答