<commit_before>#!/usr/bin/env bash

BASH_TAP_ROOT=../deps/bash-tap
. ../deps/bash-tap/bash-tap-bootstrap

PATH=../bin:$PATH # for vg

plan tests 37

vg construct -r complex/c.fa -v complex/c.vcf.gz > c.vg
cat <(vg view c.vg | grep ^S | sort) <(vg view c.vg | grep L | uniq | wc -l) <(vg paths -v c.vg -E) > c.info

vg convert c.vg -x > c.xg
vg convert c.xg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout xg conversion"

rm -f c.xg c1.vg c1.info

vg convert c.vg -a > c.hg
vg convert c.hg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout hash-graph conversion"

rm -f c.hg c1.vg c1.info

vg convert c.vg -p > c.pg
vg convert c.pg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout packed-graph conversion"

rm -f c.pg c1.vg c1.info

vg convert c.vg -f > c.gfa
vg convert -g c.gfa -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout gfa conversion"

rm -f c.gfa c1.vg c1.info

vg convert c.vg -o > c.odgi
vg convert c.odgi -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout ODGI conversion"

rm -f c.vg c.odgi c1.vg c.info c1.info

# some less rigorous tests I made without noticing that the earlier ones had already been written
vg construct -r small/x.fa -v small/x.vcf.gz > x.vg
vg view x.vg > x.gfa

is "$(vg convert -a x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "hash graph conversion looks good"
is "$(vg convert -p x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "packed graph conversion looks good"
is "$(vg convert -v x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "vg conversion looks good"
is "$(vg convert -o x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "odgi conversion looks good"
is "$(vg convert -f x.vg | vg convert -g - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa conversion looks good"
is "$(vg convert -x x.vg | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "xg conversion looks good"

is "$(vg convert -g -a x.gfa | vg view - | wc -l)" "$(wc -l < x.gfa)" "on disk gfa conversion looks good"
is "$(cat x.gfa | vg convert -g -a - | vg view - | wc -l)" "$(wc -l < x.gfa)" "streaming gfa conversion looks good"
is "$(vg convert -g -x x.gfa | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa to xg conversion looks good"
is "$(vg convert -g -f x.gfa | vg convert -g - | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa to gfa conversion looks good"

rm x.vg x.gfa
rm -f c.vg c.pg c1.vg c.info c1.info

vg construct -r small/x.fa -v small/x.vcf.gz > x.vg
vg index x.vg -g x.gcsa
vg sim -x x.vg -n 10 -s 23 -a > sim.gam
vg map -x x.vg -g x.gcsa -G sim.gam > sim-rm.gam
vg convert x.vg -G sim-rm.gam -t 1 > sim-rm.gaf
vg convert x.vg -F sim-rm.gaf -t 1 | vg convert x.vg -G - -t 1 > sim-rm2.gaf
diff sim-rm.gaf sim-rm2.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf makes same gaf twice"

vg convert x.vg -G sim-rm.gam | vg convert x.vg -F - | vg convert x.vg -G - | sort > sim-rm2-mt-sort.gaf
sort sim-rm2.gaf > sim-rm2-sort.gaf
diff sim-rm2-sort.gaf sim-rm2-mt-sort.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf gives same result multithreaded as with -t 1"

vg convert x.vg -G sim-rm.gam | bgzip | vg convert x.vg -F - | vg convert x.vg -G - | sort > sim-rm2-mtbg-sort.gaf
diff sim-rm2-sort.gaf sim-rm2-mtbg-sort.gaf
is "$?" 0 "vg convert gam -> gaf.gz -> gam -> gaf gives same result multithreaded as with -t 1"

# some snps and indels
vg map -s "TAATGGATATGTTAAGCTTTTTTTTTCTTTGATTTATTTGAAAAGACGTTTGACAATCTATCGGGTAATGTGGGGAAA" -x x.vg -g x.gcsa > mut.gam
# reverse complement of above
vg map -s "TTTCCCCACATTACCCGATAGATTGTCAAACGTCTTTTCAAATAAATCAAAGAAAAAAAAAGCTTAACATATCCATTA" -x x.vg -g x.gcsa >> mut.gam

