Conda

gapmm2: gapped alignment using minimap2

This tool is a wrapper for minimap2 to run spliced/gapped alignment, ie aligning transcripts to a genome. You are probably saying, yes minimap2 runs this with -x splice --cs option (you are correct). However, there are instances where the terminal exons from stock minimap2 alignments are missing. This tool detects those alignments that have unaligned terminal eons and uses edlib to find the terminal exon positions. The tool then updates the PAF output file with the updated information.

Rationale

We can pull out a gene model in GFF3 format that has a short 5’ terminal exon:

scaffold_9    funannotate    gene    408904    409621    .    -    .    ID=OPO1_006919;
scaffold_9    funannotate    mRNA    408904    409621    .    -    .    ID=OPO1_006919-T1;Parent=OPO1_006919;product=hypothetical protein;
scaffold_9    funannotate    exon    409609    409621    .    -    .    ID=OPO1_006919-T1.exon1;Parent=OPO1_006919-T1;
scaffold_9    funannotate    exon    409320    409554    .    -    .    ID=OPO1_006919-T1.exon2;Parent=OPO1_006919-T1;
scaffold_9    funannotate    exon    409090    409255    .    -    .    ID=OPO1_006919-T1.exon3;Parent=OPO1_006919-T1;
scaffold_9    funannotate    exon    408904    409032    .    -    .    ID=OPO1_006919-T1.exon4;Parent=OPO1_006919-T1;
scaffold_9    funannotate    CDS    409609    409621    .    -    0    ID=OPO1_006919-T1.cds;Parent=OPO1_006919-T1;
scaffold_9    funannotate    CDS    409320    409554    .    -    2    ID=OPO1_006919-T1.cds;Parent=OPO1_006919-T1;
scaffold_9    funannotate    CDS    409090    409255    .    -    1    ID=OPO1_006919-T1.cds;Parent=OPO1_006919-T1;
scaffold_9    funannotate    CDS    408904    409032    .    -    0    ID=OPO1_006919-T1.cds;Parent=OPO1_006919-T1;

If we then map this transcript against the genome, we get the following PAF alignment.

$ minimap2 -x splice --cs genome.fasta cds-transcripts.fa | grep 'OPO1_006919'
OPO1_006919-T1    543    13    543    -    scaffold_9    658044    408903    409554    530    530    60    NM:i:0    ms:i:530    AS:i:466    nn:i:0    ts:A:+    tp:A:P    cm:i:167    s1:i:510    s2:i:0    de:f:0    rl:i:0    cs:Z::129~ct57ac:166~ct64ac:235

The --cs flag in minimap2 can be used to parse the coordinates (below) and you can see we are missing the 5’ exon.

>>> cs2coords(408903, 13, 543, '-', ':129~ct57ac:166~ct64ac:235')
([(409320, 409554), (409090, 409255), (408904, 409032)],

So if we run this same alignment with gapmm2 we are able to properly align the 5’ terminal exon.

$ gapmm2 genome.fa cds-transcripts.fa | grep 'OPO1_006919'
OPO1_006919-T1    543    0    543    -    scaffold_9    658044    408903    409621    543    543    60    tp:A:P    ts:A:+    NM:i:0    cs:Z::129~ct57ac:166~ct64ac:235~ct54ac:13

>>> cs2coords(408903, 0, 543, '-', ':129~ct57ac:166~ct64ac:235~ct54ac:13')
([(409609, 409621), (409320, 409554), (409090, 409255), (408904, 409032)]

Usage:

gapmm2 can be run as a command line script:

$ gapmm2
usage: gapmm2 [-o] [-f] [-t] [-m] [-i] [-d] [-h] [--version] reference query

gapmm2: gapped alignment with minimap2. Performs minimap2/mappy alignment with splice options and refines terminal alignments with edlib.

Positional arguments:
  reference           reference genome (FASTA)
  query               transcipts in FASTA or FASTQ

Optional arguments:
  -o , --out          output in PAF format (default: stdout)
  -f , --out-format   output format [paf,gff3] (default: paf)
  -t , --threads      number of threads to use with minimap2 (default: 3)
  -m , --min-mapq     minimum map quality value (default: 1)
  -i , --max-intron   max intron length, controls terminal search space (default: 500)
  -d, --debug         write some debug info to stderr (default: False)

Help:
  -h, --help          Show this help message and exit
  --version           Show program's version number and exit

Python API

It can also be run as a python module. The module provides several functions for working with spliced alignments:

`aligner` function

The main function for aligning transcripts to a genome. It can write an output file in either PAF or GFF3. It returns a dictionary with alignment statistics.

>>> from gapmm2.align import aligner
>>> stats = aligner('genome.fa', 'transcripts.fa', out_fmt="gff3", output="output.gff3")
>>> stats
{'n': 6926, 'low-mapq': 0, 'refine-left': 409, 'refine-right': 63}

`cs2coords` function

This function parses the CIGAR string (cs) from minimap2 and converts it to genomic coordinates, identifying exons, introns, and other alignment features.

>>> from gapmm2.align import cs2coords
>>> cs2coords(408903, 0, 543, '-', ':129~ct57ac:166~ct64ac:235~ct54ac:13')
([(409609, 409621), (409320, 409554), (409090, 409255), (408904, 409032)], [(0, 13), (13, 248), (248, 414), (414, 543)], 0, 0, True)

Installation

You can install gapmm2 using pip:

pip install gapmm2

Or you can install the latest development version directly from GitHub:

pip install git+https://github.com/nextgenusfs/gapmm2.git

You can also install from conda:

conda install -c bioconda gapmm2

Dependencies

Gapmm2 requires the following Python packages:

mappy (Python bindings for minimap2)
edlib (for sequence alignment)
natsort (for natural sorting)

These dependencies will be automatically installed when you install gapmm2 using pip or conda. Note that I’ve recently seen some seqmentation faults from mappy, so as of v25.4.13 it will run minimap2 directly instead of mappy if minimap2 is installed.

Development

Testing

Gapmm2 includes a test suite that can be run using pytest. To run the tests, first install pytest:

pip install pytest pytest-cov

Then run the tests from the root directory of the repository:

python -m pytest tests/ --cov=gapmm2

Code Formatting

This project uses pre-commit to ensure code quality and consistency. The pre-commit hooks run Black (code formatter), isort (import sorter), and flake8 (linter).

To set up pre-commit:

Install pre-commit:

pip install pre-commit

Install the git hooks:

pre-commit install

(Optional) Run against all files:

pre-commit run --all-files

After installation, the pre-commit hooks will run automatically on each commit to ensure your code follows the project’s style guidelines.