2
$\begingroup$

I'm working with gvcf files, containing also non-variant positions. I have noticed some insertions that I'm not able to correctly interpret or handle. An example (selected samples from a huge vcf):

NC_013991.2 1403    .   G   GAA 176,32  .   AC=4;AF=0.014;AN=286;DP=1294;ExcessHet=0.0001;FS=0;InbreedingCoeff=0.5011;MLEAC=3;MLEAF=0.01;MQ=28.09;QD=28.2;SOR=3.912 GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:1403_G_GAA:90,6,0:1403  0/0:2,0:2:6:.:.:0,6,88:.    ./.:0,0:1:.:.:.:0,0,0:. ./.:0,0:0:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 0/0:2,0:2:3:.:.:0,3,45:.    ./.:1,0:1:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 1/1:0,3:3:9:.:.:123,9,0:.   0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:2,0:2:6:.:.:0,6,90:.    0/0:9,0:9:21:.:.:0,21,315:.
NC_013991.2 1404    .   A   .   .   .   DP=1286 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21
NC_013991.2 1405    .   A   .   .   .   DP=1282 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21

Here, the reference sequence from position 1403 to 1405 is "GAA". At site 1403, an alternative allele "GAA" is indicated, which I would interpret as an insertion of "AA". Two individuals thus would have genotype "GAAAA" here while the rest has "GAA".

I've seen multiple cases like this, here with longer alternative allele (again, selected samples, not necessarily the same as before):

NC_013991.2 23  .   T   TAATACTTATGTGTTG    66,98   .   AC=2;AF=0.009346;AN=214;DP=314;ExcessHet=0.0105;FS=0;InbreedingCoeff=0.3797;MLEAC=3;MLEAF=0.014;MQ=51.77;QD=33.49;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:85,6,0  0/0:3,0:3:9:0,9,135 0/0:3,0:3:9:0,9,95  0/0:1,0:1:3:0,3,45  ./.:1,1:2:.:0,0,0   0/0:3,1:4:0:0,0,102
NC_013991.2 24  .   A   .   .   .   DP=322  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:2:2:6   0/0:4:4:12
NC_013991.2 25  .   A   .   .   .   DP=339  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:4:4:12
NC_013991.2 26  .   T   .   .   .   DP=353  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 27  .   A   .   .   .   DP=361  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 28  .   C   .   .   .   DP=368  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 29  .   T   .   .   .   DP=386  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 30  .   T   .   .   .   DP=388  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 31  .   A   .   .   .   DP=398  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 32  .   T   .   .   .   DP=405  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 33  .   G   .   .   .   DP=419  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 34  .   T   .   .   .   DP=427  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 35  .   G   .   0,09    LowQual BaseQRankSum=-1;DP=433;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=51.42;MQRankSum=1;ReadPosRankSum=0    GT:DP:RGQ   0/0:3:9 0/0:3:9 0/0:3:9 0/0:1:3 0/0:3:9 0/0:5:0
NC_013991.2 36  .   T   .   .   .   DP=441  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 37  .   T   .   .   .   DP=447  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   0/0:5:5:12
NC_013991.2 38  .   G   .   .   .   DP=449  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   ./.:4:5:0

The alternative allele "TAATACTTATGTGTTG" is exactly the reference. Of course, the sample in question is homozygous. It seems to me these are "false" insertions: they look like insertions, but aren't. I can write a simple awk script trimming the the alternative allele, but there are cases where there are real insertions, and I don't know how to distinguish those.

I've tried the solution mentioned in this answer (double pass through bcftools norm) but it doesn't resolve these issues.

I could spend a few days coding a C++ program to do this, but surely someone with better understanding of vcf, C++, etc, has already done this?

Here's a case I really don't know how to interpret:

