Mixing of 0-based and 1-based position in the outputed vcf file of DiscoSNP++

[jdalino]

Hello Pierre,

I use the bioconda discoSNP++ version (2.5.4).

I found a bug in the given positions of the SNPs in the vcf output of discoSNP++ (except if I missunderstood something).

When there are several SNPs on the same sequence, their positions are 1-based, and where there is only one SNP, it is 0-based.

Here are an extract of the discoRes_k_31_c_3_D_100_P_3_b_0_coherent.vcf:

SNP_higher_path_99998   55      99998_1 A       G       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99998   75      99998_2 G       A       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99998   101     99998_3 A       C       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99996   684     99996   A       T       .       .       Ty=SNP;Rk=1;UL=192;UR=277;CL=654;CR=778;Genome=.;Sd=.   GT:DP:PL:AD:HQ  1/1:19:384,61,5:0,19:0,70       0/0:13:5,43,264:13,0:72,0
SNP_higher_path_99994   209     99994   A       T       .       .       Ty=SNP;Rk=1;UL=57;UR=10;CL=179;CR=252;Genome=.;Sd=.     GT:DP:PL:AD:HQ  0/0:7:5,25,144:7,0:70,0 1/1:12:244,40,5:0,12:0,72
SNP_higher_path_99991   205     99991   C       G       .       .       Ty=SNP;Rk=1;UL=115;UR=173;CL=175;CR=205;Genome=.;Sd=.   GT:DP:PL:AD:HQ  0/0:18:5,58,364:18,0:70,0       1/1:22:444,70,5:0,22:0,73
SNP_higher_path_99990   94      99990_1 A       G       .       .       Ty=SNP;Rk=1;UL=63;UR=109;CL=63;CR=287;Genome=.;Sd=.     GT:DP:PL:AD:HQ  1|1:17:344,55,5:0,17:0,69       0|0:14:5,46,284:14,0:71,0
SNP_higher_path_99990   95      99990_2 G       C       .       .       Ty=SNP;Rk=1;UL=63;UR=109;CL=63;CR=287;Genome=.;Sd=.     GT:DP:PL:AD:HQ  1|1:17:344,55,5:0,17:0,69       0|0:14:5,46,284:14,0:71,0

And the corresponding discoRes_k_31_c_3_D_100_P_3_b_0_coherent.fa:

>SNP_higher_path_99998
ataactgcaccattttcaacccaaACACTACTCCACAAGTAAAAAGTCACCCAAAACAAATCAAGTTCAAACACGAACAAGCGACGGCAACCTCAACAGTAATTATGGCTCTCATTAGCAACTAACGTTCGcacaaggttca
>SNP_lower_path_99998
ataactgcaccattttcaacccaaACACTACTCCACAAGTAAAAAGTCACCCAAGACAAATCAAGTTCAAACACAAACAAGCGACGGCAACCTCAACAGTCATTATGGCTCTCATTAGCAACTAACGTTCGcacaaggttca

