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

记录—提取合并VCF文件

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

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

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

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

import os,sys

inf=sys.argv[1] # input file
ouf=sys.argv[2] # output file

with 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.keep

vcftools --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]  # results

with open(inf1) as infile1, open(inf2) as infile2, open(ouf,"w") as oufile:
    line1=infile1.readline() ## get first line as default value
    line2=infile2.readline()
    while line1 != "" and line2 != "": 
        if line1.startswith("##fileformat"):
            oufile.write(line1)
            line1=infile1.readline() ## iterate next line
            line2=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 line
            line2=infile2.readline()
        else:    
            line1=infile1.readline() ## iterate next line
            line2=infile2.readline()
            #break
            a=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 situations
                    res= 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 ...

以上，仅为个人记录。