评估 beagle 基因型填充的准确率

小明的数据分析笔记本

最简单的一个思路，只保留vcf文件中不包含任何缺失数据的位点。然后随机把某些样本的部分位点替换成缺失，用beagle做基因型填充，比较填充后和填充前的一致性。

使用 https://www.nature.com/articles/s41586-022-04808-9 这篇论文里SNP的数据

Graph pangenome captures missing heritability and empowers tomato breeding

只用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文件的内容

每列的内容 位点行 样本列 基因型

基因型填充

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文件 上一步输出的三列真实基因型数据，输出数据

输出数据的内容

第三类是真实基因型，第四列是填充后的基因型

统计错误率

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把这个流程穿起来，然后多重复几次统计个平均值

欢迎大家关注我的公众号 小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

微信扫一扫分享文章