vg convert x.vg -G mut.gam -t 1 > mut.gaf
vg convert x.vg -F mut.gaf -t 1 > mut-back.gam
vg view -a mut.gam | jq .path > mut.path
vg view -a mut-back.gam | jq .path > mut-back.path
# Json comparison that is not order dependent: https://stackoverflow.com/a/31933234
is $(jq --argfile a mut.path --argfile b mut-back.path -n 'def post_recurse(f): def r: (f | select(. != null) | r), .; r; def post_recurse: post_recurse(.[]?); ($a | (post_recurse | arrays) |= sort) as $a | ($b | (post_recurse | arrays) |= sort) as $b | $a == $b') true "vg convert gam -> gaf -> gam produces same gam Paths with snps and indels"

vg view -a mut.gam | jq .sequence > mut.seq
vg view -a mut-back.gam | jq .sequence > mut-back.seq
diff mut.seq mut-back.seq
is "$?" 0 "vg convert gam -> gaf -> gam preserves sequence"

vg convert x.vg -G mut-back.gam -t 1 > mut-back.gaf
diff mut.gaf mut-back.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf makes same gaf twice in presence of indels and snps"

rm -f x.vg x.gcsa sim.gam sim-rm.gam sim-rm.gaf sim-rm2.gaf sim-rm2-mt-sort.gaf sim-rm2-mtbg-sort.gaf sim-rm2-sort.gaf mut.gam mut-back.gam mut.gaf mut-back.gaf mut.path mut-back.path mut.seq mut-back.seq

vg construct -r 1mb1kgp/z.fa -v 1mb1kgp/z.vcf.gz > z.vg 2> /dev/null
vg sim -n 10000 -s 23 -a -x z.vg > sim.gam
vg construct -r 1mb1kgp/z.fa > zflat.vg
vg index zflat.vg -g zflat.gcsa
vg map -x zflat.vg -g zflat.gcsa -G sim.gam > sim-map.gam
vg convert zflat.vg -G sim-map.gam | bgzip > sim-map.gaf.gz
vg convert zflat.vg -F sim-map.gaf.gz > sim-map-back.gam
vg view -a sim-map.gam | jq .sequence | sort > sim-map.sequence
vg view -a sim-map-back.gam | jq .sequence | sort > sim-map-back.sequence
diff sim-map.sequence sim-map-back.sequence
is "$?" 0 "vg convert gam -> gaf -> gam preserves sequences of 1mb1kgp simulated reads"

vg convert zflat.vg -G sim-map-back.gam | sort > sim-map-back.gaf
bgzip -dc sim-map.gaf.gz | sort > sim-map.gaf
diff sim-map-back.gaf sim-map.gaf
is "$?" 0 "vg convert gam -> gaf -> gam ->gaf makes same gaf each time on 1mb1kgp simulated reads"

printf '{"name": "split", "path": {"mapping": [{"edit": [{"from_length": 13, "to_length": 13}], "position": {"node_id": "1", "offset": "10"}}, {"edit": [{"from_length": 2, "to_length": 2}], "position": {"node_id": "3", "offset": "5"}}]}}' | vg view -JaG - > split.gam
vg convert zflat.vg -G split.gam > split.gaf
is "$(awk '{print $13}' split.gaf)" "cs:Z::13-CCAGTGCTC-GCATC:2" "split alignment converted using deletions to represent internal offsets"
vg convert zflat.vg -F split.gaf | vg convert zflat.vg -G - > split-back.gaf
diff split.gaf split-back.gaf
is "$?" 0 "vg convert gam -> gaf ->gam -> gaf makes same gaf each time for split alignment"

rm -f z.vg zflat.vg sim.gam sim-map.gam sim-map-back.gam sim-map.gaf.gz sim-map.sequence sim-map-back.sequence sim-map-back.gaf sim-map.gaf split.gam split.gaf split-back.gaf

