vcf.Writer outputs mangled vcf header #68

Closed
cooketho opened this Issue Aug 16, 2012 · 1 comment

Projects

None yet

2 participants

@cooketho

When I try out the following code the order of the lines in the original vcf file header is scrambled. The biggest problem with this is that the ##fileformat=VCFv4.1 line gets moved someplace else, and so GATK doesn't automatically recognize the file as a vcf. There may be a workaround by passing an option to GATK specifying the format, but really this is a problem with pyvcf, not GATK. Please update the next version so that the organization of the header is kept intact in the output vcf.

Python code:
"
import vcf
reader=vcf.Reader(filename='test.vcf')
writer=vcf.Writer(open('test.out.vcf', 'w'), reader)
for record in reader:
writer.write_record(record)
"

input vcf:
"

fileformat=VCFv4.1

INFO=<ID=LDAF,Number=1,Type=Float,Description="MLE Allele Frequency Accounting for LD">

INFO=<ID=AVGPOST,Number=1,Type=Float,Description="Average posterior probability from MaCH/Thunder">

INFO=<ID=RSQ,Number=1,Type=Float,Description="Genotype imputation quality from MaCH/Thunder">

INFO=<ID=ERATE,Number=1,Type=Float,Description="Per-marker Mutation rate from MaCH/Thunder">

INFO=<ID=THETA,Number=1,Type=Float,Description="Per-marker Transition rate from MaCH/Thunder">

INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">

INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">

INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">

INFO=<ID=HOMLEN,Number=.,Type=Integer,Description="Length of base pair identical micro-homology at event breakpoints">

INFO=<ID=HOMSEQ,Number=.,Type=String,Description="Sequence of base pair identical micro-homology at event breakpoints">

INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Difference in length between REF and ALT alleles">

INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">

INFO=<ID=AC,Number=.,Type=Integer,Description="Alternate Allele Count">

INFO=<ID=AN,Number=1,Type=Integer,Description="Total Allele Count">

ALT=<ID=DEL,Description="Deletion">

FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

FORMAT=<ID=DS,Number=1,Type=Float,Description="Genotype dosage from MaCH/Thunder">

FORMAT=<ID=GL,Number=.,Type=Float,Description="Genotype Likelihoods">

INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_alignments/README">

INFO=<ID=AF,Number=1,Type=Float,Description="Global Allele Frequency based on AC/AN">

INFO=<ID=AMR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from AMR based on AC/AN">

INFO=<ID=ASN_AF,Number=1,Type=Float,Description="Allele Frequency for samples from ASN based on AC/AN">

INFO=<ID=AFR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from AFR based on AC/AN">

INFO=<ID=EUR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from EUR based on AC/AN">

INFO=<ID=VT,Number=1,Type=String,Description="indicates what type of variant the line represents">

INFO=<ID=SNPSOURCE,Number=.,Type=String,Description="indicates if a snp was called when analysing the low coverage or exome alignment data">

source_20120714.1=vcf-subset(r731) -c NA18505,NA18508,NA19648,NA19704, downloads/ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA18505 NA18508 NA19648 NA19704

1 10583 rs58108140 G A 100 PASS AA=.;AC=0;AF=0.14;AFR_AF=0.04;AMR_AF=0.17;AN=8;ASN_AF=0.13;AVGPOST=0.7707;ERATE=0.0161;EUR_AF=0.21;LDAF=0.2327;RSQ=0.4319;SNPSOURCE=LOWCOV;THETA=0.0046;VT=SNP GT:DS:GL 0|0:0.000:-0.07,-0.80,-5.00 0|0:0.000:-0.01,-1.68,-5.00 0|0:0.100:-0.03,-1.15,-5.00 0|0:0.150:-0.19,-0.45,-2.42
"

output vcf:
"

ALT=<ID=DEL,Description="Deletion">

fileformat=VCFv4.1

source_20120714.1=vcf-subset(r731) -c NA18505,NA18508,NA19648,NA19704, downloads/ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz

INFO=<ID=ERATE,Number=1,Type=Float,Description="Per-marker Mutation rate from MaCH/Thunder">

INFO=<ID=AFR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from AFR based on AC/AN">

INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">

INFO=<ID=SNPSOURCE,Number=.,Type=String,Description="indicates if a snp was called when analysing the low coverage or exome alignment data">

INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Difference in length between REF and ALT alleles">

INFO=<ID=AVGPOST,Number=1,Type=Float,Description="Average posterior probability from MaCH/Thunder">

INFO=<ID=HOMSEQ,Number=.,Type=String,Description="Sequence of base pair identical micro-homology at event breakpoints">

INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_alignments/README">

INFO=<ID=AC,Number=.,Type=Integer,Description="Alternate Allele Count">

INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">

INFO=<ID=AF,Number=1,Type=Float,Description="Global Allele Frequency based on AC/AN">

INFO=<ID=ASN_AF,Number=1,Type=Float,Description="Allele Frequency for samples from ASN based on AC/AN">

INFO=<ID=AN,Number=1,Type=Integer,Description="Total Allele Count">

INFO=<ID=VT,Number=1,Type=String,Description="indicates what type of variant the line represents">

INFO=<ID=THETA,Number=1,Type=Float,Description="Per-marker Transition rate from MaCH/Thunder">

INFO=<ID=EUR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from EUR based on AC/AN">

INFO=<ID=LDAF,Number=1,Type=Float,Description="MLE Allele Frequency Accounting for LD">

INFO=<ID=HOMLEN,Number=.,Type=Integer,Description="Length of base pair identical micro-homology at event breakpoints">

INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">

INFO=<ID=AMR_AF,Number=1,Type=Float,Description="Allele Frequency for samples from AMR based on AC/AN">

INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">

INFO=<ID=RSQ,Number=1,Type=Float,Description="Genotype imputation quality from MaCH/Thunder">

FORMAT=<ID=GL,Number=.,Type=Float,Description="Genotype Likelihoods">

FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

FORMAT=<ID=DS,Number=1,Type=Float,Description="Genotype dosage from MaCH/Thunder">

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA18505 NA18508 NA19648 NA19704

1 10583 rs58108140 G A 100 . AA=.;AVGPOST=0.7707;AC=0;AF=0.14;ASN_AF=0.13;AFR_AF=0.04;AMR_AF=0.17;ERATE=0.0161;AN=8;LDAF=0.2327;VT=S;SNPSOURCE=LOWCOV;THETA=0.0046;RSQ=0.4319;EUR_AF=0.21 GT:DS:GL 0|0:0.000:-0.07,-0.8,-5.0 0|0:0.000:-0.01,-1.68,-5.0 0|0:0.100:-0.03,-1.15,-5.0 0|0:0.150:-0.19,-0.45,-2.42
"

@cooketho

Never mind. Sorry! I see that this was an issue with version 0.4. I installed version 0.6 and it seems to be working fine. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment