Different output between vg snarls and vg deconstruct -a

[AliceLaurent]

1. What were you trying to do?
I was trying to compute the snarls on a little graph by running both vg snarls and vg deconstruct -a

2. What did you want to happen?
I thought vg snarls and vg deconstruct -a would output the same snarls

3. What actually happened?
I had two different output, two snarls are "missing" in vg deconstruct output

Output of vg snarls :
{"end": {"node_id": "16"}, "end_self_reachable": true, "start": {"node_id": "1"}, "start_end_reachable": true, "start_self_reachable": true}
{"directed_acyclic_net_graph": true, "end": {"node_id": "15"}, "end_self_reachable": true, "parent": {"end": {"node_id": "16"}, "start": {"node_id": "1"}}, "start": {"node_id": "11"}, "start_end_reachable": true}
{"directed_acyclic_net_graph": true, "end": {"node_id": "8"}, "end_self_reachable": true, "parent": {"end": {"node_id": "16"}, "start": {"node_id": "1"}}, "start": {"node_id": "2"}, "start_end_reachable": true}

Ouput of vg deconstruct :
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT inv
ref 100 >1>16 CTTGACTTAGGCCAATACCTTTTTGTGTCTTGACCCCTGCAAATACGTTAGTGGTTCGTGACCGACTTCAGGGTCCCTGACCT CTCCCAAAGGTATTGGCCTAAGTCAACAAAACACAAACTAACGTATTTGCAGGGGTGTGACGAACCAGGTCAGGGACCCTGAAGTCGT 60 . AC=1;AF=1;AN=1;AT=>1>2>3>4>8>10>9>11>13>14>15>16,>1>2<5<6<4<3<2<7<9<10<12<15<14<13>15>16;NS=1;LV=0;RC=ref;RS=100;RD=183 GT 1

Also got "[vg deconstruct] Using flat processing" in stderr
I had this case for several graph, I would like to understand why if there is a reason.

5. What data and command can the vg dev team use to make the problem happen?

Here is the gfa graph I used to produce this output :

H VN:Z:1.0
S 1 CTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAACCTGGGTAAAC
S 2 T
S 3 TGACTTAGGC
S 4 CAATACCTTT
S 5 GG
S 6 G
S 7 TTGTGTTTTG
S 8 TTGTGTCTTG
S 9 ATACGTTAGT
S 10 ACCCCTGCAA
S 11 GGTTCGTGAC
S 12 GGTTCGTCAC
S 13 CGACTTCAGG
S 14 GTCCCTGACC
S 15 T
S 16 CGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCGCGAGTAGGCG
L 2 + 3 + 0M
L 9 - 10 - 0M
L 12 - 15 - 0M
L 13 + 14 + 0M
L 10 + 9 + 0M
L 4 - 3 - 0M
L 11 + 13 + 0M
L 8 + 10 + 0M
L 3 - 2 - 0M
L 4 + 8 + 0M
L 9 + 11 + 0M
L 15 - 14 - 0M
L 3 + 4 + 0M
L 13 - 15 + 0M
L 5 - 6 - 0M
L 6 - 4 - 0M
L 14 + 15 + 0M
L 10 - 12 - 0M
L 2 - 7 - 0M
L 1 + 2 + 0M
L 14 - 13 - 0M
L 7 - 9 - 0M
L 15 + 16 + 0M
L 2 + 5 - 0M
P ref 1+,2+,3+,4+,8+,10+,9+,11+,13+,14+,15+,16+ ,,,,,,,,,,,
P inv 1+,2+,5-,6-,4-,3-,2-,7-,9-,10-,12-,15-,14-,13-,15+,16+ ,,,,,,,,,,,,,,,

The command I used are :
vg convert -f -W graph.gfa > graph.vg
vg snarls graph.vg > snarl_ouput.snarls
vg view -R snarl_output.snarls > snarls.json
vg deconstruct -a -p ref graph.vg > snarls.vcf (also tried with -r vg_snarls_output.snarls)

