Экстраполирующие дисперсионные компоненты от Weir-Fst на Vcftools

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

Приведенный выше скрипт вычисляет расстояния Fst на 1000 популяционных данных генома, используя формулу Вейра и Кокерхэма 1984 года. Эта формула использует 3 компонента дисперсии, а именно, a, b, c (между популяциями; между индивидуумами внутри популяций; между гаметами внутри индивидуумов внутри популяций).

Выходные данные непосредственно предоставляют результат формулы, но не компоненты, которые программа рассчитала для достижения окончательного результата. Как я могу попросить Vcftools вывести значения для a, b, c?

Answer 1

Если вы можете получить данные в формате для hierfstat, вы можете получить компоненты дисперсии от varcomp.glob, Что я обычно делаю, это:

1. использование vcftools с --012 получить генотипы
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 у меня есть тонны 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 минорными аллелями, но это зависит от цели вашего исследования.