<p><strong>Edit</strong>: I've written a short blog post describing my use of nodeSeqs.sh to investigate the functions of sequences at the nodes of short-read assembly graphs for a set of human oral microbiome samples.</p>
<p><a href="https://albertvill.com/posts/nodeseqs_oral/" rel="nofollow noreferrer">https://albertvill.com/posts/nodeseqs_oral/</a></p>
<p><a href="https://www.merriam-webster.com/dictionary/TL%3BDR" rel="nofollow noreferrer">TL;DR</a> -- transposases and recombinases are enriched, as expected.</p>
<hr />
<p>Based on <a href="https://bioinformatics.stackexchange.com/a/18528/3967">Maximilian's excellent answer</a>, I have created a script to address my own question.</p>
<p><a href="https://github.com/acvill/nodeSeqs" rel="nofollow noreferrer">https://github.com/acvill/nodeSeqs</a></p>
<p>Here are the details on how it works:</p>
<blockquote>
<p>Link lines (beginning with <code>L</code>) are extracted from the GFA file, and segment occurrences are counted from both the <code>From</code> and <code>To</code> fields. To ensure that the final sequences capture the overlapping (linked) portions of segments, links are processed as segment-orientation pairs. Segment IDs with at least the specified minimum degree (<code>-d</code>) are written to <code>temp_XXXXX/segments.txt</code>, and the complementary segment lines (beginning with <code>S</code>) are pulled. <em>k</em>-mer counts (<code>KC</code> tag) are used to compute coverage for each segment using the approximation from <a href="https://dx.doi.org/10.7717%2Fpeerj.2681" rel="nofollow noreferrer">Gonnella &amp; Kurtz 2016</a>; rewritten:</p>
<p><a href="https://i.sstatic.net/hjwHC.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/hjwHC.png" alt="coverage estimation equation" /></a></p>
<p>where <em>K</em><sub><em>S</em></sub> is the <em>k</em>-mer count for a segment, <em>L</em><sub><em>S</em></sub> is the length of a segment, and <em>k</em> is the <em>k</em>-mer size used in de Bruijn graph construction (<code>-k</code>). For (meta)SPAdes, graphs are built iteratively using an increasing <em>k</em>-mer size, so <em>k</em> is equal to the largest <em>k</em>-mer size used. Sequences are extracted corresponding to segments whose coverage is at least the specified minimum coverage (<code>-m</code>). Sequences longer than the specified context (<code>-c</code>) are truncated to the context length with respect to the link orientation. Including the <code>-t</code> flag will automatically remove the temporary folder when the program is finished, which holds intermediate files and <code>log.txt</code>.</p>
</blockquote>
<p>The code is indeed &quot;shell stuff&quot;, though hopefully it is not too &quot;kludgey&quot;. ;-)</p>

CC BY-SA 4.0
2022-02-14T17:41:32.047Z
+1 for putting in the k-mer coverage approximation, which I neglected. Not too kludgey, unlike my answer! I think that you could accept your own answer if you liked, it is a better and more useful approach than mine, and fully answers your question.
Maximilian Press
CC BY-SA 4.0
2022-02-14T18:04:29.367Z
@MaximilianPress thanks for the feedback. Since your answer provided the framework to create my answer, I'm happy to keep the accept in place.
acvill
CC BY-SA 4.0
2022-02-12T03:52:18.133Z
2022-03-28T18:53:55.063Z
2022-03-28T18:53:55.063Z
acvill
<p>I don't have a recommendation of specific tools, but we might be able to use the GFA file itself to get at this with just bash commands.</p>
<p>For these examples I use the linked GFA <a href="https://towardsdatascience.com/visualising-assembly-graphs-fb631f46bbd1" rel="noreferrer">here</a>, which is a tutorial with more useful info.</p>
<p>GFA segment (node) lines start with <code>S</code>, link lines start with <code>L</code> (spec <a href="http://gfa-spec.github.io/GFA-spec/" rel="noreferrer">here</a>):</p>
<pre><code>S   3438628 CACAATCTAAGAAATCTGAACTACCTGAAACAGGTGGAGAAGAATCAACAAACAACGGCATGTTGTTCGGCGGATTATTTAGCATTTTAG...
L   3461020 +   3461022 +   55M
</code></pre>
<p>Links connect the numeric IDs of segments to each other.</p>
<h2>Find high-degree nodes</h2>
<p>We can just count the occurrences of node IDs in links, which should give us the degree.</p>
<pre><code>egrep &quot;\bL\b&quot; | cut -f4 assembly_graph_with_scaffolds.gfa  | sort | uniq -c | sort | tail

