search

connor-lab / ncov2019-artic-nf

A Nextflow pipeline for running the ARTIC network's fieldbioinformatics tools (https://github.com/artic-network/fieldbioinformatics), with a focus on ncov2019
GNU Affero General Public License v3.0
88 stars 89 forks source link

Frameshift variants may interfere with typing in type_vcf.py, even with low ALT_FREQ #92

Open dfornika opened 3 years ago

dfornika commented 3 years ago

Running type_vcf.py -dp 10 -af 0.75 on the following input:

REGION  POS REF ALT REF_DP  REF_RV  REF_QUAL    ALT_DP  ALT_RV  ALT_QUAL    ALT_FREQ    TOTAL_DP    PVAL    PASS    GFF_FEATURE REF_CODON   REF_AA  ALT_CODON   ALT_AA
MN908947.3  241 C   T   0   0   0   653 118 44  1   653 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  913 C   T   1   0   39  590 289 42  0.998308    591 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  3037    C   T   0   0   0   247 162 43  1   247 5.21659e-158    TRUE    NA  NA  NA  NA  NA
MN908947.3  3267    C   T   0   0   0   341 2   42  1   341 2.44084e-216    TRUE    NA  NA  NA  NA  NA
MN908947.3  5388    C   A   0   0   0   153 2   40  0.993506    154 5.24955e-95 TRUE    NA  NA  NA  NA  NA
MN908947.3  5986    C   T   0   0   0   155 70  40  1   155 6.34799e-96 TRUE    NA  NA  NA  NA  NA
MN908947.3  6954    T   C   0   0   0   250 167 40  1   250 3.29933e-154    TRUE    NA  NA  NA  NA  NA
MN908947.3  11287   G   -TCTGGTTTT  618 522 41  655 0   20  0.967504    677 1.10773e-145    TRUE    NA  NA  NA  NA  NA
MN908947.3  14408   C   T   0   0   0   260 163 43  1   260 2.14099e-166    TRUE    NA  NA  NA  NA  NA
MN908947.3  14676   C   T   0   0   0   110 1   42  1   110 7.17302e-69 TRUE    NA  NA  NA  NA  NA
MN908947.3  15279   C   T   0   0   0   127 79  41  1   127 2.25986e-79 TRUE    NA  NA  NA  NA  NA
MN908947.3  16176   T   C   0   0   0   217 110 40  1   217 6.3057e-135 TRUE    NA  NA  NA  NA  NA
MN908947.3  17066   A   T   329 163 43  195 99  43  0.372137    524 7.79486e-76 TRUE    NA  NA  NA  NA  NA
MN908947.3  17609   A   -T  346 292 44  135 0   20  0.335821    402 1.56485e-34 TRUE    NA  NA  NA  NA  NA
MN908947.3  17615   A   G   0   0   0   337 283 44  1   337 2.76713e-218    TRUE    NA  NA  NA  NA  NA
MN908947.3  21701   G   +T  21  1   51  10  0   20  0.357143    28  0.000849992 TRUE    NA  NA  NA  NA  NA
MN908947.3  21764   A   -TACATG 18  6   47  22  0   20  1   22  5.12135e-06 TRUE    NA  NA  NA  NA  NA
MN908947.3  21834   A   -T  23  8   39  7   0   20  0.291667    24  0.0132072   FALSE   NA  NA  NA  NA  NA
MN908947.3  21990   T   -TTA    19  8   41  20  0   20  0.909091    22  1.64422e-05 TRUE    NA  NA  NA  NA  NA
MN908947.3  22000   C   +A  21  10  39  9   0   20  0.409091    22  0.00470272  TRUE    NA  NA  NA  NA  NA
MN908947.3  23063   A   T   0   0   0   235 104 39  1   235 9.76466e-144    TRUE    NA  NA  NA  NA  NA
MN908947.3  23271   C   A   0   0   0   209 132 41  1   209 2.13804e-130    TRUE    NA  NA  NA  NA  NA
MN908947.3  23403   A   G   0   0   0   199 137 44  1   199 1.41954e-129    TRUE    NA  NA  NA  NA  NA
MN908947.3  23604   C   A   1   0   77  92  26  39  0.989247    93  2.18401e-54 TRUE    NA  NA  NA  NA  NA
MN908947.3  23709   C   T   0   0   0   102 25  44  1   102 5.7399e-66  TRUE    NA  NA  NA  NA  NA
MN908947.3  24506   T   G   0   0   0   120 86  41  1   120 5.37115e-76 TRUE    NA  NA  NA  NA  NA
MN908947.3  24914   G   C   0   0   0   273 95  39  1   273 9.11315e-170    TRUE    NA  NA  NA  NA  NA
MN908947.3  27972   C   T   0   0   0   513 96  44  1   513 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28048   G   T   0   0   0   744 188 43  1   744 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28111   A   G   0   0   0   839 212 42  1   839 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28270   T   -A  731 344 41  686 0   20  0.873885    785 2.59338e-153    TRUE    NA  NA  NA  NA  NA
MN908947.3  28280   G   C   0   0   0   618 293 40  1   618 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28281   A   T   0   0   0   620 297 41  1   620 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28282   T   A   0   0   0   634 309 41  1   634 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28881   G   A   1   1   33  739 446 42  0.998649    740 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28882   G   A   1   1   33  743 452 41  0.998656    744 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28883   G   C   1   1   33  747 454 42  0.998663    748 0   TRUE    NA  NA  NA  NA  NA
MN908947.3  28977   C   T   0   0   0   501 351 40  1   501 9.95898e-309    TRUE    NA  NA  NA  NA  NA

