大家好,我是邓飞,今天介绍一下基因型数据的填充,包括自填充和填充到参考集上的方法。
对基因型数据进行分析时,无论是GWAS还是GS,都需要对基因型数据进行填充,这里使用Beagle软件说明一下如何对基因型数据进行填充。
第一,定相,Phasing,第二,填充,imputation 根据定向的结构,进行填充。
❝所谓Phasing就是要把一个二倍体(甚至是多倍体)基因组上的等位基因(或者杂合位点),按照其亲本正确地定位到父亲或者母亲的染色体上,最终使得所有来自同一个亲本的等位基因都能够排列在同一条染色体里面。 ❞
LD Phasing是一个非常常用的基因定相方法,它是利用群体中大量无血缘关系的个体,依据基本的连锁不平衡(Linkage disequilibrium,LD)遗传原理和相关数学模型,推断群体中每个个体的单倍体的方法,因此它也是计算量最大的一个。
其它的定相方法:
「软件地址:」https://faculty.washington.edu/browning/beagle/beagle.html
它是一个java程序,所以在windows,Linux,Mac都可以执行,需要提前配置好java环境。
「下载:」
https://faculty.washington.edu/browning/beagle/beagle.22Jul22.46e.jar
mv beagle.22Jul22.46e.jar beagle.jar
说明文档:https://faculty.washington.edu/browning/beagle/beagle_5.4_18Mar22.pdf
「查看帮助文档:」
$ java -jar beagle.jar
beagle.22Jul22.46e.jar (version 5.4)
Copyright (C) 2014-2022 Brian L. Browning
Usage: java -jar beagle.22Jul22.46e.jar [arguments]
data parameters ...
gt=<VCF file with GT FORMAT field> (required)
ref=<bref3 or VCF file with phased genotypes> (optional)
out=<output file prefix> (required)
map=<PLINK map file with cM units> (optional)
chrom=<[chrom] or [chrom]:[start]-[end]> (optional)
excludesamples=<file with 1 sample ID per line> (optional)
excludemarkers=<file with 1 marker ID per line> (optional)
phasing parameters ...
burnin=<max burnin iterations> (default=3)
iterations=<phasing iterations> (default=12)
phase-states=<model states for phasing> (default=280)
imputation parameters ...
impute=<impute ungenotyped markers (true/false)> (default=true)
imp-states=<model states for imputation> (default=1600)
cluster=<max cM in a marker cluster> (default=0.005)
ap=<print posterior allele probabilities> (default=false)
gp=<print posterior genotype probabilities> (default=false)
general parameters ...
ne=<effective population size> (default=100000)
err=<allele mismatch probability> (default: data dependent)
em=<estimate ne and err parameters (true/false)> (default=true)
window=<window length in cM> (default=40.0)
overlap=<window overlap in cM> (default=2.0)
seed=<random seed> (default=-99999)
nthreads=<number of threads> (default: machine dependent)
-Xmx10g,10表示10Gb内存,可以修改 -jar beagle.jar,这个beagle.jar就是我们下载的填充软件
gt=[file],比如gt=a.vcf
,就是被填充文件,vcf格式的文件。如果自填充,这个参数就够了。
如果有参考群,用:ref=[file],比如ref=ref.vcf
,已经被填充好的文件,vcf格式
out=[string],输出名称,比如out=result
,结果会生成result.vcf,被填充好的文件。
map=[file],如果有map文件,单位是cM。如果没有提供,默认是1cM = 1Mb长度。
chrom=[chrom]:[start]-[end],指定染色体和区间,可以不设置
window=[positive number],窗口的大小,如果填充报错,可能区间过大,可以设置大一些
overlap=[positive number],窗口重叠的值,小于窗口数
nthreads=[positive integer],默认是使用所有的线程,也可以设置一个数字
数据:
$ ls
test.map test.ped
「数据质控」
这里,我们想先进行缺失质控,去除缺失率大于10%的位点,然后填充小于10%的位点。
这里,我们使用1-10染色体,缺失率去除大于10%的位点,并将其变为vcf格式
plink --file test --geno 0.1 --chr-set 30 --chr 1-10 --allow-extra-chr --recode vcf-iid --out b --snps-only just-acgt
「结果文件:」
$ ls b.vcf
b.vcf
「填充命令:」
java -Xmx5g -jar beagle.jar window=50 overlap=20 gt=b.vcf out=re_imputed nthreads=2
「结果文件:」
re_imputed.vcf.gz
将其变为plink格式的数据:
plink --vcf re_imputed.vcf.gz --recode --out c1
结果文件:
c1.map c1.ped
首先,参考集是已经填充过的vcf文件,候选集应该是参考集的子集(位点子集)。
「参考集:」
ref.vcf.gz
「测试集:」
test_test.vcf
「填充命令:」代码是空的。参考windows里面的操作
一般,填充需要大的算力,普通电脑内存、算力都不能满足,不建议填充。
但是,填充软件beagle确实可以在windows下运行,下面演示一下。
到官网下载java:https://www.java.com/download/ie_manual.jsp
安装:
在这里插入图片描述
安装好之后,进入windows的cmd终端,键入java
出现下面命令,说明java安装成功。
文件:
进入cmd路径,切换到文件夹中:
windows设置最大内存,报错,就去掉内存设置参数:
java -jar beagle.jar window=50 overlap=20 gt=b.vcf out=re_imputed nthreads=2
结果:
切换到文件夹,进入cmd文件
运行命令:
java -jar beagle.jar window=50 overlap=20 gt=test_test.vcf ref=ref.vcf.gz out=re_imputed nthreads=2
结果:
内存使用情况:
这里,使用测试数据,演示了基因型数据自填充和填充到参考群中的方法。
注意1:填充前,需要将indel位点删除
注意2:填充前,需要先进行缺失质控,因为缺失率大的位点,错误率会很高,所以要先去除
注意3:样本不能有重复,否则填充报错
注意4:填充到参考集时,先要对测试集处理,确保在参考集的子集中