>SNP_higher_path_99996
agcatggatggaagtgccttgaatagctgtacttggtggtgcctcctagtgtatatgtcgtctttgattatccctgagaacagaaaaagaggaacaatttatggggagtggtgttgttatttcagatgatgctgaggaagccttttctcgtgctctcaagcttttgttggtttgtgattcttctccttctcaatcttgtgtctgatgtttgtgacttgagctattatttatacgcttcaacgcacatggttaatctttttgtttcctgcatcttgtagagtaaaagcttggtgaatgatggccagcatgtcactcttgtccacagtggagcacgaccaatctggcagaaatcatgttgagtagttggcatctaggttttattttgtaatttattcttacatgcagttgttttttttagctcaatttattctttctgcgatgatgagatatagaagttccttgagcgtggtttattatctttcaaagttatctgtggatgaaagattgcagttagatatcttttttgttgcattatagtttggatttcaccgtgtctagtcatcatcatctcagttctttctggattttctagttttactagacctctgaaacaaatacatcaacattgcggataaaacaaaacaatagaaaggtATTCTCTTTTGTAGTCTCTAGTGTAGTGTTACTTTTCTCAGAATTTGCATTTGATGACCTGaatccttttcggttagtcatattttactagtaaatgataaacttatgttctttgttgaactaattttgtgaggaaaatgcttcgagctctgaattgcatatgaaattaatgagatcttttaatacatggaactttttgggtacagtagtacttttcaacagtttcctggatacaggataaaccaatccaatttacatatcaatccaaaatgcttgtcaatttagaactaatatgatggaacaaacttatgcagcaacaacgagaaacgatgaacctagaaagtaagagttgtccagttccctggagacttctgcagagctgggcattttgaattgattgagttggagcattgccgagggaggaaatccaaataccaatcataaactctactcagaatgttgttggggtggctgcacatcttcctcttattcatatggttttgcaaaacatgcctaaggtatgatccatgtcaggtacgttagacttcacacttgtcctctgcctcatgttcctgatacttctcgctatgcttaatattcttgtaccacttagcacatgctatgtatgcagccaaatctattaaggtcaaaccagctaatagaaaatagaatctgtctaaatgacctttgttaaggtttcctggtatccatccaggcatgttgtcagttgctgatattttcatcactatgctcacgagaaagctactcacgtagttccccaatgatatggatgtcatacatagtgcacttccaaagcttttaagtccatcaggtgcttg
>SNP_lower_path_99996
agcatggatggaagtgccttgaatagctgtacttggtggtgcctcctagtgtatatgtcgtctttgattatccctgagaacagaaaaagaggaacaatttatggggagtggtgttgttatttcagatgatgctgaggaagccttttctcgtgctctcaagcttttgttggtttgtgattcttctccttctcaatcttgtgtctgatgtttgtgacttgagctattatttatacgcttcaacgcacatggttaatctttttgtttcctgcatcttgtagagtaaaagcttggtgaatgatggccagcatgtcactcttgtccacagtggagcacgaccaatctggcagaaatcatgttgagtagttggcatctaggttttattttgtaatttattcttacatgcagttgttttttttagctcaatttattctttctgcgatgatgagatatagaagttccttgagcgtggtttattatctttcaaagttatctgtggatgaaagattgcagttagatatcttttttgttgcattatagtttggatttcaccgtgtctagtcatcatcatctcagttctttctggattttctagttttactagacctctgaaacaaatacatcaacattgcggataaaacaaaacaatagaaaggtATTCTCTTTTGTAGTCTCTAGTGTAGTGTTTCTTTTCTCAGAATTTGCATTTGATGACCTGaatccttttcggttagtcatattttactagtaaatgataaacttatgttctttgttgaactaattttgtgaggaaaatgcttcgagctctgaattgcatatgaaattaatgagatcttttaatacatggaactttttgggtacagtagtacttttcaacagtttcctggatacaggataaaccaatccaatttacatatcaatccaaaatgcttgtcaatttagaactaatatgatggaacaaacttatgcagcaacaacgagaaacgatgaacctagaaagtaagagttgtccagttccctggagacttctgcagagctgggcattttgaattgattgagttggagcattgccgagggaggaaatccaaataccaatcataaactctactcagaatgttgttggggtggctgcacatcttcctcttattcatatggttttgcaaaacatgcctaaggtatgatccatgtcaggtacgttagacttcacacttgtcctctgcctcatgttcctgatacttctcgctatgcttaatattcttgtaccacttagcacatgctatgtatgcagccaaatctattaaggtcaaaccagctaatagaaaatagaatctgtctaaatgacctttgttaaggtttcctggtatccatccaggcatgttgtcagttgctgatattttcatcactatgctcacgagaaagctactcacgtagttccccaatgatatggatgtcatacatagtgcacttccaaagcttttaagtccatcaggtgcttg

>SNP_higher_path_99994
ttttttcaattttggttgtgttgttttttttgagtgtagcattgcaggcaaaaacaatatacagttgggaatacatgtcaccagacggactggctttaaactccaagtggaatgaagctgagaaatatatctgcaatcctttatcaggggaagtcccattagaatgtttatctgcaaaaACACTAAGTGGAAGATCATTTCGACAATTAACCAATAAAATCACCATGTCTGCACCTTTGAtttatccttcacaatatcagtgtgctcgacgattcaatccaaaacctcttacaaaagtagtacctcatgtgcctccacatcaactacaaatacaaattccaagcaaaggtatatatacatatatacatttcctctctttccatttcctcattactacatatatatatacatactccgatttctaacttgtacgtttttattattcaatttagagagaattggtagcattacaagagatgtgggaacacagag
>SNP_lower_path_99994
ttttttcaattttggttgtgttgttttttttgagtgtagcattgcaggcaaaaacaatatacagttgggaatacatgtcaccagacggactggctttaaactccaagtggaatgaagctgagaaatatatctgcaatcctttatcaggggaagtcccattagaatgtttatctgcaaaaACACTAAGTGGAAGATCATTTCGACAATTATCCAATAAAATCACCATGTCTGCACCTTTGAtttatccttcacaatatcagtgtgctcgacgattcaatccaaaacctcttacaaaagtagtacctcatgtgcctccacatcaactacaaatacaaattccaagcaaaggtatatatacatatatacatttcctctctttccatttcctcattactacatatatatatacatactccgatttctaacttgtacgtttttattattcaatttagagagaattggtagcattacaagagatgtgggaacacagag

