import pysam
from numpy.random import choice
import math


def assign_random_gt(input_vcf, outname, sample_name="HG", default_af=0.01):
    vcf_pointer = pysam.VariantFile(filename=input_vcf)
    new_header = vcf_pointer.header.copy()
    if "GT" not in new_header.formats:
        new_header.formats.add("GT", "1", "String", "Consensus Genotype across all datasets with called genotype")
        new_header.samples.add(sample_name)

    default_probs = [1 - default_af - math.pow(default_af, 2), default_af, math.pow(default_af, 2)]
    with open(outname, 'w') as out_vcf:
        out_vcf.write(str(new_header))
        for rec in vcf_pointer.fetch():
            rec_copy = rec.copy()
            if "GT" not in rec_copy.format.keys():
                if "AF" not in rec_copy.info.keys():
                    gt_probs = default_probs
                else:
                    af = rec_copy.info["AF"]
                    gt_probs = [1 - af - math.pow(af, 2), af, math.pow(af, 2)]
                c = choice(["0/0", "0/1", "1/1"], p=gt_probs)
                out_vcf.write("\t".join([str(rec_copy)[:-1], "GT", c]) + "\n")
    vcf_pointer.close()

import pysam
from numpy.random import choice


def assign_random_gt(input_vcf, outname, sample_name="HG", default_af=0.01):
    vcf_pointer = pysam.VariantFile(filename=input_vcf)
    new_header = vcf_pointer.header.copy()
    if "GT" not in new_header.formats:
        new_header.formats.add("GT", "1", "String", "Consensus Genotype across all datasets with called genotype")
        new_header.samples.add(sample_name)

    default_probs = [1 - default_af * (1 + default_af), default_af, default_af * default_af]
    with open(outname, 'w') as out_vcf:
        out_vcf.write(str(new_header))
        for rec in vcf_pointer.fetch():
            rec_copy = rec.copy()
            if "GT" not in rec_copy.format.keys():
                if "AF" not in rec_copy.info.keys():
                    gt_probs = default_probs
                else:
                    af = rec_copy.info["AF"]
                    gt_probs = [1 - af * (1 + af), af, af * af]
                c = choice(["0/0", "0/1", "1/1"], p=gt_probs)
                out_vcf.write("\t".join([str(rec_copy)[:-1], "GT", c]) + "\n")
    vcf_pointer.close()

Add random GT to a given vcf
Add random GT to a given vcf

  import pysam
  from numpy.random import choice
- import math
  
  
  def assign_random_gt(input_vcf, outname, sample_name="HG", default_af=0.01):
      vcf_pointer = pysam.VariantFile(filename=input_vcf)
      new_header = vcf_pointer.header.copy()
      if "GT" not in new_header.formats:
          new_header.formats.add("GT", "1", "String", "Consensus Genotype across all datasets with called genotype")
          new_header.samples.add(sample_name)
  
-     default_probs = [1 - default_af - math.pow(default_af, 2), default_af, math.pow(default_af, 2)]
+     default_probs = [1 - default_af * (1 + default_af), default_af, default_af * default_af]
      with open(outname, 'w') as out_vcf:
          out_vcf.write(str(new_header))
          for rec in vcf_pointer.fetch():
              rec_copy = rec.copy()
              if "GT" not in rec_copy.format.keys():
                  if "AF" not in rec_copy.info.keys():
                      gt_probs = default_probs
                  else:
                      af = rec_copy.info["AF"]
-                     gt_probs = [1 - af - math.pow(af, 2), af, math.pow(af, 2)]
+                     gt_probs = [1 - af * (1 + af), af, af * af]
                  c = choice(["0/0", "0/1", "1/1"], p=gt_probs)
                  out_vcf.write("\t".join([str(rec_copy)[:-1], "GT", c]) + "\n")
      vcf_pointer.close()
  
Add random GT to a given vcf
## Code Before:
import pysam
from numpy.random import choice
import math


def assign_random_gt(input_vcf, outname, sample_name="HG", default_af=0.01):
    vcf_pointer = pysam.VariantFile(filename=input_vcf)
    new_header = vcf_pointer.header.copy()
    if "GT" not in new_header.formats:
        new_header.formats.add("GT", "1", "String", "Consensus Genotype across all datasets with called genotype")
        new_header.samples.add(sample_name)

    default_probs = [1 - default_af - math.pow(default_af, 2), default_af, math.pow(default_af, 2)]
    with open(outname, 'w') as out_vcf:
        out_vcf.write(str(new_header))
        for rec in vcf_pointer.fetch():
            rec_copy = rec.copy()
            if "GT" not in rec_copy.format.keys():
                if "AF" not in rec_copy.info.keys():
                    gt_probs = default_probs
                else:
                    af = rec_copy.info["AF"]
                    gt_probs = [1 - af - math.pow(af, 2), af, math.pow(af, 2)]
                c = choice(["0/0", "0/1", "1/1"], p=gt_probs)
                out_vcf.write("\t".join([str(rec_copy)[:-1], "GT", c]) + "\n")
    vcf_pointer.close()

## Instruction:
Add random GT to a given vcf
## Code After:
import pysam
from numpy.random import choice


def assign_random_gt(input_vcf, outname, sample_name="HG", default_af=0.01):
    vcf_pointer = pysam.VariantFile(filename=input_vcf)
    new_header = vcf_pointer.header.copy()
    if "GT" not in new_header.formats:
        new_header.formats.add("GT", "1", "String", "Consensus Genotype across all datasets with called genotype")
        new_header.samples.add(sample_name)

    default_probs = [1 - default_af * (1 + default_af), default_af, default_af * default_af]
    with open(outname, 'w') as out_vcf:
        out_vcf.write(str(new_header))
        for rec in vcf_pointer.fetch():
            rec_copy = rec.copy()
            if "GT" not in rec_copy.format.keys():
                if "AF" not in rec_copy.info.keys():
                    gt_probs = default_probs
                else:
                    af = rec_copy.info["AF"]
                    gt_probs = [1 - af * (1 + af), af, af * af]
                c = choice(["0/0", "0/1", "1/1"], p=gt_probs)
                out_vcf.write("\t".join([str(rec_copy)[:-1], "GT", c]) + "\n")
    vcf_pointer.close()
