compare genotypes for the same sample from different sequencing experiments (eg. control data versus 1000 Genomes data for the same HapMap samples)

VCFcompare.py

#!/usr/bin/python

import vcf
import sys

Usage="""
python VCFcompare.py ${1KG_VCF} ${CONTROL_VCF} ${CHR} ${OUTFILE}
"""

if len(sys.argv) < 5:
  print Usage
  sys.exit(1)

def get_gts(genotype):
    if '/' in genotype:
        return map(int, genotype.split("/"))
    elif '|' in genotype:
        return map(int, genotype.split("|"))
    else:
        return int(genotype), int(genotype)

def main():
    
    file_kg    = sys.argv[1]
    file_pd    = sys.argv[2]
    chromosome = sys.argv[3]
    mm_output  = open(sys.argv[4], "w")
    
    vcf_kg = vcf.Reader(filename=file_kg)
    vcf_pd = vcf.Reader(filename=file_pd)
    
    kg_indices = []
    pd_indices  = []
    kg_samples = dict([(vcf_kg.samples[x], x) for x in range(len(vcf_kg.samples))])
    for i in xrange(len(vcf_pd.samples)):
        if vcf_pd.samples[i] in kg_samples:
            kg_indices.append(kg_samples[vcf_pd.samples[i]])
            pd_indices.append(i)
            
    num_matches    = 0
    num_mismatches = 0
    
    for pd_record in vcf_pd.fetch(chromosome, start=0, end=1000000000):
        if not pd_record.is_snp:
            continue
        
        kg_record = vcf_kg.fetch(pd_record.CHROM, start=pd_record.POS)
        if kg_record is None or not kg_record.is_snp:
            continue
        if kg_record is not None:
            for i in xrange(len(kg_indices)):
                if kg_record.samples[kg_indices[i]]["GT"] is None:
                    continue
                if pd_record.samples[pd_indices[i]]["GT"] is None:
                    continue
                
                gt_a1, gt_a2 = map(lambda x: pd_record.alleles[x], get_gts(pd_record.samples[pd_indices[i]]["GT"]))
                gt_b1, gt_b2 = map(lambda x: kg_record.alleles[x], get_gts(kg_record.samples[kg_indices[i]]["GT"]))
        
                if gt_a1 == gt_b1 and gt_a2 == gt_b2:
                    num_matches += 1
                elif gt_a1 == gt_b2 and gt_a2 == gt_b1:
                    num_matches += 1
                else:
                    num_mismatches +=1
                    mm_output.write("%s\t%d\t%s\n"%(pd_record.CHROM, pd_record.POS, vcf_kg.samples[kg_indices[i]]))
                    
    print("CHR:%s\tNUM_MATCHES=%d\tNUM_MISMATCHES=%d"%(chromosome, num_matches, num_mismatches))
    mm_output.close()
    
    
    
    
if __name__ == "__main__":
    main()