This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Convert the genotype format from numeric to letters

I want to use IVAS for sQTL analysis and it accepts only allelic encoding of genotypes, so that they should be two letters of A,C,G,T

The format of my vcf file is like this:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  108 139 159 265

1   73  .   C   A   40  PASS    .   GT:DP:GQ    0|0:5:40    0|0:9:40    0|0:6:38    ./.:.:.

1   83  .   T   C,A 40  PASS    .   GT:DP:GQ    1|1:5:40    1|1:9:40    0|0:8:38    ./.:.:.

I want to convert the genotype format from numeric to letters

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  108 139 159 265

1   73  .   C   A   40  PASS    .   GT:DP:GQ    CC  CC  CA  NA

1   83  .   T   C,A 40  PASS    .   GT:DP:GQ    TC  TC  TT  NA
snp vcf
• 4,141 views • link
updated by shaunjc • shaunjc • written by waqaskhokhar999 waqaskhokhar999

2 answers 

#!/usr/bin/env python2
#!/coding:utf-8
import sys
import os
import gzip


class Letter(object):
    standard_nt = ['A', 'C', 'G', 'T']

    def __init__(self, vcf):
        self.vcf = vcf
        self.get_head()

    def _open(self):
        if os.path.splitext(self.vcf)[-1] == '.gz':
            input = gzip.open(self.vcf)
        else:
            input = open(self.vcf)
        return input

    def get_head(self):
        input = self._open()
        self.vcfhead = []
        for h in input:
            if h.startswith('#'):
                self.vcfhead.append(h)
            else:
                break
        input.close()

    def parse_vcf(self):
        input = self._open()
        for i in input:
            if i.startswith('#'):
                continue
            else:
                line = i.strip().split()
                line[7]='.'
                line[8]='GT'
                ref = line[3]
                alt = line[4]
                if len(ref) > 1 or len(alt) > 1:
                    print 'only accept BIALLELIC SNP\nremove this site\n'
                    continue
                elif ref not in self.standard_nt or alt not in self.standard_nt:
                    print 'only accept valid nucleotide\n'
                    continue
                else:
                    part1 = '\t'.join(line[0:9])
                    sample_gt = []
                    for gt in line[9:]:
                        gt = gt.split(':')[0]
                        if gt == '0|0' or gt == '0/0':
                            letter = ref*2
                        elif gt == '0|1'or gt == '0/1':
                            letter = ref+alt
                        elif gt == '1|1'or gt == '1/1':
                            letter = alt*2
                        else:
                            letter = 'NA'
                        sample_gt.append(letter)
                    sample_gt = '\t'.join(sample_gt)
                    all = part1+'\t'+sample_gt
                    yield all
        input.close()

    def write_out(self):
        with gzip.open('letter.'+self.vcf, 'w') as f:
            f.write(''.join(self.vcfhead))
            for line in self.parse_vcf():
                f.write(line+'\n')

if __name__ == '__main__':
    vcf = sys.argv[1]
    Letter(vcf).write_out()

save the script as letter.py,run `python letter.py yours.vcf(.gz format is ok too).The file with 'letter.' prefix is the result file.

• 0 views • link
written by Leon Leon

Many thanks, the genotype format of my vcf file is like this 0|0, 1|1, 2|2, ./. so I replaced

if gt == '0|0' or gt == '0/0':
                        letter = ref*2
                    elif gt == '0|1'or gt == '0/1':
                        letter = ref+alt
                    elif gt == '1|1'or gt == '1/1':
                        letter = alt*2

by 

if gt == '0|0' or gt == '0/0':
                        letter = ref*2
                    elif gt == '1|1'or gt == '1/1':
                        letter = ref+alt
                    elif gt == '2|2'or gt == '2/2':
                        letter = alt*2

It needs parenthesis in print statement so I added parenthesis

(print 'only accept BIALLELIC SNP\nremove this site\n')

But I am getting this error:

Traceback (most recent call last):
  File "letter.py", line 76, in <module>
    Letter(vcf).write_out()
  File "letter.py", line 70, in write_out
    f.write(''.join(self.vcfhead))
  File "/home/waqas/miniconda3/lib/python3.6/gzip.py", line 260, in write
    data = memoryview(data)
TypeError: memoryview: a bytes-like object is required, not 'str'
• 0 views • link
written by waqaskhokhar999 waqaskhokhar999

Actually, it is running in py2 not py3. By the way, it is a little weird. Is standard format of your vcf ? If so, I hope my code will resolve your problem without change. But if you just want to use 0 represent refHOM ,1 represent HET, 2 represent altHOM, i think you will meet some problems in data format when use the script.

• 0 views • link
written by Leon Leon

Many thanks for your response, I have downloaded vcf file from here. |In vcf file the genotype fromats are like this:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  108 139 159 265

1   73  .   C   A   40  PASS    .   GT:DP:GQ    0|0:5:40    0|0:9:40    0|0:6:38    ./.:.:.

1   83  .   T   C,A 40  PASS    .   GT:DP:GQ    1|1:5:40    1|1:9:40    0|0:8:38    ./.:.:.

I was thinking if your code tries to find 1|0 or 0|1 and tries to replace it with het as I tried to grep 0|1 and 1|0 from vcf file and did not get any count.

• 0 views • link
written by waqaskhokhar999 waqaskhokhar999

I have a similar problem, where I have downloaded a file that seems to be in phased nucleotide format (A|A, A|G, G|G, etc). How could I modify this script to convert it back based into the ref and alt of to have phased numeric values (0|0, 0|1, 1|1)?

I feel like it would need a lot of modifying to determine which of the nucleotide matches ref or alt and then convert accordingly?

• 0 views • link
written by shaunjc • shaunjc

Log in to answer this question.

More posts like this