Understanding PAF format

Question

Does anyone have any experience with PAF files (minimap2 output)? I don’t fully understand what my output means.

I’ve looked at the minimap2 explanation (https://github.com/lh3/miniasm/blob/master/PAF.md) but I want to use the data here to go back and extract regions in the genome that the gene target mapped to. I just don't know what specific columns to look at to accomplish this.

Here is a little snippet of a PAF file I have:

contig_25   522460  514241  515712  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_26   413461  403718  406719  +   GENE_TARGET 4988    17      3649    1762  3633  60
contig_44   24218   16005   17476   +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_52   244411  236199  237670  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_57   1043583 1025163 1026895 +   GENE_TARGET 4988    3046    4986    1376  2059  60

Is this saying that the gene target mapped to coordinates 514241 through 515712 on contig_25 but only coordinates 1962 through 3649 of the gene actually mapped to that contig? I guess I’m mostly confused on the descriptions for the coordinates and what they actually mean.

I think for my use case I have to use the gene target as the reference sequence as I am trying to extract homologous regions to sequences of interest. I.e., my project involves 90+ assemblies and I would like to extract the regions in the contigs where they map to the gene target. I used minimap2 to map the assemblies to the gene target, used sort and awk to get the best/high quality matches, converted PAF to BED format, and then used seqtk subseq to extract the sequences from the respective assemblies.

This is what I get when I use the gene target as a reference sequence:

GENE    4988    17      3649    +   contig_26   413461  403718  406719  1762    3633    60
...
GENE    4988    3046    4986    +   contig_57   1043583 1025163 1026895 1376    2059    60

Would I use columns 8 and 9 as the region range to go back and extract the sequences from?

Answer 1

As @user172818 mentioned in a comment to your question, it looks like you might have used the query as the reference; they should be switched round.

A minimap2 mapping command is typically represented as follows:

minimap2 reference.fa query.fa > results.paf

This follows the "indexed reference first" command line order used in the overwhelming majority of mapping tools. In the minimap2 PAF output (which you have linked), the first column represents the query sequence name, and the 6th column represents the target sequence name (i.e. the opposite order from the command line order, also following with established convention).

Based on the desired mapping operation that you have described (i.e. "I want to use the data here to go back and extract regions in the genome that the gene target mapped to"), your minimap2 command should probably look like this:

minimap2 genome_contigs.fa GENE_TARGET.fa > results.paf

This will appropriately map the target region to the regions in the genome, and then columns 8 and 9 will represent "regions in the genome that the gene target mapped to."

Minimap2 (and other mappings programs) match query sequences to a reference, with the usual expectation being that the output will contain a single primary mapping per input query sequence (sometimes with secondary mappings to other locations). In other words, they try to find one match to any reference sequence for each query. Changing the order of mapping to target-as-reference means that the program will attempt to find one match to the target sequence(s) for each input reference sequence, which is usually not what is desired.

In any case, the columns in the PAF file correspond to the closest sequence ID to the left of them (which is, confusingly, ordered the opposite way round to the command line order):

$ minimap2 [genome_contigs.fa] {GENE_TARGET.fa}
{GENE 4988 17   3649} + [contig_26 413461  403718  406719 ] 1762 3633 60
...
{GENE 4988 3046 4986} + [contig_57 1043583 1025163 1026895] 1376 2059 60

$ minimap2 [GENE_TARGET.fa] {genome_contigs.fa}
{contig_26 413461  403718  406719 } + [GENE 4988 17   3649] 1762 3633 60
...
{contig_57 1043583 1025163 1026895} + [GENE 4988 3046 4986] 1376 2059 60