printf "H\tVN:Z:1.0
S\t73333\tGGTGGGCGAGGACCTCCACACGTGTCACCA
S\t73368\tGCCCCT
S\t72943\tGGCGACTCTTCAGCAAGCCCCTCCACACGTGT
S\t72940\tC
S\t72941\tGGCCAGGT
S\t73255\tACTCTTCAGCAGGCCCCTCTGGT
S\t72942\tGGGCGAGGACCTCCACACGTGTCACCAGGCCA
S\t73318\tTCAGCA
S\t73367\tA
S\t73271\tA
S\t73289\tC
S\t73317\tC\n" | vg convert -g - -p > soft.pg
printf '{"annotation": {"fragment_length": 242, "fragment_length_distribution": "-I 561.110526 -D 141.152986", "mapq_applied_cap": 23.832780374978935, "mapq_extended_cap": 15, "mapq_uncapped": 10.325140756048304, "secondary_scores": [187.15265582416521, 182.4063994408188]}, "identity": 0.89682539682539686, "mapping_quality": 10, "name": "ERR903030.2067", "path": {"mapping": [{"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "72943", "offset": "26"}}, {"edit": [{"from_length": 32, "to_length": 32}], "position": {"is_reverse": true, "node_id": "72942"}, "rank": "1"}, {"edit": [{"from_length": 23, "to_length": 23}], "position": {"is_reverse": true, "node_id": "73255"}, "rank": "2"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73271"}, "rank": "3"}, {"edit": [{"from_length": 8, "to_length": 8}], "position": {"is_reverse": true, "node_id": "72941"}, "rank": "4"}, {"edit": [{"from_length": 7, "to_length": 7}, {"from_length": 1, "sequence": "C", "to_length": 1}, {"from_length": 2, "to_length": 2}, {"from_length": 1, "sequence": "G", "to_length": 1}, {"from_length": 19, "to_length": 19}], "position": {"is_reverse": true, "node_id": "73333"}, "rank": "5"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "72940"}, "rank": "6"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73289"}, "rank": "7"}, {"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "73368"}, "rank": "8"}, {"edit": [{"from_length": 1, "sequence": "G", "to_length": 1}], "position": {"is_reverse": true, "node_id": "73367"}, "rank": "9"}, {"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "73318"}, "rank": "10"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73317"}, "rank": "11"}, {"edit": [{"sequence": "GGGTGGCCTG", "to_length": 10}], "position": {"is_reverse": true, "node_id": "73317", "offset": "1"}, "rank": "12"}]}, "quality": "ISEhISEmJiUmJCYmJSYlJiYmJiYmJCYjJiYdIRsmIyUmJiYmJiElJg4PDx0PDyEkIhsYJCYiECQQJBAcGR0kHx0QGSQPJCImHSElHR0PJA4OGyQbIxwPIg0iHw8PGRwiIR0jAgICAgICAgICAgICAgICAgICAgICAgICAgIC", "sample_name": "HG00514_961a37c_gssw", "score": 106, "sequence": "GTCGCCTGGCCTGGTGACACGTGTGGAGGTCCTCGCCCACCAGAGGGGCCTGCTGAAGAGTTACCTGGCCTGGTGACCCGGGTGGAGGTCCTCGCCCACCGGAGGGGCGTGCTGAGGGGTGGCCTG"}' | vg view -JaG - > soft.gam
vg convert soft.pg -G soft.gam > soft.gaf
vg view -a soft.gam | jq .sequence > gam.sequence
vg convert soft.pg -F soft.gaf | vg view -a - | jq .sequence > gam2.sequence
diff gam.sequence gam2.sequence
is "$?" 0 "convert gam->gaf->gam on softclipped read preserves sequence"

vg convert soft.pg -F soft.gaf | vg convert soft.pg -G - > soft2.gaf
diff soft.gaf soft2.gaf
is "$?" 0 "convert gam->gaf->gam->gaf makes same gaf each time of soft clipped alignment" 

rm -f soft.pg soft.gam soft.gaf gam.sequence gam2.sequence soft2.gaf

printf "H\tVN:Z:1.0
S\t91194329\tAGGAAGGAGAGGGAG\n" | vg convert -g - -p > floating-ins.pg
printf '{"annotation": {"fragment_length": 1098, "fragment_length_distribution": "-I 542.973684 -D 141.206118", "mapq_applied_cap": 46.364361584299516, "mapq_extended_cap": "Infinity", "mapq_uncapped": 1.5051499783199018, "rescued": true, "secondary_scores": [99.415737573445895]}, "mapping_quality": 1, "name": "ERR903030.51990324", "path": {"mapping": [{"edit": [{"sequence": "GGGCACGGTGGCTCACAGCTGTCACCACNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN", "to_length": 126}], "position": {"is_reverse": true, "node_id": "91194329", "offset": "15"}, "rank": "1"}]}, "quality": "ISEhICEhJCIkJiYdJiUQJSYmHyMmIxAkJiICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC", "sample_name": "HG00514_961a37c", "sequence": "GGGCACGGTGGCTCACAGCTGTCACCACNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN"}' | vg view -JaG - > floating-ins.gam
vg convert floating-ins.pg -G floating-ins.gam > floating-ins.gaf
vg view -a floating-ins.gam | jq .sequence > gam.sequence
vg convert floating-ins.pg -F floating-ins.gaf | vg view -a - | jq .sequence > gam2.sequence
diff gam.sequence gam2.sequence
is "$?" 0 "convert gam->gaf->gam on read with floating insertion preserves sequence"

