1
$\begingroup$

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.

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)

$\endgroup$
0

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Perfect! It worked super good, I am so grateful to both of you! So one thing has gotten forgotten, one also needs to adjust the header (it's easy). Check my post here: biostars.org/p/9463791/#9467359 $\endgroup$ Apr 28, 2021 at 8:51
  • $\begingroup$ You're very welcome, @DavidStreuli. If one of the answers here solved your issue, please take a moment and accept it by clicking on the checkmark on the left. That is the best way to express your thanks on the Stack Exchange sites. Feel very fre to wait for a better answer, no need to accept mine now, I just wanted to let you know how the site works. $\endgroup$
    – terdon
    Apr 28, 2021 at 14:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.