vcftools 012 apparently is giving me wrong genotypes

Question

[shinken123]

Hi,

I want to convert a vcf to a numeric format. I am using vctools 012 using vcftools/0.1.16.

This is how my vcf looks like:

zcat file.vcf.gz | grep  -v "#" | cut -f 1,2,10,11,12 | head -n 4
Chr01    1076    1|1     0|0     0|0
Chr01    1114    0|1     0|0     0|0
Chr01    1128    1|1     0|0     0|0 
Chr01    1202    0|0     0|0     0|0

What I would expect about this file in a numeric format is:

    2  1  2  0
    0  0  0  0
    0  0  0  0

With individuals in rows and SNPs in columns

However, using vcftools like this:

vcftools --gzvcf file.vcf.gz --012 --out file_num

I have

cut -f 1,2,3,4 file_num | head -n 3
0       2       1       2
1       0       0       0
2       0       0       0

I am wondering if I am understanding wrong the output from vcftools or if it has a bug. Many thanks-

Answer 1

[sbstevenlee]

If you are a Python user, below is a Python API solution using the pyvcf submodule from the fuc package I wrote:

Assume you have the following data:

>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr1', 'chr1', 'chr1'],
...     'POS': [1076, 1114, 1128, 1202],
...     'ID': ['.', '.', '.', '.'],
...     'REF': ['G', 'T', 'A', 'C'],
...     'ALT': ['A', 'C', 'T', 'G'],
...     'QUAL': ['.', '.', '.', '.'],
...     'FILTER': ['.', '.', '.', '.'],
...     'INFO': ['.', '.', '.', '.'],
...     'FORMAT': ['GT', 'GT', 'GT', 'GT'],
...     'A': ['1|1', '0|1', '1|1', '0|0'],
...     'B': ['0|0', '0|0', '0|0', '0|0'],
...     'C': ['0|0', '0|0', '0|0', '0|0'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> vf.df
  CHROM   POS ID REF ALT QUAL FILTER INFO FORMAT    A    B    C
0  chr1  1076  .   G   A    .      .    .     GT  1|1  0|0  0|0
1  chr1  1114  .   T   C    .      .    .     GT  0|1  0|0  0|0
2  chr1  1128  .   A   T    .      .    .     GT  1|1  0|0  0|0
3  chr1  1202  .   C   G    .      .    .     GT  0|0  0|0  0|0

We define a custom method (one_row) to be applied to each row:

>>> def one_row(r):
...     
...     def one_gt(g):
...         return sum([0 if x == '0' else 1 for x in g.split('|')])
...         
...     return r[9:].apply(one_gt)
...

We apply the custom method and transpose the resulting dataframe:

>>> df = vf.df.apply(one_row, axis=1).T
>>> df
   0  1  2  3
A  2  1  2  0
B  0  0  0  0
C  0  0  0  0

You can easily read and write VCF files:

vf = pyvcf.VcfFrame.from_file('in.vcf')
vf.to_file('out.vcf')