>SNP_higher_path_99991
cattctccctaagcgcctctcttggcccagtcaccatcggttcattaccagcaagcatatttgttatcaattgcttgatttcatccatggcaccaacaactttcgaatccacacttctgctcagttccgagatggatccaacaacaccttccatggatcctttcaaattgtcaatCTCCCTGCGCATTTCCTGCAAGCTAGTGTTCTGATCACCCACACTCACAGGTCGAAAGAGCtctgaagctctgataccaacattacaactagcttcagagaacttctaattgattaagacttgctaattctttacttttctcaacaacgaacttcaactcattccttccaatctttttatagagtattacaagtaataggaaatttggaggaaaaaacagtaacaactaccttccgaaatgaacagtaaacttctaacagctttcc
>SNP_lower_path_99991
cattctccctaagcgcctctcttggcccagtcaccatcggttcattaccagcaagcatatttgttatcaattgcttgatttcatccatggcaccaacaactttcgaatccacacttctgctcagttccgagatggatccaacaacaccttccatggatcctttcaaattgtcaatCTCCCTGCGCATTTCCTGCAAGCTAGTGTTGTGATCACCCACACTCACAGGTCGAAAGAGCtctgaagctctgataccaacattacaactagcttcagagaacttctaattgattaagacttgctaattctttacttttctcaacaacgaacttcaactcattccttccaatctttttatagagtattacaagtaataggaaatttggaggaaaaaacagtaacaactaccttccgaaatgaacagtaaacttctaacagctttcc

>SNP_higher_path_99990
ccttacttttaaacttttgacaattcatactcgtgatacgttttaattccatatggctatcaaATCCATCAATAGTATATGGACTCTAGTGTTAGCTTTTATCTCACCTTTTGAATTTCTCAACATcttctacagctactgaagatggagtatagatattgctttcatctcaatatctataaagtaaatattgaagtattgagtatcgtatcgaaataatttcttatatatacaactaatacattatacacctgcaaacgcttaagctcatatagtcttcacatataataattagctatatactatttgtcaataaccaaaatattctgaaactcacattttcaattatttcatattacttgcgtttgcctttccaattatgttcatgattccaacactttaaaaaaaaaaaa
>SNP_lower_path_99990
ccttacttttaaacttttgacaattcatactcgtgatacgttttaattccatatggctatcaaATCCATCAATAGTATATGGACTCTAGTGTTGCTTTTTATCTCACCTTTTGAATTTCTCAACATcttctacagctactgaagatggagtatagatattgctttcatctcaatatctataaagtaaatattgaagtattgagtatcgtatcgaaataatttcttatatatacaactaatacattatacacctgcaaacgcttaagctcatatagtcttcacatataataattagctatatactatttgtcaataaccaaaatattctgaaactcacattttcaattatttcatattacttgcgtttgcctttccaattatgttcatgattccaacactttaaaaaaaaaaaa

The real position are:

SNP_higher_path_99998   55      99998_1 A       G       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99998   75      99998_2 G       A       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99998   101     99998_3 A       C       .       .       Ty=SNP;Rk=1;UL=24;UR=11;CL=24;CR=11;Genome=.;Sd=.       GT:DP:PL:AD:HQ  1|1:18:364,58,5:0,18:0,67       0|0:37:6,116,744:37,0:72,0
SNP_higher_path_99996   685     99996   A       T       .       .       Ty=SNP;Rk=1;UL=192;UR=277;CL=654;CR=778;Genome=.;Sd=.   GT:DP:PL:AD:HQ  1/1:19:384,61,5:0,19:0,70       0/0:13:5,43,264:13,0:72,0
SNP_higher_path_99994   210     99994   A       T       .       .       Ty=SNP;Rk=1;UL=57;UR=10;CL=179;CR=252;Genome=.;Sd=.     GT:DP:PL:AD:HQ  0/0:7:5,25,144:7,0:70,0 1/1:12:244,40,5:0,12:0,72
SNP_higher_path_99991   206     99991   C       G       .       .       Ty=SNP;Rk=1;UL=115;UR=173;CL=175;CR=205;Genome=.;Sd=.   GT:DP:PL:AD:HQ  0/0:18:5,58,364:18,0:70,0       1/1:22:444,70,5:0,22:0,73
SNP_higher_path_99990   94      99990_1 A       G       .       .       Ty=SNP;Rk=1;UL=63;UR=109;CL=63;CR=287;Genome=.;Sd=.     GT:DP:PL:AD:HQ  1|1:17:344,55,5:0,17:0,69       0|0:14:5,46,284:14,0:71,0
SNP_higher_path_99990   95      99990_2 G       C       .       .       Ty=SNP;Rk=1;UL=63;UR=109;CL=63;CR=287;Genome=.;Sd=.     GT:DP:PL:AD:HQ  1|1:17:344,55,5:0,17:0,69       0|0:14:5,46,284:14,0:71,0

