vcftools 012 apparently is giving me wrong genotypes
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-
• 2,039 views
•
link
1 answer
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')
• 0 views
•
link
Log in to answer this question.