...results in the following .csq.vcf:

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##source=iVar
##contig=<ID=MN908947.3>
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FILTER=<ID=FAIL,Description="Result of p-value > 0.05">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=REF_DP,Number=1,Type=Integer,Description="Depth of reference base">
##FORMAT=<ID=REF_RV,Number=1,Type=Integer,Description="Depth of reference base on reverse reads">
##FORMAT=<ID=REF_QUAL,Number=1,Type=Integer,Description="Mean quality of reference base">
##FORMAT=<ID=ALT_DP,Number=1,Type=Integer,Description="Depth of alternate base">
##FORMAT=<ID=ALT_RV,Number=1,Type=Integer,Description="Depth of alternate base on reverse reads">
##FORMAT=<ID=ALT_QUAL,Number=1,Type=String,Description="Mean quality of alternate base">
##FORMAT=<ID=ALT_FREQ,Number=1,Type=String,Description="Frequency of alternate base">
##bcftools/csqVersion=1.11+htslib-1.11
##bcftools/csqCommand=csq -v 0 -f MN908947.3.fasta -g MN908947.3.gff; Date=Fri Jan 15 15:08:49 2021
##INFO=<ID=BCSQ,Number=.,Type=String,Description="Haplotype-aware consequence annotation from BCFtools/csq, see http://samtools.github.io/bcftools/howtos/csq-calling.html for details. Format: Consequence|gene|transcript|biotype|strand|amino_acid_change|dna_change">
##FORMAT=<ID=BCSQ,Number=.,Type=Integer,Description="Bitmask of indexes to INFO/BCSQ, with interleaved first/second haplotype. Use \"bcftools query -f'[%CHROM\t%POS\t%SAMPLE\t%TBCSQ\n]'\" to translate.">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  SAMPLE
MN908947.3  241 .   C   T   .   PASS    DP=653  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ   1:0:0:0:653:118:44:1
MN908947.3  913 .   C   T   .   PASS    DP=591;BCSQ=synonymous|ORF1ab|ENSSAST00005000003|protein_coding|+|216S|913C>T,synonymous|ORF1ab|ENSSAST00005000002|protein_coding|+|216S|913C>T GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:0:39:590:289:42:0.998308:5
MN908947.3  3037    .   C   T   .   PASS    DP=247;BCSQ=synonymous|ORF1ab|ENSSAST00005000003|protein_coding|+|924F|3037C>T,synonymous|ORF1ab|ENSSAST00005000002|protein_coding|+|924F|3037C>T   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:247:162:43:1:5
MN908947.3  3267    .   C   T   .   PASS    DP=341;BCSQ=missense|ORF1ab|ENSSAST00005000003|protein_coding|+|1001T>1001I|3267C>T,missense|ORF1ab|ENSSAST00005000002|protein_coding|+|1001T>1001I|3267C>T GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:341:2:42:1:5
MN908947.3  5388    .   C   A   .   PASS    DP=154;BCSQ=missense|ORF1ab|ENSSAST00005000003|protein_coding|+|1708A>1708D|5388C>A,missense|ORF1ab|ENSSAST00005000002|protein_coding|+|1708A>1708D|5388C>A GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:153:2:40:0.993506:5
MN908947.3  5986    .   C   T   .   PASS    DP=155;BCSQ=synonymous|ORF1ab|ENSSAST00005000003|protein_coding|+|1907F|5986C>T,synonymous|ORF1ab|ENSSAST00005000002|protein_coding|+|1907F|5986C>T GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:155:70:40:1:5
MN908947.3  6954    .   T   C   .   PASS    DP=250;BCSQ=missense|ORF1ab|ENSSAST00005000003|protein_coding|+|2230I>2230T|6954T>C,missense|ORF1ab|ENSSAST00005000002|protein_coding|+|2230I>2230T|6954T>C GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:250:167:40:1:5
MN908947.3  11287   .   GTCTGGTTTT  G   .   PASS    DP=677;BCSQ=inframe_deletion|ORF1ab|ENSSAST00005000003|protein_coding|+|3674LSGF>3674L|11287GTCTGGTTTT>G,inframe_deletion|ORF1ab|ENSSAST00005000002|protein_coding|+|3674LSGF>3674L|11287GTCTGGTTTT>G   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:618:522:41:655:0:20:0.967504:5
MN908947.3  14408   .   C   T   .   PASS    DP=260;BCSQ=synonymous|ORF1ab|ENSSAST00005000002|protein_coding|+|4715L|14408C>T    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:260:163:43:1:1
MN908947.3  14676   .   C   T   .   PASS    DP=110;BCSQ=missense|ORF1ab|ENSSAST00005000002|protein_coding|+|4804P>4801L|14676C>T    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:110:1:42:1:1
MN908947.3  15279   .   C   T   .   PASS    DP=127;BCSQ=missense|ORF1ab|ENSSAST00005000002|protein_coding|+|5005T>5002I|15279C>T    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:127:79:41:1:1
MN908947.3  16176   .   T   C   .   PASS    DP=217;BCSQ=missense|ORF1ab|ENSSAST00005000002|protein_coding|+|5304L>5301P|16176T>C    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:217:110:40:1:1
MN908947.3  17066   .   A   T   .   PASS    DP=524;BCSQ=missense|ORF1ab|ENSSAST00005000002|protein_coding|+|5601I>5598F|17066A>T    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:329:163:43:195:99:43:0.372137:1
MN908947.3  17609   .   AT  A   .   PASS    DP=402;BCSQ=frameshift|ORF1ab|ENSSAST00005000002|protein_coding|+|5782IISLKHIKTNQLNALKCFIRVLSRMMFHLQLTGHK*>5779K*|17609AT>A+17615A>G    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:346:292:44:135:0:20:0.335821:1
MN908947.3  17615   .   A   G   .   PASS    DP=337;BCSQ=@17609  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:337:283:44:1:1
MN908947.3  21701   .   G   GT  .   PASS    DP=28;BCSQ=inframe_deletion|S|ENSSAST00005000004|protein_coding|+|47VLHSTQDLFLPFFSNVTWFHAIHVSGTNGTKRFDNPVLPFNDGVY>47VFTFNSGLVLTFLFQCYLVPCYLWDQWY*|21701G>GT+21764ATACATG>A+21834AT>A    GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:21:1:51:10:0:20:0.357143:1
MN908947.3  21764   .   ATACATG A   .   PASS    DP=22;BCSQ=@21701   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:18:6:47:22:0:20:1:1
MN908947.3  21834   .   AT  A   .   FAIL    DP=24;BCSQ=@21701   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:23:8:39:7:0:20:0.291667:1
MN908947.3  21990   .   TTTA    T   .   PASS    DP=22;BCSQ=*inframe_deletion|S|ENSSAST00005000004|protein_coding|+|143VY>141V|21990TTTA>T   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:19:8:41:20:0:20:0.909091:1
MN908947.3  22000   .   C   CA  .   PASS    DP=22;BCSQ=*frameshift|S|ENSSAST00005000004|protein_coding|+|146HKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLEGKQGNFKNLREFVFKNIDGYFKIYSKHTPINLVRDLPQGFSALEPLVDLPIGINITRFQTLLALHRSYLTPGDSSSGWTAGAAAYYVGYLQPRTFLLKYNENGTITDAVDCALDPLSETKCTLKSFTVEKGIYQTSNFRVQPTESIVRFPNITNLCPFGEVFNATRFASVYAWNRKRISNCVADYSVLYNSASFSTFKCYGVSPTKLNDLCFTNVYADSFVIRGDEVRQIAPGQTGKIADYNYKLPDDFTGCVIAWNSNNLDSKVGGNYNYLYRLFRKSNLKPFERDISTEIYQAGSTPCNGVEGFNCYFPLQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVCGPKKSTNLVKNKCVNFNFNGLTGTGVLTESNKKFLPFQQFGRDIADTTDAVRDPQTLEILDITPCSFGGVSVITPGTNTSNQVAVLYQDVNCTEVPVAIHADQLTPTWRVYSTGSNVFQTRAGCLIGAEHVNNSYECDIPIGAGICASYQTQTNSPRRARSVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSMTKTSVDCTMYICGDSTECSNLLLQYGSFCTQLNRALTGIAVEQDKNTQEVFAQVKQIYKTPPIKDFGGFNFSQILPDPSKPSKRSFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVLPPLLTDEMIAQYTSALLAGTITSGWTFGAGAALQIPFAMQMAYRFNGIGVTQNVLYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVEAEVQIDRLITGRLQSLQTYVTQQLIRAAEIRASANLAATKMSECVLGQSKRVDFCGKGYHLMSFPQSAPHGVVFLHVTYVPAQEKNFTTAPAICHDGKAHFPREGVFVSNGTHWFVTQRNFYEPQIITTDNTFVSGNCDVVIGIVNNTVYDPLQPELDSFKEELDKYFKNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIKWPWYIWLGFIAGLIAIVMVTIMLCCMTSCCSCLKGCCSCGSCCKFDEDDSEPVLKGVKLHYT*>143HKKQQKLDGK*|22000C>CA+23063A>T+23271C>A+23403A>G+23604C>A+23709C>T+24506T>G+24914G>C   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:21:10:39:9:0:20:0.409091:1
MN908947.3  23063   .   A   T   .   PASS    DP=235;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:235:104:39:1:1
MN908947.3  23271   .   C   A   .   PASS    DP=209;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:209:132:41:1:1
MN908947.3  23403   .   A   G   .   PASS    DP=199;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:199:137:44:1:1
MN908947.3  23604   .   C   A   .   PASS    DP=93;BCSQ=@22000   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:0:77:92:26:39:0.989247:1
MN908947.3  23709   .   C   T   .   PASS    DP=102;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:102:25:44:1:1
MN908947.3  24506   .   T   G   .   PASS    DP=120;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:120:86:41:1:1
MN908947.3  24914   .   G   C   .   PASS    DP=273;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:273:95:39:1:1
MN908947.3  27972   .   C   T   .   PASS    DP=513;BCSQ=stop_gained|ORF8|ENSSAST00005000008|protein_coding|+|27Q>27*|27972C>T   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:513:96:44:1:1
MN908947.3  28048   .   G   T   .   PASS    DP=744;BCSQ=*missense|ORF8|ENSSAST00005000008|protein_coding|+|52R>52I|28048G>T GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:744:188:43:1:1
MN908947.3  28111   .   A   G   .   PASS    DP=839;BCSQ=*missense|ORF8|ENSSAST00005000008|protein_coding|+|73Y>73C|28111A>G GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:839:212:42:1:1
MN908947.3  28270   .   TA  T   .   PASS    DP=785  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ   1:731:344:41:686:0:20:0.873885
MN908947.3  28280   .   G   C   .   PASS    DP=618;BCSQ=missense|N|ENSSAST00005000005|protein_coding|+|3D>3L|28280G>C+28281A>T+28282T>A GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:618:293:40:1:1
MN908947.3  28281   .   A   T   .   PASS    DP=620;BCSQ=@28280  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:620:297:41:1:1
MN908947.3  28282   .   T   A   .   PASS    DP=634;BCSQ=@28280  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:634:309:41:1:1
MN908947.3  28881   .   G   A   .   PASS    DP=740;BCSQ=missense|N|ENSSAST00005000005|protein_coding|+|203R>203K|28881G>A+28882G>A  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:1:33:739:446:42:0.998649:1
MN908947.3  28882   .   G   A   .   PASS    DP=744;BCSQ=@28881  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:1:33:743:452:41:0.998656:1
MN908947.3  28883   .   G   C   .   PASS    DP=748;BCSQ=missense|N|ENSSAST00005000005|protein_coding|+|204G>204R|28883G>C   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:1:33:747:454:42:0.998663:1
MN908947.3  28977   .   C   T   .   PASS    DP=501;BCSQ=missense|N|ENSSAST00005000005|protein_coding|+|235S>235F|28977C>T   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:501:351:40:1:1