vg convert floating-ins.pg -F floating-ins.gaf | vg convert floating-ins.pg -G - > floating-ins2.gaf
diff floating-ins.gaf floating-ins2.gaf
is "$?" 0 "convert gam->gaf->gam->gaf makes same gaf each time of floating insertion alignment" 

rm -f floating-ins.pg floating-ins.gam floating-ins.gaf gam.sequence gam2.sequence floating-ins2.gaf

vg convert -g tiny/tiny.gfa | vg view - | sort > tiny.gfa.gfa
vg convert -g tiny/tiny.rgfa | vg view - | sort > tiny.rgfa.gfa
diff tiny.gfa.gfa tiny.rgfa.gfa
is "$?" 0 "converting gfa and equivalent rgfa produces same output"

rm -f tiny.gfa.gfa tiny.rgfa.gfa

is "$(vg convert -g tiny/tiny.rgfa -r 1 | vg view - | grep y | awk '{print $1","$2","$3}')" "P,y[35],2+" "rank-1 rgfa subpath found"

vg convert -g tiny/tiny.rgfa -T gfa-id-mapping.tsv > /dev/null
grep ^S tiny/tiny.rgfa  | awk '{print $2}' | sort > rgfa_nodes
grep ^S tiny/tiny.gfa  | awk '{print $2}' | sort > gfa_nodes
awk '{print $2}' gfa-id-mapping.tsv | sort > rgfa_translated_nodes
awk '{print $3}' gfa-id-mapping.tsv | sort > gfa_translated_nodes
diff rgfa_nodes rgfa_translated_nodes
is "$?" 0 "2nd column of gfa id translation file contains all rgfa nodes"
diff gfa_nodes gfa_translated_nodes
is "$?" 0 "3rd column of gfa id translation file contains all gfa nodes"

rm -f  gfa-id-mapping.tsv rgfa_nodes gfa_nodes rgfa_translated_nodes gfa_translated_nodes

vg convert -g tiny/tiny.gfa -v | vg convert - -f -P x > tiny.gfa.rgfa
is "$(grep ^P tiny.gfa.rgfa | wc -l)" 0 "rgfa output wrote no P-lines"
vg convert -g tiny/tiny.gfa -v | vg convert - -f | sort > tiny.gfa.gfa
vg convert -g tiny.gfa.rgfa -f | sort > tiny.gfa.rgfa.gfa
diff tiny.gfa.gfa tiny.gfa.rgfa.gfa
is "$?" 0 "rgfa export roundtrips back to normal P-lines"

rm -f tiny.gfa.rgfa tiny.gfa.gfa tiny.gfa.rgfa.gfa

# GFA to GBWTGraph to HashGraph to GFA
vg gbwt -o components.gbwt -g components.gg -G graphs/components_walks.gfa
vg convert -b components.gbwt -a components.gg > components.hg
is $? 0 "GBWTGraph to HashGraph conversion"
grep "^S" graphs/components_walks.gfa | sort > sorted.gfa
vg view components.hg | grep "^S" | sort > converted.gfa
cmp sorted.gfa converted.gfa
is $? 0 "GFA -> GBWTGraph -> HashGraph -> GFA conversion maintains segments"

rm -f components.gbwt components.gg
rm -f components.hg
rm -f sorted.gfa converted.gfa<commit_msg>Tests for GBWTGraph to HandleGraph conversion with reference paths<commit_after>#!/usr/bin/env bash

BASH_TAP_ROOT=../deps/bash-tap
. ../deps/bash-tap/bash-tap-bootstrap

PATH=../bin:$PATH # for vg

plan tests 41

vg construct -r complex/c.fa -v complex/c.vcf.gz > c.vg
cat <(vg view c.vg | grep ^S | sort) <(vg view c.vg | grep L | uniq | wc -l) <(vg paths -v c.vg -E) > c.info

vg convert c.vg -x > c.xg
vg convert c.xg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout xg conversion"

rm -f c.xg c1.vg c1.info

vg convert c.vg -a > c.hg
vg convert c.hg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout hash-graph conversion"

rm -f c.hg c1.vg c1.info

vg convert c.vg -p > c.pg
vg convert c.pg -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout packed-graph conversion"

rm -f c.pg c1.vg c1.info

vg convert c.vg -f > c.gfa
vg convert -g c.gfa -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout gfa conversion"

rm -f c.gfa c1.vg c1.info