NC_013991.2 3321    .   A   ATTTCATT    638,04  .   AC=14;AF=0.044;AN=320;DP=1873;ExcessHet=0;FS=0;InbreedingCoeff=0.5521;MLEAC=10;MLEAF=0.031;MQ=29.5;QD=34.42;SOR=1.739   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:2,0:2:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:33,0:33:96:.:.:0,96,1440:.  0/0:15,0:15:39:.:.:0,39,585:.
NC_013991.2 3322    .   T   .   .   .   DP=1870 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:2:2:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:33:33:96    0/0:14:14:36
NC_013991.2 3323    .   T   .   .   .   DP=1859 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   ./.:1:2:0   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:36
NC_013991.2 3324    .   T   .   .   .   DP=1849 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:33
NC_013991.2 3325    .   TGG T   639,64  .   AC=14;AF=0.044;AN=318;DP=1843;ExcessHet=0;FS=0;InbreedingCoeff=0.5416;MLEAC=11;MLEAF=0.035;MQ=29.53;QD=28.53;SOR=1.739  GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3326    .   G   T   1294,06 .   AC=15;AF=0.055;AN=274;BaseQRankSum=0;DP=1833;ExcessHet=0;FS=2.451;InbreedingCoeff=0.5364;MLEAC=17;MLEAF=0.062;MQ=56.66;MQRankSum=0;QD=16.38;ReadPosRankSum=0;SOR=0.828  GT:AD:DP:GQ:PGT:PID:PL:PS   ./.:0,2:2:.:.:.:0,0,0:. ./.:0,4:4:.:.:.:0,0,0:. ./.:0,2:2:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 1|1:0,1:1:3:1|1:3309_G_A:45,3,0:3309    1/1:0,4:4:12:.:.:145,12,0:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3327    .   G   .   0,02    LowQual DP=1832;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=29   GT:DP:RGQ   0/0:2:0 0/0:4:0 0/0:2:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:3 0/0:8:18    0/0:34:99   0/0:14:33
NC_013991.2 3328    .   G   A   650,66  .   AC=15;AF=0.048;AN=312;BaseQRankSum=0;DP=1826;ExcessHet=0;FS=3.31;InbreedingCoeff=0.5371;MLEAC=12;MLEAF=0.038;MQ=34.34;MQRankSum=0;QD=27.23;ReadPosRankSum=0;SOR=0.324   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:18:.:.:0,18,270:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3329    .   G   .   18,49   LowQual BaseQRankSum=-1;DP=1818;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=57.63;MQRankSum=0;ReadPosRankSum=1   GT:DP:RGQ   0/0:2:6 0/0:4:12    0/0:2:6 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:7:18    0/0:35:99   0/0:13:33

It seems like only the "CATT" part of the alternative allele is really inserted...

The mapping was done either with bwa aln or bwa mem, and we've observed this behaviour both when using GATK4 and bcftools for genotyping.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ "It seems to me these are "false" insertions: they look like insertions, but aren't.": could you explain why they seem false to you? And if we don't need to know about all of the genotypes, please only show us the relevant ones. You say "of course the sample in question is homozygous", but we can't see that since there are so many samples shown. Can you focus only on the relevant ones so we can understand better? $\endgroup$
    – terdon ♦
    Dec 23, 2022 at 11:57
  • $\begingroup$ They look false to me because they exactly recapitulate the reference sequence. I removed samples in the examples and could add a few lines illustrating a longer pseudo-insertion. $\endgroup$
    – Jos
    Dec 23, 2022 at 15:18
  • 1
    $\begingroup$ Sorry if I'm being dense, this sort of thing is always complicated to understand, despite working with this sort of thing every day. Focusing on your first example, the VCF with three lines, we have an insertion of AA, and then two non-variant positions. So some of your samples have GAAAA and some have GAA, the reference, as you say. What am I missing? What's wrong with that scenario? Isn't this just a case of a duplication indel? $\endgroup$
    – terdon ♦
    Dec 23, 2022 at 17:52
  • $\begingroup$ Thanks for sharing your thoughts. I started with the first example as it is short, but most of these "insertions" are rather long such as the second one. I found it suspicious that the large majority of insertions /exactly/ recapitulates the reference sequence. But maybe that's common. Also, we don't seem to get these when mapping RNAseq data with STAR, but with BWA, we have them also in the coding sequences. If the second example really means there are individuals with TAATACTTATGTGTTG and others with TAATACTTATGTGTTGAATACTTATGTGTTG, I invite you to post an answer that I'll accept. $\endgroup$
    – Jos
    Dec 23, 2022 at 21:00

2 Answers 2

Reset to default
2
$\begingroup$

Microsatellites, minisatellites and other type of repetitive variation are a known feature of genomes.

I agree with previous comments in the sense that what you are observing here might as well be real and not just an artefact caused by real repetitive variation in your genome of study, plus read mapping together with variant calling issues at the site.

The latter might be the case, but before getting rid of these alleles I suggest that you inspect the sites and regions around these alleles more carefully: GATK HaplotypeCaller has an option (-- output-bam or similar) that writes a BAM file with the reads as they were used for variant calling after the tool itself performs local re-alignment of reads around the variants. I would use this tool together with the original BAM file in order to inspect the region visually and try to understand it better. For this you can use samtools tview or IGV. The length of your reads becomes somewhat relevant here, but even if you are sequencing around 100 bp, given the small size of these observed alleles you should be able to draw some more solid insights about the variants/errors and how to handle them.

Improve this answer
$\endgroup$
0
$\begingroup$

INDELs tend to occur tandemly. Your observation makes them more likely to be real.

Improve this answer
$\endgroup$

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.