4
$\begingroup$

I am using Minimap2 (v2.26-r1175) in Linux to generate a sequence alignment between the Streptomyces coelicolor A3(2) chromosome (ref.fa) and the Mycobacterium tuberculosis chromosome (query.fa). My desired output is a PAF (Pairwise mApping Format) file.

The general way to align reference and query sequences with Minimap2 is the following:

minimap2 ref.fa query.fa > approx_mapping.paf
# OR equivalently...
minimap2 -k15 -w10 ref.fa query.fa > approx_mapping.paf

You can get Minimap2 to generate custom CIGAR tags (cg:) in the PAF by adding the -c flag:

minimap2 -c ref.fa query.fa > approx_mapping_cigar.paf
# OR equivalently...
minimap2 -c -k15 -w10 ref.fa query.fa > approx_mapping_cigar.paf

Moving to an R environment [with the intention of downstream data visualization], you can quickly assess the number of observations in each object (a.k.a. rows in each dataframe) with dim():

file1 <- "approx_mapping.paf"
paf_basic <- read.table(file1, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file1, sep = "\t")))))
dim(paf_basic)
[1] 205  18

file2 <- "approx_mapping_cigar.paf"
paf_cigar <- read.table(file2, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file2, sep = "\t")))))
dim(paf_cigar)
[1] 200  24

Note the unqiue way we have to read the PAF files into R using read.table() to ensure no columns are cut off and added to new rows. This is because read.table() is designed to examine only the first five lines of an input file to determine the number of columns in the dataframe. Any columns on lines beyond the fifth line that exceed the initially determined count are appended to a new row. In the case of PAF files generated by Minimap2 with CIGAR tags, the number of columns can vary per line, which makes this necessary.

With the default settings, I get 205 observations and with the addition of the CIGAR flag, I get 200 observations.

Why do I get fewer computed alignments when the CIGAR tags are added? I would expect the alignments between the reference and query to be identical between these runs.


Comparing the PAF files (directly generated by minimap2) after sorting (sort -k3,3 -k4,4) using vimdiff reveals that the alignments computed between the default and CIGAR setting runs are largely similar, but not identical:

vimdiff_of_default_and_cigar_minimap2_output_files

Note that the image above only shows the first 44 rows of both files.

Due to the small differences seen in alignment coordinates between these two dataframes (see fields $3 and $4 in the vimdiff screenshot), it is NOT trivial to determine which five observations (rows) are new alignments computed when using the CIGAR flag -c.


Plotting these dataframes confirms that the paf_basic and paf_cigar objects do, in fact, contain some different sequence alignments from each other:

library("ggplot2")
ggplot(data = paf) +
  geom_segment(aes(x = V8, xend = V9, y = V3, yend = V4), size = 2) +
  labs(title = "paf_name", x = "Reference (Mb)", y = "Query (Mb)") +
  coord_fixed(ratio = 1) +
  theme_bw()

paf_basic_dot_plot paf_cigar_dot_plot

So, it appears as though minimap2 -c produces a slightly different alignment than minimap2, which is surprising to me and begs the question: Which alignment is better or preferable?

$\endgroup$
6
  • 1
    $\begingroup$ You've shown that you have 5 extra rows when you read it into R. However, you haven't actually shown what the extra rows are. If we could see them, we might be able to guess what's going on. I can't reproduce the issue on my own. Can you rule out the trivial explanation that R's read.table is parsing the PAF incorrectly? That seems most likely given your observation about tag-only rows. $\endgroup$ Dec 12, 2023 at 18:11
  • $\begingroup$ Great point! Unfortunately, there is no straightforward way to determine what the 5 extra rows are because the two data frames have no matching observations between them (checked by dplyr::semi_join(paf_basic, paf_cigar, by=c("V3","V4","V8","V9"))). Based on the graphs included in the question (a recent edit), it seems that many of the alignment coordinates are close, but not identical. Unfortunately, the dput() outputs for these data frames are too large for SO question character limits. However, I have verified that the *.paf files match the paf R objects produced by read.table(). $\endgroup$
    – Gawain
    Dec 13, 2023 at 5:45
  • 1
    $\begingroup$ have you looked at the two PAF files side by side in e.g. vimdiff? It should not be too hard to find the differences leading to 2 extra rows. But IDK exactly what is going on. If you've demonstrated that adding the -c option changes the alignments themselves, that certainly sounds wrong. I'd suggest posting an issue on minimap2 github, Heng is pretty responsive. (If you do so, be sure to link to the issue here!) $\endgroup$ Dec 13, 2023 at 18:51
  • 1
    $\begingroup$ Thanks for the idea Maximilian -- I hadn't heard of vimdiff before. I gave it a go and found it immensely useful. I've updated the question to include a snapshot of the two files sorted and side-by-side in the program which reveals just how similar (but not identical) the data frames are. I'll look into posting the issue on the minimap2 github. Thanks for taking time out of your day to ponder this issue and help troubleshoot with me. Cheers! $\endgroup$
    – Gawain
    Dec 14, 2023 at 4:04
  • 1
    $\begingroup$ Note also that the scores are very different with and without -c. I made a short post about this here gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea $\endgroup$
    – Colin D
    Dec 15, 2023 at 3:57

2 Answers 2

5
$\begingroup$

The program's author pointed me to the Minimap2 FAQ:

1. Alignment different with option -a or -c?

Without -a, -c or --cs, minimap2 only finds approximate mapping locations without detailed base alignment. In particular, the start and end positions of the alignment are impricise. With one of those options, minimap2 will perform base alignment, which is generally more accurate but is much slower.


So, to specifically address the main questions asked here:

  1. Why are there fewer and non-identical computed alignments when the CIGAR tags are added?

minimap2 alignment accuracy differs when running minimap2 -k15 -w10 vs. when running minimap2 -c -k15 -w10. It is actually not unexpected that the alignment coordinates in each PAF file are only similar, not identical, and an increase or decrease in the number of alignments between the outputs is simply an artefact of this accuracy discrepancy.

  1. Which alignment (default vs. with CIGAR tags) is better or preferable?

Although both alignments report very similar sequence relationships, adding the -c flag for CIGAR tag generation is preferable to get more accurate base alignments.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks for posting the answer and accepting it. Kudos to you @Gawain $\endgroup$
    – M__
    Dec 18, 2023 at 1:12
0
$\begingroup$

Note also that the "score" can vary dramatically depending on whether -c is used or not. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea it is also alluded to briefly in the PAF.md where score is a required field but "If alignment is not available, column 10 and 11 are still required but may be highly inaccurate."

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.