vg convert c.vg -o > c.odgi
vg convert c.odgi -v > c1.vg
cat <(vg view c1.vg | grep ^S | sort) <(vg view c1.vg | grep L | uniq | wc -l) <(vg paths -v c1.vg -E) > c1.info
diff c.info c1.info
is "$?" 0 "vg convert maintains same nodes throughout ODGI conversion"

rm -f c.vg c.odgi c1.vg c.info c1.info

# some less rigorous tests I made without noticing that the earlier ones had already been written
vg construct -r small/x.fa -v small/x.vcf.gz > x.vg
vg view x.vg > x.gfa

is "$(vg convert -a x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "hash graph conversion looks good"
is "$(vg convert -p x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "packed graph conversion looks good"
is "$(vg convert -v x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "vg conversion looks good"
is "$(vg convert -o x.vg | vg view - | wc -l)" "$(wc -l < x.gfa)" "odgi conversion looks good"
is "$(vg convert -f x.vg | vg convert -g - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa conversion looks good"
is "$(vg convert -x x.vg | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "xg conversion looks good"

is "$(vg convert -g -a x.gfa | vg view - | wc -l)" "$(wc -l < x.gfa)" "on disk gfa conversion looks good"
is "$(cat x.gfa | vg convert -g -a - | vg view - | wc -l)" "$(wc -l < x.gfa)" "streaming gfa conversion looks good"
is "$(vg convert -g -x x.gfa | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa to xg conversion looks good"
is "$(vg convert -g -f x.gfa | vg convert -g - | vg find -n 1 -c 300 -x - | vg view - | wc -l)" "$(wc -l < x.gfa)" "gfa to gfa conversion looks good"

rm x.vg x.gfa
rm -f c.vg c.pg c1.vg c.info c1.info

vg construct -r small/x.fa -v small/x.vcf.gz > x.vg
vg index x.vg -g x.gcsa
vg sim -x x.vg -n 10 -s 23 -a > sim.gam
vg map -x x.vg -g x.gcsa -G sim.gam > sim-rm.gam
vg convert x.vg -G sim-rm.gam -t 1 > sim-rm.gaf
vg convert x.vg -F sim-rm.gaf -t 1 | vg convert x.vg -G - -t 1 > sim-rm2.gaf
diff sim-rm.gaf sim-rm2.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf makes same gaf twice"

vg convert x.vg -G sim-rm.gam | vg convert x.vg -F - | vg convert x.vg -G - | sort > sim-rm2-mt-sort.gaf
sort sim-rm2.gaf > sim-rm2-sort.gaf
diff sim-rm2-sort.gaf sim-rm2-mt-sort.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf gives same result multithreaded as with -t 1"

vg convert x.vg -G sim-rm.gam | bgzip | vg convert x.vg -F - | vg convert x.vg -G - | sort > sim-rm2-mtbg-sort.gaf
diff sim-rm2-sort.gaf sim-rm2-mtbg-sort.gaf
is "$?" 0 "vg convert gam -> gaf.gz -> gam -> gaf gives same result multithreaded as with -t 1"

# some snps and indels
vg map -s "TAATGGATATGTTAAGCTTTTTTTTTCTTTGATTTATTTGAAAAGACGTTTGACAATCTATCGGGTAATGTGGGGAAA" -x x.vg -g x.gcsa > mut.gam
# reverse complement of above
vg map -s "TTTCCCCACATTACCCGATAGATTGTCAAACGTCTTTTCAAATAAATCAAAGAAAAAAAAAGCTTAACATATCCATTA" -x x.vg -g x.gcsa >> mut.gam

vg convert x.vg -G mut.gam -t 1 > mut.gaf
vg convert x.vg -F mut.gaf -t 1 > mut-back.gam
vg view -a mut.gam | jq .path > mut.path
vg view -a mut-back.gam | jq .path > mut-back.path
# Json comparison that is not order dependent: https://stackoverflow.com/a/31933234
is $(jq --argfile a mut.path --argfile b mut-back.path -n 'def post_recurse(f): def r: (f | select(. != null) | r), .; r; def post_recurse: post_recurse(.[]?); ($a | (post_recurse | arrays) |= sort) as $a | ($b | (post_recurse | arrays) |= sort) as $b | $a == $b') true "vg convert gam -> gaf -> gam produces same gam Paths with snps and indels"

vg view -a mut.gam | jq .sequence > mut.seq
vg view -a mut-back.gam | jq .sequence > mut-back.seq
diff mut.seq mut-back.seq
is "$?" 0 "vg convert gam -> gaf -> gam preserves sequence"

