最简单的一个思路，只保留vcf文件中不包含任何缺失数据的位点。然后随机把某些样本的部分位点替换成缺失，用beagle做基因型填充，比较填充后和填充前的一致性。
使用 Nature tomato 这篇论文里SNP的数据
Graph pangenome captures missing heritability and empowers tomato breeding
https://www.nature.com/articles/s41586-022-04808-9
只用3号染色体的数据
首先是把含有缺失基因型的位点过滤掉
代码语言：javascript复制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%
代码语言：javascript复制python randomvcf2NA.py chr3.snp.332.recode.vcf output.snp.vcf truth.snp.sites
randomvcf2NA.py脚本的内容
代码语言：javascript复制# 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("%dt%dt%sn"%(row,i 1,GT))
else:
new_list.append(temp_list[i])
fw.write('%sn'%("t".join(new_list)))
else:
fw.write('%sn'%("t".join(temp_list)))
fw.close()
fw02.close()
fr.close()
三个位置参数是 输入vcf 输出vcf 和随机选中位点的基因型 truth.sites文件的内容
image.png
每列的内容 位点行 样本列 基因型
基因型填充
代码语言：javascript复制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分钟左右。不知道这个运行很慢的印象是怎么来的了
提取填充过后的基因型
代码语言：javascript复制python getImputeSites.py output.snp.impute.vcf.gz truth.snp.sites call.snp.sites
getImputeSites.py脚本的内容是
代码语言：javascript复制# 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("%dt%dt%st%sn"%(row,col,trueGT,imputeGT))
fw.close()
三个参数 输入填充后的vcf文件 上一步输出的三列真实基因型数据，输出数据
输出数据的内容
image.png
第三类是真实基因型，第四列是填充后的基因型
统计错误率
代码语言：javascript复制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把这个流程穿起来，然后多重复几次统计个平均值