# prints:
      
      5 3461076
      5 3461246
      5 3461286
      5 591392
      5 606626
      5 634272
      5 721568
      5 766622
      5 822586
      6 3455660
</code></pre>
<p>So we can see that there are a lot of links of segment <code>3455660</code>. We can confirm this in the GFA:</p>
<pre><code>grep 3455660 assembly_graph_with_scaffolds.gfa
# prints:
S       3455660 CTGATAAACGGCGGGAACAGAACCAACACTACGCGTTGCTCCGATCTCAACAACTTCTAAG   KC:i:2936
L       3455660 +       56698   -       55M
L       3455660 +       3455666 +       55M
L       3455660 +       3456290 -       55M
L       3455660 -       3455652 -       55M
L       3455660 -       3456398 -       55M
P       NODE_6_length_300278_cov_129.624456_1   3448698+,3456292-,3448662+,3448670+,3451896-,3460544-,3460388+,3434104-,3434124-,56698+,3455660-,3455652-,3460374+,3460806-,3460824-,3460806-,3460824-,3460806-,3460530-,3449022-,3454210+,3459744-,3459742-    *
P       NODE_8_length_282610_cov_136.723533_1   3460440-,3434818+,3276488+,3461076-,3461070-,3461062-,3461054-,3461046-,3461038-,3460426-,3461020+,3461022+,3461020+,3460516-,3455644+,3455652+,3455660+,3456290-,3448662+,3460546-,3275518-,3387950+,3456304+,3367628+,523826-,3418048-,3456388+,3460376-,3327686+,3455382-,3460586+,655086-,3275436+,3460670+,3278936+,3460218+,3460346+,3460498+,3459664+,3460496+,3459664+,3459200+,3459424-,3458356+,3455426-,3455420-     *
P       NODE_39_length_23557_cov_173.180283_1   3456398+,3455660+,3455666+,3418048-,3273044-,3274776+,3433484-,3433478-,3433470-,3433462-,3456812-,3460134+,3455826-,3319248+,3454700+,3421874+,3453038+,3459348-,3453046-,3234810-,3460596+,3460980+,3460654+  *
</code></pre>
<h2>Extract sequences proximal to a node</h2>
<p>Having extracted the high-degree node, we can use the output of <code>grep</code> above to find links to other sequences:</p>
<pre><code>grep 3455660 assembly_graph_with_scaffolds.gfa | egrep &quot;\bL&quot;
# prints:
L   3455660 +   56698   -   55M
L   3455660 +   3455666 +   55M
L   3455660 +   3456290 -   55M
L   3455660 -   3455652 -   55M
L   3455660 -   3456398 -   55M
</code></pre>
<p>We can then go lookup <code>S</code> lines for those segments:</p>
<pre><code>egrep &quot;\bS\W56698&quot; assembly_graph_with_scaffolds.gfa
# prints:
S   56698   GCATGTCGAGAAAACACCTTAGAAGTTGTTGAGATCGGAGCAACGCGTAGTGTTGGTTCTGTTCCCGCCGTT    KC:i:2067
</code></pre>
<h2>Coverage filtering</h2>
<p>Read count (<code>RC</code> tag) is an optional field for GFA, which does not exist in the linked file.</p>
<p>However, we do have k-mer count (<code>KC</code>), which you could use. Presumably high k-mer count would correlate well to coverage. I don't have a suggestion of a nifty tool that filters on this, I'd use a custom script.</p>
<h1>Summary</h1>
<p>Probably there exists some command line tool that does this stuff, but I don't know about it. It has to all exist inside assemblers, but it might be pretty hard to pull out.</p>
<p>If I had a little more time, I might write a simple tool that uses knowledge of the GFA structure to perform these kinds of data slicing operations, it seems pretty useful. It really shouldn't be hard to do this stuff!</p>
<p>A related approach, that I believe has a somewhat different focus, can be found <a href="https://www.biorxiv.org/content/10.1101/2020.05.21.109587v2.full.pdf" rel="noreferrer">here</a>.</p>