vg convert x.vg -G mut-back.gam -t 1 > mut-back.gaf
diff mut.gaf mut-back.gaf
is "$?" 0 "vg convert gam -> gaf -> gam -> gaf makes same gaf twice in presence of indels and snps"

rm -f x.vg x.gcsa sim.gam sim-rm.gam sim-rm.gaf sim-rm2.gaf sim-rm2-mt-sort.gaf sim-rm2-mtbg-sort.gaf sim-rm2-sort.gaf mut.gam mut-back.gam mut.gaf mut-back.gaf mut.path mut-back.path mut.seq mut-back.seq

vg construct -r 1mb1kgp/z.fa -v 1mb1kgp/z.vcf.gz > z.vg 2> /dev/null
vg sim -n 10000 -s 23 -a -x z.vg > sim.gam
vg construct -r 1mb1kgp/z.fa > zflat.vg
vg index zflat.vg -g zflat.gcsa
vg map -x zflat.vg -g zflat.gcsa -G sim.gam > sim-map.gam
vg convert zflat.vg -G sim-map.gam | bgzip > sim-map.gaf.gz
vg convert zflat.vg -F sim-map.gaf.gz > sim-map-back.gam
vg view -a sim-map.gam | jq .sequence | sort > sim-map.sequence
vg view -a sim-map-back.gam | jq .sequence | sort > sim-map-back.sequence
diff sim-map.sequence sim-map-back.sequence
is "$?" 0 "vg convert gam -> gaf -> gam preserves sequences of 1mb1kgp simulated reads"

vg convert zflat.vg -G sim-map-back.gam | sort > sim-map-back.gaf
bgzip -dc sim-map.gaf.gz | sort > sim-map.gaf
diff sim-map-back.gaf sim-map.gaf
is "$?" 0 "vg convert gam -> gaf -> gam ->gaf makes same gaf each time on 1mb1kgp simulated reads"

printf '{"name": "split", "path": {"mapping": [{"edit": [{"from_length": 13, "to_length": 13}], "position": {"node_id": "1", "offset": "10"}}, {"edit": [{"from_length": 2, "to_length": 2}], "position": {"node_id": "3", "offset": "5"}}]}}' | vg view -JaG - > split.gam
vg convert zflat.vg -G split.gam > split.gaf
is "$(awk '{print $13}' split.gaf)" "cs:Z::13-CCAGTGCTC-GCATC:2" "split alignment converted using deletions to represent internal offsets"
vg convert zflat.vg -F split.gaf | vg convert zflat.vg -G - > split-back.gaf
diff split.gaf split-back.gaf
is "$?" 0 "vg convert gam -> gaf ->gam -> gaf makes same gaf each time for split alignment"

rm -f z.vg zflat.vg sim.gam sim-map.gam sim-map-back.gam sim-map.gaf.gz sim-map.sequence sim-map-back.sequence sim-map-back.gaf sim-map.gaf split.gam split.gaf split-back.gaf

