4
$\begingroup$

Let's say I have a pandas dataframe with fields CHROM, POS, ALT, REF. In this special case, I also wouldn't care about IDor FILTER, INFO could be blank (or meaningless) and we'll write QUAL as each 40.

import pandas as pd
example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)
print(df)

   CHROM      POS   REF ALT
0     20    14370     G   A
1     20    17330     T   A
2     20  1110696     A   G
3     20  1230237     T   C
4     20  1234567  GTCT   G

In this case, saving the output with the extension .vcf via df.to_csv("myfile.vcf") with the other required fields would created the tab-delimited format of the VCF, but it wouldn't include the mandatory header. https://samtools.github.io/hts-specs/VCFv4.1.pdf

##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
...
#CHROM POS ID REF ALT QUAL FILTER INFO

The most popular package to work with VCF's appears to be PyVCF: https://pypi.python.org/pypi/PyVCF

However, one needs to begin with a VCF as a template, and add to this (from the documentation):

vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
for record in vcf_reader:
     vcf_writer.write_record(record)

It appears the only information necessary from vcf.Reader(filename='vcf/test/tb.vcf.gz') is the metadata?

print(vcf_reader.metadata)
OrderedDict([('fileformat', 'VCFv4.0'), ('fileDate', '20090805'), ('source', ['myImputationProgramV3.1']), ('reference', '1000GenomesPilot-NCBI36'), ('phasing', ['partial'])])

Otherwise, tb.vcf.gz already has variants, which I think means I would need to feed in a VCF with no variants?:

for record in vcf_reader:
     print(record)

Record(CHROM=20, POS=14370, REF=G, ALT=[A])
Record(CHROM=20, POS=17330, REF=T, ALT=[A])
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Record(CHROM=20, POS=1234567, REF=GTCT, ALT=[G, GTACT])
$\endgroup$

1 Answer 1

8
$\begingroup$

You will have to define the header either from an existing VCF file or hardcoded into you python script. Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.

You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ID, FILTER, INFO and QUAL columns

import pandas as pd 

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']] 

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

Then myfile.vcf looks like

##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40      
20  17330       T   A   40      
20  1110696     A   G   40      
20  1230237     T   C   40      
20  1234567     GTCT    G   40
$\endgroup$
2
  • $\begingroup$ Thanks for the help. BTW, is fileDate year-month-day or year-day-month? I'm guessing there's no standard. $\endgroup$ Dec 15, 2017 at 20:01
  • 1
    $\begingroup$ fileDate is not a required header field you can remove it it you want $\endgroup$
    – Bioathlete
    Dec 15, 2017 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.