Note the variant at position 22000:

REGION  POS REF ALT REF_DP  REF_RV  REF_QUAL    ALT_DP  ALT_RV  ALT_QUAL    ALT_FREQ    TOTAL_DP    PVAL    PASS    GFF_FEATURE REF_CODON   REF_AA  ALT_CODON   ALT_AA
MN908947.3  22000   C   +A  21  10  39  9   0   20  0.409091    22  0.00470272  TRUE    NA  NA  NA  NA  NA

...is included in the .csq.vcf file despite having ALT_FREQ of 0.409091. This introduces a frameshift in the S gene that interferes with consequence assignment downstream:

MN908947.3  22000   .   C   CA  .   PASS    DP=22;BCSQ=*frameshift|S|ENSSAST00005000004|protein_coding|+|146HKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLEGKQGNFKNLREFVFKNIDGYFKIYSKHTPINLVRDLPQGFSALEPLVDLPIGINITRFQTLLALHRSYLTPGDSSSGWTAGAAAYYVGYLQPRTFLLKYNENGTITDAVDCALDPLSETKCTLKSFTVEKGIYQTSNFRVQPTESIVRFPNITNLCPFGEVFNATRFASVYAWNRKRISNCVADYSVLYNSASFSTFKCYGVSPTKLNDLCFTNVYADSFVIRGDEVRQIAPGQTGKIADYNYKLPDDFTGCVIAWNSNNLDSKVGGNYNYLYRLFRKSNLKPFERDISTEIYQAGSTPCNGVEGFNCYFPLQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVCGPKKSTNLVKNKCVNFNFNGLTGTGVLTESNKKFLPFQQFGRDIADTTDAVRDPQTLEILDITPCSFGGVSVITPGTNTSNQVAVLYQDVNCTEVPVAIHADQLTPTWRVYSTGSNVFQTRAGCLIGAEHVNNSYECDIPIGAGICASYQTQTNSPRRARSVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSMTKTSVDCTMYICGDSTECSNLLLQYGSFCTQLNRALTGIAVEQDKNTQEVFAQVKQIYKTPPIKDFGGFNFSQILPDPSKPSKRSFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVLPPLLTDEMIAQYTSALLAGTITSGWTFGAGAALQIPFAMQMAYRFNGIGVTQNVLYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVEAEVQIDRLITGRLQSLQTYVTQQLIRAAEIRASANLAATKMSECVLGQSKRVDFCGKGYHLMSFPQSAPHGVVFLHVTYVPAQEKNFTTAPAICHDGKAHFPREGVFVSNGTHWFVTQRNFYEPQIITTDNTFVSGNCDVVIGIVNNTVYDPLQPELDSFKEELDKYFKNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIKWPWYIWLGFIAGLIAIVMVTIMLCCMTSCCSCLKGCCSCGSCCKFDEDDSEPVLKGVKLHYT*>143HKKQQKLDGK*|22000C>CA+23063A>T+23271C>A+23403A>G+23604C>A+23709C>T+24506T>G+24914G>C   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:21:10:39:9:0:20:0.409091:1
MN908947.3  23063   .   A   T   .   PASS    DP=235;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:235:104:39:1:1
MN908947.3  23271   .   C   A   .   PASS    DP=209;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:209:132:41:1:1
MN908947.3  23403   .   A   G   .   PASS    DP=199;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:199:137:44:1:1
MN908947.3  23604   .   C   A   .   PASS    DP=93;BCSQ=@22000   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:1:0:77:92:26:39:0.989247:1
MN908947.3  23709   .   C   T   .   PASS    DP=102;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:102:25:44:1:1
MN908947.3  24506   .   T   G   .   PASS    DP=120;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:120:86:41:1:1
MN908947.3  24914   .   G   C   .   PASS    DP=273;BCSQ=@22000  GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:273:95:39:1:1
MN908947.3  27972   .   C   T   .   PASS    DP=513;BCSQ=stop_gained|ORF8|ENSSAST00005000008|protein_coding|+|27Q>27*|27972C>T   GT:REF_DP:REF_RV:REF_QUAL:ALT_DP:ALT_RV:ALT_QUAL:ALT_FREQ:BCSQ  1:0:0:0:513:96:44:1:1
MN908947.3  28048   .   G   T   .   PASS    DP=744;

