git如何合并单个文件_如何利用公共位点合并VCF文件？

欢迎阅读本文，敬请批评指正。

如有问题请邮箱联系(chenshh@shanghaitech.edu.cn)。

VCF(variant calling format)是一种文件格式，可以存储生物信息学研究中多个样本的SNP/indel/SV等信息。

有时候，我们想同时分析来自多个数据集的样本，此时便需要合并存储了样本信息的VCF文件。我们希望合并的文件中包含所有需要的样本，以及这些样本共有的位点。

方法1：bcftools

使用软件bcftools可以达到这一目的。

在合并之前，需要先用bcftools的index功能为VCF建立索引(以下展示的文件均为bgzip压缩过后的VCF)，例如：

 --bash-4.1$ bcftools index A.vcf.gz

需要注意的是，表头前没有“##”注释的VCF是无法建立索引的。这里提供一个注释的例子，仅供参考。

-bash-4.1$ cat vcf.header.txt ##fileformat=VCFv4.1##FILTER=##FORMAT=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=##contig=

方法1.1：先合并，再过滤基因型空缺率高的位点

我们可以首先使用bcftools的merge功能来合并VCF文件。merge会保留原始输入文件中的所有样本和位点。假如某个原始文件中某些样本缺少在另一个原始文件中出现的位点，那么在合并后的文件中，这些样本在这一位点的基因型就是空缺的。例如：

为了去除这些空缺的基因型，可以使用vcftools软件，设置max-missing参数(为一个0到1之间的浮点数，基因型的非空缺率低于该浮点数的位点会被过滤，0表示保留所有空缺基因型，1表示不允许任何缺失基因型出现)。这一参数的具体数值应当根据各个样本文件的具体大小来确定。

比如，如果需要合并两个数据集的VCF，其中一个数据集样本量为2000，另一个数据集样本量为200，我们希望第二个数据集在合并后不包含由于合并造成的空缺基因型。

动用shell自带的bc计算器(or Python，or R，as you like)，

-bash-4.1$ echo "scale = 4; 200/(200+2000)" | bc .0909

我们可以知道，在该数据集中所有样本于某一位点的基因型全部缺失的情况下，整个合并的数据集中，该位点基因型的空缺率(genotype missing rate)至少是0.0909，相应地，非空缺率至多是0.9091，那么过滤缺失基因型的阈值最好设定为0.9091。

如果不考虑具体的样本数量，笼统地设一个阈值，例如0.85，容易导致样本数量较少的数据集中在合并及过滤后出现大量空缺的位点，这样既增大了文件的容量，浪费了存储空间，也容易给后续的数据分析造成一定的困扰。比如，在这样的数据集中随机选取位点，可能选到的大量位点是某些样本并不包含的位点，这些样本自然就难以用于分析了。

综上，我们可以用以下的命令合并VCF文件，然后过滤掉空缺率较高的位点(此处的0.95只是举例)。

####Method1:merge and then filtrate-bash-4.1$ bcftools merge -m id -O z -o merge.vcf.gz A.vcf.gz B.vcf.gz C.vcf.gz -bash-4.1$ vcftools --gzvcf merge.vcf.gz --max-missing 0.95 --recode --out merge.fil -bash-4.1$ bgzip merge.fil.recode.vcf

方法1.2：先取交集，再依据交集合并

我们还可以利用bcftools的isec命令先获取所有数据集中位点的交集，再在每一个数据集的VCF中选取交集包含的位点，然后合并这些含有共同位点的VCF，就像以下的命令：

####Method2: get intersection and then merge-bash-4.1$ bcftools isec A.vcf.gz B.vcf.gz C.vcf.gz -p isec_dir -n =3 -c none --threads 8 -O z-bash-4.1$ bcftools merge ./isec_dir/0000.vcf.gz ./isec_dir/0001.vcf.gz ./isec_dir/0002.vcf.gz -o isec.merge.phased.chrA.vcf.gz -O z -m id

isec的-c参数指定了定义“相同位点”的条件，例如若将参数设置为“all”，则所有位置相同的SNP或indel，无论allele情况如何，都会被视作同一个位点；若设置为“none"，则在allele完全相同的SNP或indel才可以视为同一位点。

通过isec命令的-p参数可以设置含有共同位点的VCF输出的路径。如果不指定输出前缀，输出的VCF会以四位数字编号。

方法2：Python Pandas

如果不用软件，也可以自己动手写脚本解决这个问题。以下是一个可以合并两个VCF文件的Python脚本，取两个VCF文件位点ID的交集进行合并，合并后删去INFO信息，并将所有的FILTER值设为“PASS”。

运行时，通过设置-r1，-r2参数去除VCF文件前的注释行。

import pandas as pd;import argparse;parser = argparse.ArgumentParser()parser.add_argument("-vcf1","--first_vcf",help = "vcf files",\                    required = True,type = str)parser.add_argument("-r1","--rows1",help = "rows to skip in vcf1",\                    required = False,type = int, default = 0)parser.add_argument("-vcf2","--second_vcf",help = "vcf files",\                    required = True,type = str)parser.add_argument("-r2","--rows2",help = "rows to skip in vcf2",\                    required = False,type = int, default = 0)parser.add_argument("-out","--output_vcf",help = "merged vcf to output",\                    required = True,type = str)args = parser.parse_args()def read_fun(fl,skip):    if skip > 0:        vcf = pd.read_csv(fl, sep = '\t', header = 0, skiprows = skip)    elif skip == 0:        vcf = pd.read_csv(fl, sep = '\t', header = 0)    else:        pass    return vcfvcf1 = read_fun(args.first_vcf, args.rows1)vcf2 = read_fun(args.second_vcf, args.rows2)vcf1.index = vcf1['ID']vcf2.index = vcf2['ID']co_id = list(set(vcf1['ID']).intersection(set(vcf2['ID'])))merged_vcf = pd.concat(objs = [vcf1.loc[co_id,:],vcf2.loc[co_id,:].iloc[:,9:]],\                       axis = 1,ignore_index = False)merged_vcf['FILTER'] ='PASS'merged_vcf['INFO'] ='.'merged_vcf.sort_values(by = ["#CHROM","POS"], ascending= True,inplace = True)merged_vcf.to_csv(args.output_vcf,sep = '\t',index = False,header = True,\                  compression = "gzip")

文档读完了！希望对你有所帮助！