最简单的一个思路，只保留vcf文件中不包含任何缺失数据的位点。然后随机把某些样本的部分位点替换成缺失，用beagle做基因型填充，比较填充后和填充前的一致性。
使用 Nature tomato 这篇论文里SNP的数据
Graph pangenome captures missing heritability and empowers tomato breeding
https://www.nature.com/articles/s41586-022-04808-9
只用3号染色体的数据
首先是把含有缺失基因型的位点过滤掉

time vcftools –gzvcf ../chr3.vcf.gz \
–keep ../332.samples \
–max-missing 1 \
–min-alleles 2 \
–max-alleles 2 \
–recode –recode-INFO-all –out chr3.snp.332

7m42.853s

随机把位点替换成缺失

这里每个位点被选中的概率是10%
每个样本被选中的概率是20%
python randomvcf2NA.py chr3.snp.332.recode.vcf output.snp.vcf truth.snp.sites
randomvcf2NA.py脚本的内容

1 input vcf

2 output vcf

3 truth sites

import sys
import random

fr = open(sys.argv[1],’r’)
fw = open(sys.argv[2],’w’)
fw02 = open(sys.argv[3],’w’)

row = 0
col = 0

for line in fr:
if line.startswith(“#”):
fw.write(line)
else:
row += 1
temp_list = line.strip().split(“\t”)

    if random.random() <= 0.1:
        new_list = []
        for i in range(0,len(temp_list)):
            if i < 9:
                new_list.append(temp_list[i])
            elif i >= 9:
                if random.random() <= 0.2:
                    new_list.append("./.")
                    GT = temp_list[i][0:3]

                    fw02.write("%d\t%d\t%s\n"%(row,i+1,GT))
                else:
                    new_list.append(temp_list[i])

        fw.write('%s\n'%("\t".join(new_list)))



    else:
        fw.write('%s\n'%("\t".join(temp_list)))

fw.close()
fw02.close()
fr.close()
三个位置参数是 输入vcf 输出vcf 和随机选中位点的基因型 truth.sites文件的内容

image.png
每列的内容 位点行 样本列 基因型
基因型填充

time java -Xmx96g -jar \
~/anaconda3/envs/syri/share/beagle-5.2_21Apr21.304-0/beagle.jar \
gt=output.snp.vcf \
out=output.snp.impute \
nthreads=48
这个填充我之前印象里运行是非常慢的。但实际是很快的。一条染色体50万个位点2分钟左右。不知道这个运行很慢的印象是怎么来的了
提取填充过后的基因型

python getImputeSites.py output.snp.impute.vcf.gz truth.snp.sites call.snp.sites
getImputeSites.py脚本的内容是

1 impute vcf gz

2 truth.sites

3 output

import sys
import pandas as pd

df = pd.read_csv(sys.argv[1],comment=”#”,sep=”\t”,header=None)

fw = open(sys.argv[3],’w’)

with open(sys.argv[2],’r’) as fr:
for line in fr:
row = int(line.strip().split()[0])
col = int(line.strip().split()[1])
trueGT = line.strip().split()[2]
imputeGT = df.iloc[row-1,col-1]

    fw.write("%d\t%d\t%s\t%s\n"%(row,col,trueGT,imputeGT))

fw.close()
三个参数 输入填充后的vcf文件 上一步输出的三列真实基因型数据，输出数据
输出数据的内容

image.png
第三类是真实基因型，第四列是填充后的基因型
统计错误率

library(tidyverse)
read_tsv(“call.snp.sites”,col_names = FALSE) %>%
mutate(X5=str_count(X3,’1′),
X6=str_count(X4,’1′)) %>%
mutate(X7=case_when(
X5 == X6 ~ “TRUE”,
TRUE ~ “FALSE”
)) %>%
pull(X7) %>% table() -> x

x

x[1]/sum(x)
可以用snakemake把这个流程穿起来，然后多重复几次统计个平均值
欢迎大家关注我的公众号
小明的数据分析笔记本