Why is there no GT-entry in my .vcf file from bacteria

Question

I mapped raw reads obtained from an E. coli genome (.fastq file) to my reference genome. Next, I exported the .vcf (variant calling file), which looks something like the header example shown below (total >50 lines):

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   DP:AO   231:89

How can I back-fill or quickly assign the missing GTs to the .vcf ?

For my downstream process (plink --pca, via bfctools merge), I need that there are GT entries in every .vcf, so I can merge them and the respective GT is shown in the merged .vcf. I only found solutions that need the .bam file, but that is laborious and I don't think I need the .bam in the case of haploid data.

• This method needs .bam: https://www.biostars.org/p/319304/
• This method needs .fastq: https://training.galaxyproject.org/training-material/topics/variant-analysis/tutorials/non-dip/tutorial.html#examining-the-results

However, the.fastq is actually from a mixed culture, so I wonder do I really need to follow the .bam-strategy? Suggestion for a tool? But in that case, why did it not export the GTs in the fist place (I exported from "geneious" sequence analysis program)

Answer 1

The representation for a variant present in a haploid genome should be 1. A heterozygote in a diploid genome would be 0/1 (or 0|1 or 1|0 if phased). So, to add a 1 as the genotype, you would need to first insert a GT to the beginning of each INFO line and then the 1 at the beginning of the last field. But, this should happen only on lines that are not headers (do not start with #). This should do it:

$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf 
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   GT:DP:AO    1:231:89

So just redirect the output to a fixed file:

$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf > fixed.vcf

Now, in this specific case, you could also do something like this:

$ sed -i.bak '/^[^#]/{s/\<DP:AO\t/GT:DP:AO\t1:/}' file.vcf 

But that assumes that all lines will have the exact same INFO field of DP:AO which is not a safe assumption since each variant can have different INFO fields.