Editing GSCAN dbGaP (section)

= UK Biobank =
More information about the files used for [[UK Biobank|UKBiobank are here]]. In brief, we used the UK10K + 1kgp3 imputed vcfs provided by UKBionank and added in dosages w/ this python script:
<syntaxhighlight lang="python">
import gzip, argparse, re, os, datetime
from subprocess import Popen, PIPE

def add_dosage(pair):
        a, b = pair
        probs = b.split(b',')
        dose = float(probs[1]) + (float(probs[2]) * 2)
        return a + b':' + str(dose).encode('ascii') + b':' + b

def gziplines(fname):
  f = Popen(['zcat', fname], stdout=PIPE)
  for line in f.stdout:
      yield line

parser = argparse.ArgumentParser()
parser.add_argument('inputVCF', help = 'The path to the VCF')
args = parser.parse_args()

flag = False

for line in gziplines(args.inputVCF):
        if line.startswith(b'#'):
                os.write(1, line.rstrip() + b'\n')
                if not flag:
                        os.write(1, b'##FORMAT=<ID=DS,Number=1,Type=Float,Description="Genotype Dosages">\n')
                        os.write(1, b'##Dosages added using the script add.dosages.subprocess.py at ' +
                                str(datetime.datetime.now()).encode('ascii') + b'\n')
                        flag = True
        else:
                elements = re.split(b'\t|:', line.rstrip())
                first8 = elements[:8]
                genotypes = elements[10:]
                form = b'GT:DS:GP'

                genotypes_split = zip(genotypes[::2], genotypes[1::2])
                try:
                        dose_genos = [add_dosage(pair) for pair in genotypes_split]
                except (ValueError, IndexError) as e:
                        os.write(2, "\n" + line)
                        os.write(2, line + "\n" + args.inputVCF + "\n\n")
                        raise e
                os.write(1, b'\t'.join(first8) + b'\t' + form + b'\t' + b'\t'.join(dose_genos) + b'\n')
</syntaxhighlight>