resulting SAMPLE.variants.csv file is:

sampleID,gene,aa_var,dna_var
SAMPLE,ORF1ab,Syn.216S,C913T
SAMPLE,ORF1ab,Syn.924F,C3037T
SAMPLE,ORF1ab,T1001I,C3267T
SAMPLE,ORF1ab,A1708D,C5388A
SAMPLE,ORF1ab,Syn.1907F,C5986T
SAMPLE,ORF1ab,I2230T,T6954C
SAMPLE,ORF1ab,LSGF3674L,GTCTGGTTTT11287G
SAMPLE,ORF1ab,Syn.4715L,C14408T
SAMPLE,ORF1ab,P4804L,C14676T
SAMPLE,ORF1ab,T5005I,C15279T
SAMPLE,ORF1ab,L5304P,T16176C
SAMPLE,S,VY143V,TTTA21990T
SAMPLE,ORF8,Q27*,C27972T
SAMPLE,ORF8,R52I,G28048T
SAMPLE,ORF8,Y73C,A28111G
SAMPLE,N,D3L,G28280C
SAMPLE,N,R203K,G28881A
SAMPLE,N,G204R,G28883C
SAMPLE,N,S235F,C28977T

...as a result, the VOC-202012/01 variant is not called.

Should variants be excluded from the .csq.vcf file if their ALT_FREQ is below --allele_freq_threshold? I've implemented that in another repo where I've extracted out the typing portion of this pipeline and it does recover the correct variant calling for this example.