printf "H\tVN:Z:1.0
S\t73333\tGGTGGGCGAGGACCTCCACACGTGTCACCA
S\t73368\tGCCCCT
S\t72943\tGGCGACTCTTCAGCAAGCCCCTCCACACGTGT
S\t72940\tC
S\t72941\tGGCCAGGT
S\t73255\tACTCTTCAGCAGGCCCCTCTGGT
S\t72942\tGGGCGAGGACCTCCACACGTGTCACCAGGCCA
S\t73318\tTCAGCA
S\t73367\tA
S\t73271\tA
S\t73289\tC
S\t73317\tC\n" | vg convert -g - -p > soft.pg
printf '{"annotation": {"fragment_length": 242, "fragment_length_distribution": "-I 561.110526 -D 141.152986", "mapq_applied_cap": 23.832780374978935, "mapq_extended_cap": 15, "mapq_uncapped": 10.325140756048304, "secondary_scores": [187.15265582416521, 182.4063994408188]}, "identity": 0.89682539682539686, "mapping_quality": 10, "name": "ERR903030.2067", "path": {"mapping": [{"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "72943", "offset": "26"}}, {"edit": [{"from_length": 32, "to_length": 32}], "position": {"is_reverse": true, "node_id": "72942"}, "rank": "1"}, {"edit": [{"from_length": 23, "to_length": 23}], "position": {"is_reverse": true, "node_id": "73255"}, "rank": "2"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73271"}, "rank": "3"}, {"edit": [{"from_length": 8, "to_length": 8}], "position": {"is_reverse": true, "node_id": "72941"}, "rank": "4"}, {"edit": [{"from_length": 7, "to_length": 7}, {"from_length": 1, "sequence": "C", "to_length": 1}, {"from_length": 2, "to_length": 2}, {"from_length": 1, "sequence": "G", "to_length": 1}, {"from_length": 19, "to_length": 19}], "position": {"is_reverse": true, "node_id": "73333"}, "rank": "5"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "72940"}, "rank": "6"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73289"}, "rank": "7"}, {"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "73368"}, "rank": "8"}, {"edit": [{"from_length": 1, "sequence": "G", "to_length": 1}], "position": {"is_reverse": true, "node_id": "73367"}, "rank": "9"}, {"edit": [{"from_length": 6, "to_length": 6}], "position": {"is_reverse": true, "node_id": "73318"}, "rank": "10"}, {"edit": [{"from_length": 1, "to_length": 1}], "position": {"is_reverse": true, "node_id": "73317"}, "rank": "11"}, {"edit": [{"sequence": "GGGTGGCCTG", "to_length": 10}], "position": {"is_reverse": true, "node_id": "73317", "offset": "1"}, "rank": "12"}]}, "quality": "ISEhISEmJiUmJCYmJSYlJiYmJiYmJCYjJiYdIRsmIyUmJiYmJiElJg4PDx0PDyEkIhsYJCYiECQQJBAcGR0kHx0QGSQPJCImHSElHR0PJA4OGyQbIxwPIg0iHw8PGRwiIR0jAgICAgICAgICAgICAgICAgICAgICAgICAgIC", "sample_name": "HG00514_961a37c_gssw", "score": 106, "sequence": "GTCGCCTGGCCTGGTGACACGTGTGGAGGTCCTCGCCCACCAGAGGGGCCTGCTGAAGAGTTACCTGGCCTGGTGACCCGGGTGGAGGTCCTCGCCCACCGGAGGGGCGTGCTGAGGGGTGGCCTG"}' | vg view -JaG - > soft.gam
vg convert soft.pg -G soft.gam > soft.gaf
vg view -a soft.gam | jq .sequence > gam.sequence
vg convert soft.pg -F soft.gaf | vg view -a - | jq .sequence > gam2.sequence
diff gam.sequence gam2.sequence
is "$?" 0 "convert gam->gaf->gam on softclipped read preserves sequence"

vg convert soft.pg -F soft.gaf | vg convert soft.pg -G - > soft2.gaf
diff soft.gaf soft2.gaf
is "$?" 0 "convert gam->gaf->gam->gaf makes same gaf each time of soft clipped alignment" 

rm -f soft.pg soft.gam soft.gaf gam.sequence gam2.sequence soft2.gaf

printf "H\tVN:Z:1.0
S\t91194329\tAGGAAGGAGAGGGAG\n" | vg convert -g - -p > floating-ins.pg
printf '{"annotation": {"fragment_length": 1098, "fragment_length_distribution": "-I 542.973684 -D 141.206118", "mapq_applied_cap": 46.364361584299516, "mapq_extended_cap": "Infinity", "mapq_uncapped": 1.5051499783199018, "rescued": true, "secondary_scores": [99.415737573445895]}, "mapping_quality": 1, "name": "ERR903030.51990324", "path": {"mapping": [{"edit": [{"sequence": "GGGCACGGTGGCTCACAGCTGTCACCACNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN", "to_length": 126}], "position": {"is_reverse": true, "node_id": "91194329", "offset": "15"}, "rank": "1"}]}, "quality": "ISEhICEhJCIkJiYdJiUQJSYmHyMmIxAkJiICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIC", "sample_name": "HG00514_961a37c", "sequence": "GGGCACGGTGGCTCACAGCTGTCACCACNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN"}' | vg view -JaG - > floating-ins.gam
vg convert floating-ins.pg -G floating-ins.gam > floating-ins.gaf
vg view -a floating-ins.gam | jq .sequence > gam.sequence
vg convert floating-ins.pg -F floating-ins.gaf | vg view -a - | jq .sequence > gam2.sequence
diff gam.sequence gam2.sequence
is "$?" 0 "convert gam->gaf->gam on read with floating insertion preserves sequence"

vg convert floating-ins.pg -F floating-ins.gaf | vg convert floating-ins.pg -G - > floating-ins2.gaf
diff floating-ins.gaf floating-ins2.gaf
is "$?" 0 "convert gam->gaf->gam->gaf makes same gaf each time of floating insertion alignment" 

rm -f floating-ins.pg floating-ins.gam floating-ins.gaf gam.sequence gam2.sequence floating-ins2.gaf

