diff --git a/src/minimizer_mapper.cpp b/src/minimizer_mapper.cpp index 56c663933f..c12367ea7c 100644 --- a/src/minimizer_mapper.cpp +++ b/src/minimizer_mapper.cpp @@ -2803,14 +2803,15 @@ pair, vector> MinimizerMapper::map_paired(Alignment } } - //Annotate top pair with its fragment distance, properly-paired-ness, fragment length distrubution, and secondary scores + // Annotate top pair with its fragment distance, properly-paired-ness, fragment length distrubution, and secondary scores bool properly_paired = distances.front() == std::numeric_limits::max() ? false : (std::abs(distances.front()-fragment_length_distr.mean()) <= 6.0*fragment_length_distr.std_dev()) ; string distribution = "-I " + to_string(fragment_length_distr.mean()) + " -D " + to_string(fragment_length_distr.std_dev()); for (auto r : {0, 1}) { set_annotation(mappings[r].front(), "fragment_length", distance_to_annotation(distances.front())); set_annotation(mappings[r].front(), "proper_pair", properly_paired); - set_annotation(mappings[r].front(),"fragment_length_distribution", distribution); + // GAM actually has a fragment_length_distribution field. + mappings[r].front().set_fragment_length_distribution(distribution); set_annotation(mappings[r].front(),"secondary_scores", scores); } } diff --git a/src/subcommand/surject_main.cpp b/src/subcommand/surject_main.cpp index 54adfb0c17..209185d5d4 100644 --- a/src/subcommand/surject_main.cpp +++ b/src/subcommand/surject_main.cpp @@ -46,8 +46,6 @@ void help_surject(char** argv) { << " -F, --into-paths FILE surject into path names listed in" << endl << " HTSlib sequence dictionary or path list FILE" << endl << " -n, --into-ref NAME surject into this reference assembly" << endl - << " -i, --interleaved GAM is interleaved paired-ended, so pair reads" << endl - << " when outputting HTS formats" << endl << " -M, --multimap include secondary alignments to all" << endl << " overlapping paths instead of just primary" << endl << " -G, --gaf-input input file is GAF instead of GAM" << endl @@ -55,6 +53,12 @@ void help_surject(char** argv) { << " -c, --cram-output write CRAM instead of GAM to stdout" << endl << " -b, --bam-output write BAM instead of GAM to stdout" << endl << " -s, --sam-output write SAM instead of GAM to stdout" << endl + << " -i, --interleaved input is interleaved paired-ended, so pair reads" << endl + << " when outputting HTS formats (SAM/BAM/CRAM)" << endl + << " -U, --force-unpaired surject reads as unpaired even if they appear paired" << endl + << " -f, --max-frag-len N reads with fragment lengths greater than N won't be" << endl + << " marked properly paired in HTS formats." << endl + << " 0 for unlimited. (default: unlimited)" << endl << " -u, --supplementary divide into supplementary alignments as necessary" << endl << " -B, --left-align attempt to left-align indels" << endl << " -l, --subpath-local let the multipath mapping surjection produce local" << endl @@ -78,8 +82,6 @@ void help_surject(char** argv) { << " of the 10x-scaled scoring parameters" << endl << " -N, --sample NAME set this sample name for all reads" << endl << " -R, --read-group NAME set this read group for all reads" << endl - << " -f, --max-frag-len N reads with fragment lengths greater than N won't be" << endl - << " marked properly paired in SAM/BAM/CRAM" << endl << " -L, --list-all-paths annotate SAM records with a list of all attempted" << endl << " re-alignments to paths in SS tag" << endl << " -H, --graph-aln annotate SAM records with cs-style difference string" << endl @@ -136,10 +138,48 @@ static void ensure_alignment_is_for_graph(const Logger& logger, const MultipathA } } +/// Returns true if the given alignment seems like it was aligned paired-end, +/// and false otherwise. Should not detect paired-end reads that were aligned +/// single-ended. +static bool smells_paired(const Alignment& aln) { + // The vg way to find this out is the fragment_prev/fragment_next fields. + // vg GAF has tags for them + return aln.has_fragment_prev() || + aln.has_fragment_next() || + // If an alignment just came in from BAM somehow, this would be set. + aln.read_paired() || + // vg GAF will record a tag that carries this. + // Maybe some other GAF writer supports this but not next/prev? + has_annotation(aln, "proper_pair") || + // Supplementary alignments might have this instead. + has_annotation(aln, "mate_info"); +} + +/// Returns true if the given multipath alignment seems like it was aligned +/// paired-end, and false otherwise. Should not detect paired-end reads that +/// were aligned single-ended. +static bool smells_paired(const MultipathAlignment& mpaln) { + return has_annotation(mpaln, "fragment_length_distribution") || + has_annotation(mpaln, "proper_pair") || + has_annotation(mpaln, "mate_info"); +} + +/// Helper for producing an error message about a read that smells paired when +/// running in single-ended mode. +static void smells_paired_error(const Logger& logger, const string& name) { + #pragma omp critical (cerr) + logger.error() << "Read " << name << " seems to have been aligned as a " + << "paired-end read. Provide either the -i/--interleaved " + << "option (and, ideally, a value for -f/--max-frag-len) " + << "to surject as paired-end, or -U/--force-unpaired " + << "to surject as single-ended anyway." << endl; +} + /// This is a helper for the cases where we find two alignments /// that are adjacent in the input but not paired. static void adjacent_but_not_paired_error(const Logger& logger, const string& name1, const string& name2) { - + + #pragma omp critical (cerr) logger.error() << "alignments " << name1 << " and " << name2 << " are adjacent but not paired" << endl; } @@ -163,9 +203,10 @@ int main_surject(int argc, char** argv) { string input_format = "GAM"; bool spliced = false; bool interleaved = false; + bool force_unpaired = false; + std::optional max_frag_len; string sample_name; string read_group; - int32_t max_frag_len = 0; int compress_level = 9; int64_t min_splice_length = 20; size_t watchdog_timeout = 10; @@ -204,13 +245,15 @@ int main_surject(int argc, char** argv) { {"subpath-local", no_argument, 0, 'l'}, {"max-tail-len", required_argument, 0, 'T'}, {"max-graph-scale", required_argument, 0, 'g'}, - {"interleaved", no_argument, 0, 'i'}, {"multimap", no_argument, 0, 'M'}, {"gaf-input", no_argument, 0, 'G'}, {"gamp-input", no_argument, 0, 'm'}, {"cram-output", no_argument, 0, 'c'}, {"bam-output", no_argument, 0, 'b'}, {"sam-output", no_argument, 0, 's'}, + {"interleaved", no_argument, 0, 'i'}, + {"force-unpaired", no_argument, 0, 'U'}, + {"max-frag-len", required_argument, 0, 'f'}, {"supplementary", no_argument, 0, 'u'}, {"off-ref-position", no_argument, 0, OPT_OFF_REF_POS}, {"left-align", no_argument, 0, 'B'}, @@ -224,7 +267,6 @@ int main_surject(int argc, char** argv) { {"extra-gap-cost", required_argument, 0, 'E'}, {"sample", required_argument, 0, 'N'}, {"read-group", required_argument, 0, 'R'}, - {"max-frag-len", required_argument, 0, 'f'}, {"list-all-paths", no_argument, 0, 'L'}, {"graph-aln", no_argument, 0, 'H'}, {"compression", required_argument, 0, 'C'}, @@ -235,7 +277,7 @@ int main_surject(int argc, char** argv) { }; int option_index = 0; - c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:iGmcbsuBN:R:f:C:t:D:SPI:a:AE:LHMVw:r", + c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:GmcbsiUf:uBN:R:C:t:D:SPI:a:AE:LHMVw:r", long_options, &option_index); // Detect the end of the options. @@ -273,10 +315,6 @@ int main_surject(int argc, char** argv) { max_graph_scale.reset(new double(parse(optarg))); break; - case 'i': - interleaved = true; - break; - case 'M': multimap = true; break; @@ -302,6 +340,18 @@ int main_surject(int argc, char** argv) { output_format = "SAM"; break; + case 'i': + interleaved = true; + break; + + case 'U': + force_unpaired = true; + break; + + case 'f': + max_frag_len = parse(optarg); + break; + case 'u': report_supplementary = true; break; @@ -350,10 +400,6 @@ int main_surject(int argc, char** argv) { read_group = optarg; break; - case 'f': - max_frag_len = parse(optarg); - break; - case 'C': compress_level = parse(optarg); break; @@ -407,6 +453,21 @@ int main_surject(int argc, char** argv) { prune_anchors = (read_length == "long"); } + // Validate configuration + if (!interleaved && max_frag_len.has_value()) { + logger.error() << "-f/--max-frag-len can only be used with paired-end reads, " + << "but -i/--interleaved was not provided." << endl; + } + if (interleaved && !max_frag_len.has_value()) { + // TODO: Once we get fragment distribution learning, this will be less of a problem. + logger.warn() << "Running in paired-end mode without -f/--max-frag-len. " + << "Reads will be assumed to be properly paired at any distance!" << endl; + } + if (interleaved && force_unpaired) { + logger.error() << "Cannot use -U/--force-unpaired and -i/--interleaved at the same time. " + << "To surject interleaved reads as unpaired, use just -U/--force-unpaired." << endl; + } + string file_name = get_input_file_name(optind, argc, argv); if (have_input_file(optind, argc, argv)) { @@ -562,7 +623,6 @@ int main_surject(int argc, char** argv) { !src2.has_fragment_prev() || src2.fragment_prev().name() != src1.name()) { -#pragma omp critical (cerr) adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } @@ -572,13 +632,11 @@ int main_surject(int argc, char** argv) { !src1.has_fragment_prev() || src1.fragment_prev().name() != src2.name()) { -#pragma omp critical (cerr) adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } } else { // Alignments aren't paired up at all -#pragma omp critical (cerr) adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } @@ -616,7 +674,7 @@ int main_surject(int argc, char** argv) { auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse())); if (!is_supplementary(surjected1[i]) && it != strand_idx2.end()) { // the alignments are paired on this strand - alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len); + alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len.value_or(0)); } else { // supplementary or unpaired @@ -695,6 +753,10 @@ int main_surject(int argc, char** argv) { ensure_alignment_is_for_graph(logger, src, *xgidx); } + if (!force_unpaired && smells_paired(src)) { + smells_paired_error(logger, src.name()); + } + // Preprocess read to set metadata before surjection set_metadata(src); @@ -750,13 +812,11 @@ int main_surject(int argc, char** argv) { // TODO: Integrate into for_each_interleaved_pair_parallel when running on Alignments. if (src1.paired_read_name() != src2.name() || src2.paired_read_name() != src1.name()) { -#pragma omp critical (cerr) adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } else if (src1.paired_read_name().empty() || src2.paired_read_name().empty()) { // Alignments aren't paired up at all -#pragma omp critical (cerr) adjacent_but_not_paired_error(logger, src1.name(), src2.name()); } @@ -861,7 +921,7 @@ int main_surject(int argc, char** argv) { } // write to output - vector tlen_limits(surjected.size(), max_frag_len); + vector tlen_limits(surjected.size(), max_frag_len.value_or(0)); mp_alignment_emitter.emit_pairs(src1.name(), src2.name(), std::move(surjected), &positions, &tlen_limits); mp_alignment_emitter.emit_paired_independent(src1.name(), src2.name(), std::move(surjected_unpaired1), std::move(surjected_unpaired2), @@ -891,6 +951,10 @@ int main_surject(int argc, char** argv) { ensure_alignment_is_for_graph(logger, src, *xgidx); } + if (!force_unpaired && smells_paired(src)) { + smells_paired_error(logger, src.name()); + } + multipath_alignment_t mp_src; from_proto_multipath_alignment(src, mp_src); diff --git a/src/surjector.cpp b/src/surjector.cpp index dc19d9f76f..d89224f9d6 100644 --- a/src/surjector.cpp +++ b/src/surjector.cpp @@ -38,7 +38,7 @@ using namespace std; Surjector::Surjector(const PathPositionHandleGraph* graph) : graph(graph), choose_band_padding(algorithms::pad_band_constant(1)) { if (!graph) { - cerr << "error:[Surjector] Failed to provide an graph to the Surjector" << endl; + cerr << "error:[Surjector] Failed to provide a graph to the Surjector" << endl; exit(1); } diff --git a/test/t/15_vg_surject.t b/test/t/15_vg_surject.t index 6d25d8cab3..0c4819567b 100644 --- a/test/t/15_vg_surject.t +++ b/test/t/15_vg_surject.t @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap PATH=../bin:$PATH # for vg -plan tests 78 +plan tests 81 vg construct -r small/x.fa >j.vg vg index -x j.xg j.vg @@ -192,11 +192,11 @@ is "$(vg surject -x j.vg -b --graph-aln r.gam | samtools view | grep 'GR:Z:' | w rm -f h.vg h.gcsa r.gam r.sam x.sub.fa j.vg j.gcsa j.gcsa.lcp j.sub.vg j.sub.gcsa j.sub.gcsa.lcp r.sub.gam r.sub.sam r.sub.sam path_info.tsv -vg surject -s -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam +vg surject -U -s -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam is "$?" 0 "vg surject does not crash when surjecting a read that grazes the reference with a deletion" is "$(cat perpendicular.sam | grep -v "^@" | cut -f2)" "4" "vg surject leaves a read that grazes the reference with a deletion unmapped" -vg surject -s --prune-low-cplx -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam +vg surject -U -s --prune-low-cplx -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam is "$?" 0 "vg surject does not crash when surjecting a read that grazes the reference with a deletion and pruning low complexity anchors" is "$(cat perpendicular.sam | grep -v "^@" | cut -f2)" "4" "vg surject leaves a read that grazes the reference with a deletion unmapped when pruning low complexity anchors" @@ -209,12 +209,22 @@ vg sim -x x.xg -n 20 -l 40 -p 60 -v 10 -a --random-seed 123 > x.gam vg mpmap -x x.xg -g x.gcsa -n dna --suppress-mismapping -B -G x.gam -i -F GAM -I 60 -D 10 -t 1 > mapped.gam vg mpmap -x x.xg -g x.gcsa -n dna --suppress-mismapping -B -G x.gam -i -F GAMP -I 60 -D 10 -t 1 > mapped.gamp -is "$(vg surject -x x.xg -s -t 1 mapped.gam | grep -v '@' | wc -l)" 40 "GAM surject can return only primaries" -is "$(vg surject -x x.xg -M -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return multimappings" +vg surject -x x.xg -s -t 1 mapped.gam >/dev/null 2>/dev/null +is "$?" "1" "GAM surject produces an error when given paired reads when not in interleaved mode" + +vg surject -x x.xg -m -s -t 1 mapped.gamp >/dev/null 2>/dev/null +is "$?" "1" "GAMP surject produces an error when given paired reads when not in interleaved mode" + +vg convert x.xg -G mapped.gam -t 1 >mapped.gaf +vg surject -x x.xg -G -s -t 1 mapped.gaf >/dev/null 2>/dev/null +is "$?" 1 "GAF surject produces an error when given paired reads when not in interleaved mode" + +is "$(vg surject -x x.xg -U -s -t 1 mapped.gam | grep -v '@' | wc -l)" 40 "GAM surject can return only primaries" +is "$(vg surject -x x.xg -M -U -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return multimappings" is "$(vg surject -x x.xg -M -i -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return paired multimappings" -is "$(vg surject -x x.xg -s -m -t 1 mapped.gamp | grep -v '@' | wc -l)" 40 "GAMP surject can return only primaries" -is "$(vg surject -x x.xg -M -m -s -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings" -is "$(vg surject -x x.xg -M -m -s -i -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings" +is "$(vg surject -x x.xg -U -s -m -t 1 mapped.gamp | grep -v '@' | wc -l)" 40 "GAMP surject can return only primaries" +is "$(vg surject -x x.xg -M -U -m -s -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings" +is "$(vg surject -x x.xg -M -i -m -s -i -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return paired multimappings" vg construct -r tiny/tiny.fa > tiny.vg vg surject -x tiny.vg -s -t 1 mapped.gam >/dev/null 2>err.txt diff --git a/test/t/48_vg_convert.t b/test/t/48_vg_convert.t index 67f006aef0..ae46da4dd4 100644 --- a/test/t/48_vg_convert.t +++ b/test/t/48_vg_convert.t @@ -156,7 +156,7 @@ 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 +printf '{"annotation": {"fragment_length": 242, "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 @@ -171,7 +171,7 @@ 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 +printf '{"annotation": {"fragment_length": 1098, "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