记录---提取合并VCF文件

记录—提取合并VCF文件

最近有个需求，需要合并两个VCF集合，两个文件中材料名完全不同，SNP名部分相同，想要合并为相同SNP下不同Sample的一个VCF文件。

整体思路是：
（1）找到两个VCF文件的共有SNP
（2）合并两个VCF文件（SNP相同，Sample列不同）

0. 简化GATK结果生成的VCF文件

生成的GATK的VCF文件中包含很多信息，文件特别大，想要简化一下，保留基因型信息，剔除不想要的信息：

import os,sysinf=sys.argv[1] # input file
ouf=sys.argv[2] # output filewith open(inf) as infile, open(ouf,"w") as oufile:for line in infile:if line.startswith("##fileformat"):oufile.write(line)elif line.startswith("#CHROM"):oufile.write(line)else:num=["Chr"+str(i) for i in list(range(1,11))]a=line.strip().split()if a[0] in num:gt=list(map(lambda x: x.split(":")[0], a[9:]))oufile.write("%s\t%s\t%s\t%s\t%s\t%s\n"\%(a[0],a[1],"S"+a[0]+"_"+a[1], "\t".join(a[3:7] + ["."]),"GT","\t".join(gt) ))

1. 两个VCF基因集合的操作和合并

一个VCF的SNP数据是自测序结果（A.vcf），另一个是下载的HapMap数据结果(B.vcf)。因为A1的SNP数量少，所以提取A的SNP作为参考，通过VCFtools提取B.vcf，再将结果的SNP提取出来：

（1）第一步先使用grep和awk等提取A.vcf中的SNP(A.vcf标记少)

grep -v "#" A.vcf | cut -f 3 > A.vcf.snps.txt

（2）第二步利用上面的SNP去提取B.vcf，再次提取共有SNP文件

vcftools  --vcf B.vcf --snps \./A.vcf.snps.txt  --recode -c | \grep -v "#"| cut -f 3 > A.overlap.B.snps.txt

（3）第三步使用第二步产生的SNP再次过滤A.vcf和B.vcf

vcftools --vcf A.vcf \--snps A.overlap.B.snps.txt \--recode --out A.overlap.keepvcftools --vcf B.vcf \--snps A.overlap.B.snps.txt \--recode --out B.overlap.keep

（4）第四步合并两个vcf文件
此脚本使用时，注意两个vcf的行数相同，列中Sample不同，要合并。
VCFtools据说也能做，我用的时候出现问题，所以使用下面的脚本：

#!/usr/bin/env python
import os,sys
inf1=sys.argv[1] # sequenced res
inf2=sys.argv[2] # cau hapmap res
ouf=sys.argv[3]  # resultswith open(inf1) as infile1, open(inf2) as infile2, open(ouf,"w") as oufile:line1=infile1.readline() ## get first line as default valueline2=infile2.readline()while line1 != "" and line2 != "": if line1.startswith("##fileformat"):oufile.write(line1)line1=infile1.readline() ## iterate next lineline2=infile2.readline()elif line1.startswith("#CHROM"):a=line1.strip()b=line2.strip().split()bb="\t".join(b[9:])oufile.write(a+"\t"+bb+"\n")line1=infile1.readline() ## iterate next lineline2=infile2.readline()else:    line1=infile1.readline() ## iterate next lineline2=infile2.readline()#breaka=line1.strip().split()b=line2.strip().split()if a != [] and b != [] and a[3]==b[3]:if a[4]==b[4]:res= a + b[9:]oufile.write("%s\n"%("\t".join(res)))else:continue ## skip the Alt multiple code situationsres= a[:4] + [a[4] + ";" + b[4]] + a[5:] + b[9:]oufile.write("%s\n"%("\t".join(res)))print("Well Done !!")

2. VCF文件的文件头（VCF Header）

通常VCF文件的Header中包含了VCF中注释的信息，为了VCFtools识别且方便合并，我们仅仅使用两行作为表头即可，这样既可以识别VCF文件类型，也保存了Sample的样本信息，足够后面的使用了。

##fileformat=VCFv4.1
#CHROM ...

以上，仅为个人记录。