Each SNP where the ID is XXXXX is 0-based, and when it's XXXXX_X it's 1-based.

I confirmed this on a 1 million SNPs vcf file.
My command line was
run_discoSnp++.sh --max_threads 16 -T -r ../file.fof

Do you have an idea of where this come from?

Best regards,
Jordi

[pierrepeterlongo]

Dear Jordi,

Thanks for this issue. I'm sorry I'm seeing it only now.
I'll have a deeper look in a few days, as soon I find a few minutes.

Best,
Pierre

[pierrepeterlongo]

Dear Jordi,

I tried to reproduce the problem you mention with a toy example. With no success.

I created a file with two sequences identical but with 3 mutations, two closes and one isolated. These are lower case characters in the second sequence.
sequences.fa:

>tibo
TGTATCCGCATTTGATGCTACCGTGGATGAGTCAGCGTCGAGCACGCGGCACTTATTGCATGAGTAGGGTTGACTAAGAGCCGTTAGATGCCTCGCTGTACTAATAGTTGTCGACAGATCGTCAAGATTAGAAAACGGTACCAGCATTTTCGGAGGTTCTCTAACTAGTATGGATAGCCGTGTCTTCACTGTGCTGCGGC
>tibo2
TGTATCCGCATTTGATGCTACCGTGGATGAGTgAGCGaCGAGCACGCGGCACTTATTGCATGAGTAGGGTTGACTAAGAGCCGTTAGATGCCTCGCTGTACTAATAGTTGTCGACAGATCGTCAAGATTAGAAAACGGTACCAGCATTTTCGGAGGTTCgCTAACTAGTATGGATAGCCGTGTCTTCACTGTGCTGCGGC

I also used "tibo" as the reference genome:
ref.fa

>tibo
TGTATCCGCATTTGATGCTACCGTGGATGAGTCAGCGTCGAGCACGCGGCACTTATTGCATGAGTAGGGTTGACTAAGAGCCGTTAGATGCCTCGCTGTACTAATAGTTGTCGACAGATCGTCAAGATTAGAAAACGGTACCAGCATTTTCGGAGGTTCTCTAACTAGTATGGATAGCCGTGTCTTCACTGTGCTGCGGC

Mutations are located at positions 33 38 and 160.
I then ran the following commands with disco commit version "1be346567efc4db5cdda6c19d1d1fb07a3310385" (actual version).

ls sequences.fa > fof
run_discoSnp++.sh -r fof -c 1 -G ref.f

The resulting vcf is

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	G1
tibo	33	1_2	C	G	.	PASS	TySNP;Rk=0;UL=91;UR=2;CL=.;CR=.;Genome=.;Sd=-1	GT:DP:PL:AD:HQ	0|1:2:21,7,21:1,1:0,0
tibo	38	1_1	T	A	.	PASS	TySNP;Rk=0;UL=91;UR=2;CL=.;CR=.;Genome=.;Sd=-1	GT:DP:PL:AD:HQ	0|1:2:21,7,21:1,1:0,0
tibo	160	2	T	G	.	PASS	Ty=SNP;Rk=0;UL=91;UR=10;CL=.;CR=.;Genome=T;Sd=1	GT:DP:PL:AD:HQ	0/1:2:21,7,21:1,1:0,0

We retrieve positions 33, 38, and 160 with no difference between close SNPs and isolated SNPs.
Could you please check if this problem persists with the current version?
Feel free to re-open this issue if this is not the case.

I hope this helps.

Pierre

[jdalino]

Hello Pierre,

Thanks for your answer.

Indeed it works when using a reference sequence, but in my analysis I don't have one so I tried your data without it with:
run_discoSnp++.sh -c 1 -r fof

And then the problem appears again. (I used the latest github version)

The unitigs found by discosnp in the discoRes_k_31_c_1_D_100_P_3_b_0_coherent.fa file are:

>SNP_higher_path_2|P_1:30_T/G|high|nb_pol_1|left_unitig_length_91|right_unitig_length_10|C1_1|Q1_0|G1_0/1:21,7,21|rank_0
cgagcacgcggcacttattgcatgagtagggttgactaagagccgttagatgcctcgctgtactaatagttgtcgacagatcgtcaagattAGAAAACGGTACCAGCATTTTCGGAGGTTCTCTAACTAGTATGGATAGCCGTGTCTTCACTgtgctgcggc

>SNP_lower_path_2|P_1:30_T/G|high|nb_pol_1|left_unitig_length_91|right_unitig_length_10|C1_1|Q1_0|G1_0/1:21,7,21|rank_0
cgagcacgcggcacttattgcatgagtagggttgactaagagccgttagatgcctcgctgtactaatagttgtcgacagatcgtcaagattAGAAAACGGTACCAGCATTTTCGGAGGTTCGCTAACTAGTATGGATAGCCGTGTCTTCACTgtgctgcggc

>SNP_higher_path_1|P_1:30_A/T,P_2:35_G/C|high|nb_pol_2|left_unitig_length_91|right_unitig_length_2|C1_1|Q1_0|G1_0/1:21,7,21|rank_0
gaacctccgaaaatgctggtaccgttttctaatcttgacgatctgtcgacaactattagtacagcgaggcatctaacggctcttagtcaacCCTACTCATGCAATAAGTGCCGCGTGCTCGACGCTGACTCATCCACGGTAGCATCAAATGCGGATAca

>SNP_lower_path_1|P_1:30_A/T,P_2:35_G/C|high|nb_pol_2|left_unitig_length_91|right_unitig_length_2|C1_1|Q1_0|G1_0/1:21,7,21|rank_0
gaacctccgaaaatgctggtaccgttttctaatcttgacgatctgtcgacaactattagtacagcgaggcatctaacggctcttagtcaacCCTACTCATGCAATAAGTGCCGCGTGCTCGTCGCTCACTCATCCACGGTAGCATCAAATGCGGATAca

And the vcf discoRes_k_31_c_1_D_100_P_3_b_0_coherent.vcf is giving 3 snps:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  G1
SNP_lower_path_2        121     2       G       T       .       .       Ty=SNP;Rk=0;UL=91;UR=10;CL=.;CR=.;Genome=.;Sd=. GT:DP:PL:AD:HQ  0/1:2:21,7,21:1,1:0,0
SNP_higher_path_1       122     1_1     A       T       .       .       TySNP;Rk=0;UL=91;UR=2;CL=.;CR=.;Genome=.;Sd=.   GT:DP:PL:AD:HQ  0|1:2:21,7,21:1,1:0,0
SNP_higher_path_1       127     1_2     G       C       .       .       TySNP;Rk=0;UL=91;UR=2;CL=.;CR=.;Genome=.;Sd=.   GT:DP:PL:AD:HQ  0|1:2:21,7,21:1,1:0,0

Identically as I observed before with my dataset, the isolated SNP position is incorrect by one base (SNP_lower_path_2 121). If you look at the sequence of the SNP_lower_path_2 unitig, the nucleotide in position 121 is a C and not a G, the G being positioned at position 122 (or 121 in a 0-based system). The positions of the 2 others SNPs are correct (in a 1-based system).

From what I have seen in the code of VCF_creator.py I think you may need to add a +1 to the mappingPositionCouple value somewhere (maybe in the WhichPathIsTheRef() function of the SNP class (and maybe to the INDEL one too)). I have the impression that the calculation are correct for the SNPCLOSE class because it don't use this value to determine the position and only add the corrected position and the unitig length, but I didn't had the time to look further.

Tell me if you can reproduce the problem or if you want me to run some tests and maybe propose a pull request.

(And for the issue I do no have the right to reopen it myself so I will leave it to you if you found it appropriate)

Best regards,
Jordi

[pierrepeterlongo]

Thanks again Jordi for your detailed and precious feedback.

I found and corrected the bug for close SNPs. Everything should now be zero-based. You may pull and retry (no need to recompile).

I let the issue open as the bug also exists for INDELS. I have to explore this code.
(When I have a few days free I'll recode from scratch VCF_creator.)

[pierrepeterlongo]

The bug is also fixed for unmapped indels

[jdalino]

Thanks a lot Pierre for fixing this.

I see that you have made everything 0-based, but the vcf convention state that vcf are 1-based. To avoid confusion with software working on vcf files do you think you could produce 1-based vcf instead? I think this could prevent many mistakes in the future.

Jordi