Multipath alignments and vg mpmap

This wiki page describes how to produce multipath alignments using vg mpmap. For a general introduction to the concept of a multipath alignment, see Multipath alignments and the GAMP format.

Producing multipath alignments with vg mpmap

The vg toolkit has a subcommand implemented that can produce multipath alignments against a graph. It is best-optimized for short read sequencing data. However, it is functional for other sequencing types as well.

The command line interface for vg mpmap is mostly similar to vg map. Like vg map, it requires an XG index of the graph and a GCSA index for substring search (see Index Construction).

# -x for XG index, -g for GCSA index, -n for genomic DNA, -f for FASTQ reads
vg mpmap -x graph.xg -g graph.gcsa -n DNA -f reads.fq > mapped_reads.gamp

Paired end data

Paired end reads are also supported. Read pairs can be given either as two matched FASTQ files or as a FASTQ file with interleaved pairs.

# second -f for the matched FASTQ file containing the read pairs
vg mpmap -x graph.xg -g graph.gcsa -n DNA -f reads_1.fq -f reads_2.fq > mapped_reads.gamp

# -i indicates interleaved read pairs
vg mpmap -x graph.xg -g graph.gcsa -n DNA -f read_pairs.fq -i > mapped_reads.gamp

Producing better multipath alignment topologies

vg mpmap always tries to identify "neighboring" parts of the graph that it should align the same part of the read to. However, this is a challenging problem without any additional information. It can do a much better job with some global information about the structure of the graph. In particular, it wants to know the graph's "snarl decomposition", which is essentially the structure of "bubbles" in the graph (see our paper or preprint for a more rigorous explanation).

The snarl decomposition should first be computed using the vg snarls command:

vg snarls graph.vg > graph.snarls

The snarl decomposition can then be passed to vg mpmap with the -s flag:

# -s for the snarl decomposition
vg mpmap -x graph.xg -g graph.gcsa -s graph.snarls -f reads.fq > mapped_reads.gamp

Using vg mpmap as a single path aligner

Since there is a simple conversion from GAMP to GAM, it is possible to use vg mpmap as a sequence-to-path mapper as well. To do so, there is an option -F GAM to execute in "single path mode", which produces GAM output and implements some heuristics that are more appropriate for single path alignments.

# -F GAM to produce GAM instead of GAMP
vg mpmap -x graph.xg -g graph.gcsa -f reads.fq -F GAM > mapped_reads.gam

The -F argument can also project alignments onto the reference path with the SAM and BAM options.

Configuring vg mpmap for different types of data

Like all read mappers, vg mpmap uses some computational heuristics to speed up read mapping. Invariably, tuning these heuristics requires making some assumptions about the kind of data the mapper is likely to see. By default, vg mpmap expects to see short-read NGS data for genomic DNA. However, the heuristics can be easily tuned using three intuitive presets that can be selected using command line options:

• -n/--nt-type: RNA for transcriptomic data (default), DNA for genomic DNA.
• -l/--read-length: very-short, short (default), or long for reads that are approximately <50 bp, 50-500 bp, and >500 bp, respectively.
• -e/--error-rate: low (default) or high for reads that have a base-wise accuracy of approximately PHRED >20 or PHRED <20, respectively.

Additional considerations for using vg mpmap for to map RNA-seq reads to splice-variation graphs

vg mpmap is best optimized for mapping RNA sequencing data to a splice-variation graph, such as those produced by vg rna. For optimal performance, one additional index is needed (for efficiently measuring minimum distances, see our preprint or paper). The total set of indexes can be constructed using default best-practices with vg autoindex --workflow mpmap. Alternatively, the indexing pipeline can be built manually at the command line.

To build the distance index, first compute the snarl decomposition of the graph, and then use the vg index.

# compute snarls, note: use of -T differs from above
vg snarls -T graph.xg > graph.trivial.snarls
# compute the distance index
vg index -x graph.xg -s graph.trivial.snarls -j graph.dist

A somewhat annoying point (as of April 2020): the snarls computed for the distance index must include "trivial snarls" using the -T option, which consist of a single edge along a non-branching path in the graph. However, the snarls we supply to vg mpmap (for producing good multipath alignments, see above) ideally should not include the trivial snarls, as they slow down the mapping somewhat. We hope to make this process simpler in a future release.

Now that the distance index is computed, we supply it to vg mpmap using the -d flag. The distance index encodes all of the information contained in the snarl decomposition. If -d is supplied, -s is unnecessary.

vg mpmap -x graph.xg -g graph.gcsa -d graph.dist -n rna -f reads.fq > mapped_reads.gamp

As an aside, all of the modifications to the vg mpmap algorithm discussed above are also compatible with the distance index option.

Using vg mpmap to map reads to graphs with structural variants

It turns out that this problem is very similar to mapping RNA-seq reads. After all, an intron is a type of large-scale deletion (in the transcript). Accordingly, it is possible to effectively map reads to graphs of structural variants using the approach described in the previous section. The only difference is that -n DNA should be provided as an option.