vg convert -g tiny/tiny.gfa | vg view - | sort > tiny.gfa.gfa
vg convert -g tiny/tiny.rgfa | vg view - | sort > tiny.rgfa.gfa
diff tiny.gfa.gfa tiny.rgfa.gfa
is "$?" 0 "converting gfa and equivalent rgfa produces same output"

rm -f tiny.gfa.gfa tiny.rgfa.gfa

is "$(vg convert -g tiny/tiny.rgfa -r 1 | vg view - | grep y | awk '{print $1","$2","$3}')" "P,y[35],2+" "rank-1 rgfa subpath found"

vg convert -g tiny/tiny.rgfa -T gfa-id-mapping.tsv > /dev/null
grep ^S tiny/tiny.rgfa  | awk '{print $2}' | sort > rgfa_nodes
grep ^S tiny/tiny.gfa  | awk '{print $2}' | sort > gfa_nodes
awk '{print $2}' gfa-id-mapping.tsv | sort > rgfa_translated_nodes
awk '{print $3}' gfa-id-mapping.tsv | sort > gfa_translated_nodes
diff rgfa_nodes rgfa_translated_nodes
is "$?" 0 "2nd column of gfa id translation file contains all rgfa nodes"
diff gfa_nodes gfa_translated_nodes
is "$?" 0 "3rd column of gfa id translation file contains all gfa nodes"

rm -f  gfa-id-mapping.tsv rgfa_nodes gfa_nodes rgfa_translated_nodes gfa_translated_nodes

vg convert -g tiny/tiny.gfa -v | vg convert - -f -P x > tiny.gfa.rgfa
is "$(grep ^P tiny.gfa.rgfa | wc -l)" 0 "rgfa output wrote no P-lines"
vg convert -g tiny/tiny.gfa -v | vg convert - -f | sort > tiny.gfa.gfa
vg convert -g tiny.gfa.rgfa -f | sort > tiny.gfa.rgfa.gfa
diff tiny.gfa.gfa tiny.gfa.rgfa.gfa
is "$?" 0 "rgfa export roundtrips back to normal P-lines"

rm -f tiny.gfa.rgfa tiny.gfa.gfa tiny.gfa.rgfa.gfa


# GFA to GBWTGraph to HashGraph to GFA
vg gbwt -o components.gbwt -g components.gg -G graphs/components_walks.gfa
vg convert -b components.gbwt -a components.gg > components.hg 2> /dev/null
is $? 0 "GBWTGraph to HashGraph conversion"
grep "^S" graphs/components_walks.gfa | sort > sorted.gfa
vg view components.hg | grep "^S" | sort > converted.gfa
cmp sorted.gfa converted.gfa
is $? 0 "GFA -> GBWTGraph -> HashGraph -> GFA conversion maintains segments"

rm -f components.gbwt components.gg
rm -f components.hg
rm -f sorted.gfa converted.gfa


# GFA to GBWTGraph with paths and walks
vg gbwt -o components.gbwt -g components.gg -G graphs/components_walks_paths.gfa
vg convert -g -a graphs/components_walks_paths.gfa > direct.hg
vg paths -v direct.hg -A > correct_paths.gaf

# GBWTGraph to HashGraph with paths and walks
vg convert -b components.gbwt -a components.gg > components.hg
is $? 0 "GBWTGraph to HashGraph conversion with reference paths"
vg paths -A -v components.hg > hg_paths.gaf
cmp hg_paths.gaf correct_paths.gaf
is $? 0 "GBWTGraph to HashGraph conversion creates the correct reference paths"

# GBWTGraph to XG with paths and walks
vg convert -b components.gbwt -x components.gg > components.xg
is $? 0 "GBWTGraph to XG conversion with reference paths"
vg paths -A -v components.xg > xg_paths.gaf
cmp xg_paths.gaf correct_paths.gaf
is $? 0 "GBWTGraph to XG conversion creates the correct reference paths"

rm -f components.gbwt components.gg
rm -f direct.hg correct_paths.gaf
rm -f components.hg hg_paths.gaf
rm -f components.xg xg_paths.gaf
{'fraction_non_alphanumeric': 0.16629756392188444, 'fraction_numerical': 0.04197704852023354, 'mean_word_length': 3.861034089055619, 'pattern_counts': {'":': 314, '<': 59, '<?xml version=': 0, '>': 236, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 8, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '775300', 'n_tokens_mistral': 13587, 'n_tokens_neox': 12088, 'n_words': 3622}