6. What does running vg version say?

vg version v1.73.0 "Ducky"
Compiled with g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0 on Linux
Linked against libstd++ 20230528
Using HTSlib headers 101990, library 1.19.1-29-g3cfe8769
Built by anovak@courtyard.gi.ucsc.edu

[xchang1]

Hi Alice, I asked @jmonlong about this and he says that it might be because vg deconstruct only defines snarls based on alternate paths relative to the reference path. It looks like for both of the missing snarls, the alternate path doesn't traverse the snarl from start to end but from start to start or from end to end. I'm not sure if this is expected behavior or not, maybe @glennhickey would know?

[AliceLaurent]

Hi, thanks for your answer. I think it could be the reason, since it works also for the other examples that I have where I see a difference between vg snarls and vg deconstruct.

[jmonlong]

I think that deconstruct will do its best to represent the variants relative to the reference. In those cases, the missing snarls don't really have alternate "alleles" compared to the reference allele, say >11>13>14>15, that deconstruct could report. The other path (inv) enters and exits the snarl by the same boundary node 15 so it's not really doable to represent that in a typical VCF. I'm not sure what we would want as ALT sequence for example.

If we were to add an edge to make it possible to exit by the other boundary, say by adding this to your GFA:

L	14	+	14	-	0M
P	inv_sm	1+,2+,3+,4+,8+,10+,9+,11+,13+,14+,14-,13-,15+,16+	,,,,,,,,,,,

Then deconstruct would output that snarl because it now has another path from one boundary node to the other.

##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=CONFLICT,Number=.,Type=String,Description="Sample names for which there are multiple paths in the graph with conflicting alleles">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of alternate alleles in called genotypes">
##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=LV,Number=1,Type=Integer,Description="Level in the snarl tree (0=top level)">
##INFO=<ID=PS,Number=1,Type=String,Description="ID of variant corresponding to parent snarl">
##INFO=<ID=AT,Number=R,Type=String,Description="Allele Traversal as path in graph">
##contig=<ID=ref,length=292>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	inv	inv_sm
ref	100	>1>16	CTTGACTTAGGCCAATACCTTTTTGTGTCTTGACCCCTGCAAATACGTTAGTGGTTCGTGACCGACTTCAGGGTCCCTGACCT	CTCCCAAAGGTATTGGCCTAAGTCAACAAAACACAAACTAACGTATTTGCAGGGGTGTGACGAACCAGGTCAGGGACCCTGAAGTCGT,CTTGACTTAGGCCAATACCTTTTTGTGTCTTGACCCCTGCAAATACGTTAGTGGTTCGTGACCGACTTCAGGGTCCCTGACCGGTCAGGGACCCTGAAGTCGT	60	.	AC=1,1;AF=0.5,0.5;AN=2;AT=>1>2>3>4>8>10>9>11>13>14>15>16,>1>2<5<6<4<3<2<7<9<10<12<15<14<13>15>16,>1>2>3>4>8>10>9>11>13>14<14<13>15>16;NS=2;LV=0	GT	1	2
ref	161	>11>15	CCGACTTCAGGGTCCCTGACC	CCGACTTCAGGGTCCCTGACCGGTCAGGGACCCTGAAGTCG	60	.	AC=1;AF=1;AN=1;AT=>11>13>14>15,>11>13>14<14<13>15;NS=1;LV=1;PS=>1>16	GT	.1

By the way, I noticed that the snarls that are not deconstructable are missing one of the "end_self_reachable": true, "start_end_reachable": true, "start_self_reachable": true. Maybe a way to identify them in advance.

[AliceLaurent]

Thanks for this detailled explanation. I understand better why some snarls are missing.
For the "end_self_reachable": true, "start_end_reachable": true, "start_self_reachable": true, I tried it on my large examples and it works too. However, I have not seen any example were the "start_end_reachable":true is missing, only the two others are missing.