https://github.com/BCCDC-PHL/type-variants-nf/pull/4/files

m-bull commented 3 years ago

Good catch - I guess the problem is that I run bcftools csq on the unfiltered vcf (which already breaks calls downstream of a low-AF frameshift), then filter the broken calls for AF later - by this point the damage is already done.

I think I'd like to output the un-AF-filtered vcf somewhere, but as in https://github.com/BCCDC-PHL/type-variants-nf/pull/4/files, i'll add the AF filter to the vcf before calling bcftools csq - does that seem like an OK solution?

dfornika commented 3 years ago

That sounds good to me.

ahmedmagds commented 3 years ago

@dfornika It is my first time to use the pipeline and I tried to annotate using the -gff option but I got an error. I simply downloaded the reference.gff file from NCBI so not sure if this is the right file or there is something else that I need to use. It would be great if you can share the gff file that you use or guide me where I can find it. Thanks

dfornika commented 3 years ago

@ahmedmagds look in the typing directory of this repository

ahmedmagds commented 3 years ago

@dfornika Thank you!

pmenzel commented 3 years ago

Having the same issue with a low-frequency frameshift upstream of N501Y etc, which in turn are not called.

pmenzel commented 3 years ago

Good catch - I guess the problem is that I run bcftools csq on the unfiltered vcf (which already breaks calls downstream of a low-AF frameshift), then filter the broken calls for AF later - by this point the damage is already done.

I think I'd like to output the un-AF-filtered vcf somewhere, but as in https://github.com/BCCDC-PHL/type-variants-nf/pull/4/files, i'll add the AF filter to the vcf before calling bcftools csq - does that seem like an OK solution?

Yeah, that would be great to have the original vcf file with all detected variants kept somewhere.