CC BY-SA 4.0
2022-02-08T21:55:16.040Z
Thanks Max! This is exactly the kind of insight I am looking for. I read the GFA specification but was honestly confused by it. This makes everything clear. Excellent answer
acvill
CC BY-SA 4.0
2022-02-09T20:11:16.013Z
In my searching, I've also found [RGFATools](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5103826/), which provides a framework for parsing and editing GFA files. They include a method for computing the coverage of a segment given the length of a segment, the *k*-mer count of the segment, and the average *k*-mer count of all segments. I think I'll try my hand at implementing this for filtering
acvill
CC BY-SA 4.0
2022-02-09T20:13:51.103Z
@acvill that is interesting- I knew about gfapy (https://academic.oup.com/bioinformatics/article/33/19/3094/3883940) though I have never used it. If you find a better way to do it with that rather than all the kludgey shell stuff, consider posting your own answer.
Maximilian Press
CC BY-SA 4.0
2022-02-08T21:29:34.880Z
2022-02-08T21:35:41.433Z
2022-02-08T21:35:41.433Z
Maximilian Press
<p>I often use <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411777/" rel="noreferrer">metaSPAdes</a> to assemble short reads from human microbiomes. My simplified understanding of short-read de Bruijn graph assemblers is that they fail where ambiguous paths cannot be resolved. While it can be said that these points of failure may be due to certain <strong>structural</strong> features of the underlying sequences -- long repetitive sequences or highly conserved sequences found at multiple loci -- I'm interested in the <strong>functional</strong> classification of sequences that are proximal to nodes on assembly graphs. Put another way, what kinds of genes are poorly represented in metagenomes because they are poorly assembled? Naively, I expect mobile elements like transposons and integrons to be in such a list of genes, since they are often flanked by repeated sequences and can be present in many copies across multiple genomic contexts, though I want to take an empirical approach.</p>
<p>I've been dabbling with <a href="http://rrwick.github.io/Bandage/" rel="noreferrer">Bandage</a>, which allows you navigate assembly graphs in a GUI, though I'm having trouble extracting features of the underlying network that Bandage builds from the input (.fastg or .gfa). I'd like to know if there are any tools or workflows that would allow me to, given an assembly graph in <a href="http://gfa-spec.github.io/GFA-spec/GFA1.html" rel="noreferrer">GFA format</a> ...</p>
<ol>
<li>filter nodes by degree,</li>
<li>extract sequences proximal to high-degree nodes, and</li>
<li>filter pulled sequences based on sequencing depth, so as to preclude analyzing node-adjacent sequences in the case where the assembler failed due to low depth.</li>
</ol>
<p>I'm open to the possibility that I'm taking the wrong approach to this problem, so feel free to answer with a frame challenge.</p>

<assembly><metagenome><graphs><de-bruijn-graphs>
Extract sequence context of high-degree nodes in assembly graphs
# Extract sequence context of high-degree nodes in assembly graphs

(Asked by: username_0 on 2022-02-02)

I often use [metaSPAdes](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411777/) to assemble short reads from human microbiomes. My simplified understanding of short-read de Bruijn graph assemblers is that they fail where ambiguous paths cannot be resolved. While it can be said that these points of failure may be due to certain **structural** features of the underlying sequences -- long repetitive sequences or highly conserved sequences found at multiple loci -- I'm interested in the **functional** classification of sequences that are proximal to nodes on assembly graphs. Put another way, what kinds of genes are poorly represented in metagenomes because they are poorly assembled? Naively, I expect mobile elements like transposons and integrons to be in such a list of genes, since they are often flanked by repeated sequences and can be present in many copies across multiple genomic contexts, though I want to take an empirical approach.

I've been dabbling with [Bandage](http://rrwick.github.io/Bandage/), which allows you navigate assembly graphs in a GUI, though I'm having trouble extracting features of the underlying network that Bandage builds from the input (.fastg or .gfa). I'd like to know if there are any tools or workflows that would allow me to, given an assembly graph in [GFA format](http://gfa-spec.github.io/GFA-spec/GFA1.html) ...

1. filter nodes by degree,
2. extract sequences proximal to high-degree nodes, and
3. filter pulled sequences based on sequencing depth, so as to preclude analyzing node-adjacent sequences in the case where the assembler failed due to low depth.

I'm open to the possibility that I'm taking the wrong approach to this problem, so feel free to answer with a frame challenge.

**Tags:** assembly, metagenome, graphs, de-bruijn-graphs

---

## 2 answers

### Answer by username_1 on 2022-02-08. Score: 4 | Accepted answer

I don't have a recommendation of specific tools, but we might be able to use the GFA file itself to get at this with just bash commands.
For these examples I use the linked GFA here, which is a tutorial with more useful info.
GFA segment (node) lines start with S, link lines start with L (spec here):
S   3438628 CACAATCTAAGAAATCTGAACTACCTGAAACAGGTGGAGAAGAATCAACAAACAACGGCATGTTGTTCGGCGGATTATTTAGCATTTTAG...
L   3461020 +   3461022 +   55M

Links connect the numeric IDs of segments to each other.
Find high-degree nodes
We can just count the occurrences of node IDs in links, which should give us the degree.
egrep "\bL\b" | cut -f4 assembly_graph_with_scaffolds.gfa  | sort | uniq -c | sort | tail

# prints:
      
      5 3461076
      5 3461246
      5 3461286
      5 591392
      5 606626
      5 634272
      5 721568
      5 766622
      5 822586
      6 3455660

So we can see that there are a lot of links of segment 3455660. We can confirm this in the GFA:
grep 3455660 assembly_graph_with_scaffolds.gfa
# prints:
S       3455660 CTGATAAACGGCGGGAACAGAACCAACACTACGCGTTGCTCCGATCTCAACAACTTCTAAG   KC:i:2936
L       3455660 +       56698   -       55M
L       3455660 +       3455666 +       55M
L       3455660 +       3456290 -       55M
L       3455660 -       3455652 -       55M
L       3455660 -       3456398 -       55M
P       NODE_6_length_300278_cov_129.624456_1   3448698+,3456292-,3448662+,3448670+,3451896-,3460544-,3460388+,3434104-,3434124-,56698+,3455660-,3455652-,3460374+,3460806-,3460824-,3460806-,3460824-,3460806-,3460530-,3449022-,3454210+,3459744-,3459742-    *
P       NODE_8_length_282610_cov_136.723533_1   3460440-,3434818+,3276488+,3461076-,3461070-,3461062-,3461054-,3461046-,3461038-,3460426-,3461020+,3461022+,3461020+,3460516-,3455644+,3455652+,3455660+,3456290-,3448662+,3460546-,3275518-,3387950+,3456304+,3367628+,523826-,3418048-,3456388+,3460376-,3327686+,3455382-,3460586+,655086-,3275436+,3460670+,3278936+,3460218+,3460346+,3460498+,3459664+,3460496+,3459664+,3459200+,3459424-,3458356+,3455426-,3455420-     *
P       NODE_39_length_23557_cov_173.180283_1   3456398+,3455660+,3455666+,3418048-,3273044-,3274776+,3433484-,3433478-,3433470-,3433462-,3456812-,3460134+,3455826-,3319248+,3454700+,3421874+,3453038+,3459348-,3453046-,3234810-,3460596+,3460980+,3460654+  *

Extract sequences proximal to a node
Having extracted the high-degree node, we can use the output of grep above to find links to other sequences:
grep 3455660 assembly_graph_with_scaffolds.gfa | egrep "\bL"
# prints:
L   3455660 +   56698   -   55M
L   3455660 +   3455666 +   55M
L   3455660 +   3456290 -   55M
L   3455660 -   3455652 -   55M
L   3455660 -   3456398 -   55M

We can then go lookup S lines for those segments:
egrep "\bS\W56698" assembly_graph_with_scaffolds.gfa
# prints:
S   56698   GCATGTCGAGAAAACACCTTAGAAGTTGTTGAGATCGGAGCAACGCGTAGTGTTGGTTCTGTTCCCGCCGTT    KC:i:2067

Coverage filtering
Read count (RC tag) is an optional field for GFA, which does not exist in the linked file.
However, we do have k-mer count (KC), which you could use. Presumably high k-mer count would correlate well to coverage. I don't have a suggestion of a nifty tool that filters on this, I'd use a custom script.
Summary
Probably there exists some command line tool that does this stuff, but I don't know about it. It has to all exist inside assemblers, but it might be pretty hard to pull out.
If I had a little more time, I might write a simple tool that uses knowledge of the GFA structure to perform these kinds of data slicing operations, it seems pretty useful. It really shouldn't be hard to do this stuff!
A related approach, that I believe has a somewhat different focus, can be found here.

> **Comments:**
> **username_0**: Thanks Max! This is exactly the kind of insight I am looking for. I read the GFA specification but was honestly confused by it. This makes everything clear. Excellent answer  
> **username_0**: In my searching, I've also found [RGFATools](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5103826/), which provides a framework for parsing and editing GFA files. They include a method for computing the coverage of a segment given the length of a segment, the *k*-mer count of the segment, and the average *k*-mer count of all segments. I think I'll try my hand at implementing this for filtering  
> **username_1**: @acvill that is interesting- I knew about gfapy (https://academic.oup.com/bioinformatics/article/33/19/3094/3883940) though I have never used it. If you find a better way to do it with that rather than all the kludgey shell stuff, consider posting your own answer.  

---

### Answer by username_0 on 2022-02-12. Score: 3

**Edit**: I've written a short blog post describing my use of nodeSeqs.sh to investigate the functions of sequences at the nodes of short-read assembly graphs for a set of human oral microbiome samples.

<https://albertvill.com/posts/nodeseqs_oral/>

[TL;DR](https://www.merriam-webster.com/dictionary/TL%3BDR) -- transposases and recombinases are enriched, as expected.

---

Based on [Maximilian's excellent answer](https://bioinformatics.stackexchange.com/a/18528/3967), I have created a script to address my own question.

<https://github.com/acvill/nodeSeqs>

Here are the details on how it works:

> Link lines (beginning with `L`) are extracted from the GFA file, and segment occurrences are counted from both the `From` and `To` fields. To ensure that the final sequences capture the overlapping (linked) portions of segments, links are processed as segment-orientation pairs. Segment IDs with at least the specified minimum degree (`-d`) are written to `temp_XXXXX/segments.txt`, and the complementary segment lines (beginning with `S`) are pulled. *k*-mer counts (`KC` tag) are used to compute coverage for each segment using the approximation from [Gonnella & Kurtz 2016](https://dx.doi.org/10.7717%2Fpeerj.2681); rewritten:
>
> [![coverage estimation equation](https://i.sstatic.net/hjwHC.png)](https://i.sstatic.net/hjwHC.png)
>
> where *K**S* is the *k*-mer count for a segment, *L**S* is the length of a segment, and *k* is the *k*-mer size used in de Bruijn graph construction (`-k`). For (meta)SPAdes, graphs are built iteratively using an increasing *k*-mer size, so *k* is equal to the largest *k*-mer size used. Sequences are extracted corresponding to segments whose coverage is at least the specified minimum coverage (`-m`). Sequences longer than the specified context (`-c`) are truncated to the context length with respect to the link orientation. Including the `-t` flag will automatically remove the temporary folder when the program is finished, which holds intermediate files and `log.txt`.

The code is indeed "shell stuff", though hopefully it is not too "kludgey". ;-)

> **Comments:**
> **username_1**: +1 for putting in the k-mer coverage approximation, which I neglected. Not too kludgey, unlike my answer! I think that you could accept your own answer if you liked, it is a better and more useful approach than mine, and fully answers your question.  
> **username_0**: @MaximilianPress thanks for the feedback. Since your answer provided the framework to create my answer, I'm happy to keep the accept in place.  

---

