Extract sequence context of high-degree nodes in assembly graphs

Question

I often use metaSPAdes 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, 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 ...

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.

Answer 1

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.

Answer 2

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 -- transposases and recombinases are enriched, as expected.

---

Based on Maximilian's excellent answer, 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; rewritten:

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". ;-)