Answers.Body: <p>Note also that the &quot;score&quot; can vary dramatically depending on whether -c is used or not. I made a short post about this here <a href="https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea" rel="nofollow noreferrer">https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea</a> it is also alluded to briefly in the PAF.md where score is a required field but &quot;If alignment is not available, column 10 and 11 are still required but may be highly inaccurate.&quot;</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-20T05:34:38.300Z
Answers.LastActivityDate: 2023-12-20T05:34:38.300Z
Answers.UserDisplayName: Colin D
Answers.Body: <p>The program's author pointed me to the <a href="https://github.com/lh3/minimap2/blob/master/FAQ.md" rel="noreferrer">Minimap2 FAQ</a>:</p>
  <blockquote>
  <p><strong>1. Alignment different with option <code>-a</code> or <code>-c?</code></strong></p>
  <p>Without <code>-a</code>, <code>-c</code> or <code>--cs</code>, minimap2 only finds <em>approximate</em>
  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.</p>
  </blockquote>
  <hr />
  <p>So, to specifically address the main questions asked here:</p>
  <ol>
  <li><strong><em>Why are there fewer and non-identical computed alignments when the CIGAR tags are added?</em></strong></li>
  </ol>
  <p><code>minimap2</code> alignment accuracy differs when running <code>minimap2 -k15 -w10</code> <em>vs.</em> when running <code>minimap2 -c -k15 -w10</code>. It is actually <em>not unexpected</em> that the alignment coordinates in each PAF file are only <em>similar</em>, not identical, and an increase or decrease in <em>the number of alignments</em> between the outputs is simply an artefact of this accuracy discrepancy.</p>
  <ol start="2">
  <li><strong><em>Which alignment (default vs. with CIGAR tags) is better or preferable?</em></strong></li>
  </ol>
  <p>Although both alignments report very similar sequence relationships, adding the <code>-c</code> flag for CIGAR tag generation is preferable to get more accurate base alignments.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-18T01:12:34.073Z
Answers.Text: Thanks for posting the answer and accepting it. Kudos to you @Gawain
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-15T01:15:09.790Z
Answers.LastActivityDate: 2023-12-15T01:15:09.790Z
Answers.UserDisplayName: Gawain
Body: <p>I am using <a href="https://github.com/lh3/minimap2" rel="nofollow noreferrer">Minimap2</a> (v2.26-r1175) in Linux to generate a sequence alignment between the <em>Streptomyces coelicolor</em> A3(2) chromosome (<code>ref.fa</code>) and the <em>Mycobacterium tuberculosis</em> chromosome (<code>query.fa</code>). My desired output is a <a href="https://github.com/lh3/miniasm/blob/master/PAF.md" rel="nofollow noreferrer">PAF (Pairwise mApping Format)</a> file.</p>
  <p>The general way to align reference and query sequences with Minimap2 is the following:</p>
  <pre><code>minimap2 ref.fa query.fa &gt; approx_mapping.paf
  # OR equivalently...
  minimap2 -k15 -w10 ref.fa query.fa &gt; approx_mapping.paf
  </code></pre>
  <p>You can get Minimap2 to generate custom CIGAR tags (<code>cg:</code>) in the PAF by adding the <code>-c</code> flag:</p>
  <pre><code>minimap2 -c ref.fa query.fa &gt; approx_mapping_cigar.paf
  # OR equivalently...
  minimap2 -c -k15 -w10 ref.fa query.fa &gt; approx_mapping_cigar.paf
  </code></pre>
  <hr />
  <p>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 <code>dim()</code>:</p>
  <pre class="lang-r prettyprint-override"><code>file1 &lt;- &quot;approx_mapping.paf&quot;
  paf_basic &lt;- read.table(file1, sep = &quot;\t&quot;, fill = TRUE,
                          col.names = paste0(&quot;V&quot;,seq_len(max(count.fields(file1, sep = &quot;\t&quot;)))))
  dim(paf_basic)
  [1] 205  18
  
  file2 &lt;- &quot;approx_mapping_cigar.paf&quot;
  paf_cigar &lt;- read.table(file2, sep = &quot;\t&quot;, fill = TRUE,
                          col.names = paste0(&quot;V&quot;,seq_len(max(count.fields(file2, sep = &quot;\t&quot;)))))
  dim(paf_cigar)
  [1] 200  24
  </code></pre>
  <p><em>Note the unqiue way we have to read the PAF files into R using <code>read.table()</code> to <a href="https://stackoverflow.com/a/18922914">ensure no columns are cut off and added to new rows</a>. This is because <code>read.table()</code> 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 <strong>with CIGAR tags</strong>, the number of columns can vary per line, which makes this necessary.</em></p>
  <p>With the <strong>default settings</strong>, I get <strong>205 observations</strong> and with the <strong>addition of the CIGAR flag</strong>, I get <strong>200 observations</strong>.</p>
  <p><strong>Why do I get fewer computed alignments when the CIGAR tags are added?</strong> I would expect the alignments between the reference and query to be identical between these runs.</p>
  <hr />
  <p>Comparing the PAF files (directly generated by <code>minimap2</code>) after sorting (<code>sort -k3,3 -k4,4</code>) using <code>vimdiff</code> reveals that the alignments computed between the default and CIGAR setting runs are largely <em>similar</em>, but not identical:</p>
  <p><a href="https://i.sstatic.net/bGQze.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/bGQze.png" alt="vimdiff_of_default_and_cigar_minimap2_output_files" /></a></p>
  <p><em>Note that the image above only shows the first 44 rows of both files.</em></p>
  <p>Due to the small differences seen in alignment coordinates between these two dataframes (see fields <code>$3</code> and <code>$4</code> in the <code>vimdiff</code> screenshot), it is <strong>NOT</strong> trivial to determine which five observations (rows) are new alignments computed when using the CIGAR flag <code>-c</code>.</p>
  <hr />
  <p>Plotting these dataframes confirms that the <code>paf_basic</code> and <code>paf_cigar</code> objects do, in fact, contain some different sequence alignments from each other:</p>
  <pre class="lang-r prettyprint-override"><code>library(&quot;ggplot2&quot;)
  ggplot(data = paf) +
    geom_segment(aes(x = V8, xend = V9, y = V3, yend = V4), size = 2) +
    labs(title = &quot;paf_name&quot;, x = &quot;Reference (Mb)&quot;, y = &quot;Query (Mb)&quot;) +
    coord_fixed(ratio = 1) +
    theme_bw()
  </code></pre>
  <p><a href="https://i.sstatic.net/gu3ON.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gu3ON.png" alt="paf_basic_dot_plot" /></a>
  <a href="https://i.sstatic.net/dfB3N.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/dfB3N.png" alt="paf_cigar_dot_plot" /></a></p>
  <p>So, it appears as though <code>minimap2 -c</code> produces a slightly different alignment than <code>minimap2</code>, which is surprising to me and begs the question: <strong>Which alignment is better or preferable?</strong></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-12T18:11:33.280Z
Comments.Text: 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.
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-13T05:45:26.447Z
Comments.Text: 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()`.
Comments.UserDisplayName: Gawain
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-13T18:51:43.517Z
Comments.Text: 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](https://github.com/lh3/minimap2/issues) on minimap2 github, Heng is pretty responsive. (If you do so, be sure to link to the issue here!)
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-14T04:04:12.467Z
Comments.Text: 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!
Comments.UserDisplayName: Gawain
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-15T03:57:25.260Z
Comments.Text: Note also that the scores are very different with and without -c. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea
Comments.UserDisplayName: Colin D
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-15T04:38:13.077Z
Comments.Text: Awesome analysis Colin! Interesting to see how much of an effect `-c` has on identity scores of different length alignments. Thanks for adding it to the conversation here.
Comments.UserDisplayName: Gawain
Tags: <phylogenetics><minimap2>
Title: Why are Minimap2 alignments different with CIGAR generation flag?
ThreadText: # Why are Minimap2 alignments different with CIGAR generation flag?
  
  (Asked by: username_0 on 2023-12-12)
  
  I am using [Minimap2](https://github.com/lh3/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)](https://github.com/lh3/miniasm/blob/master/PAF.md) 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()`:
  
  ```r
  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](https://stackoverflow.com/a/18922914). 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](https://i.sstatic.net/bGQze.png)](https://i.sstatic.net/bGQze.png)
  
  *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:
  
  ```r
  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](https://i.sstatic.net/gu3ON.png)](https://i.sstatic.net/gu3ON.png)
  [![paf_cigar_dot_plot](https://i.sstatic.net/dfB3N.png)](https://i.sstatic.net/dfB3N.png)
  
  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?**
  
  **Tags:** phylogenetics, minimap2
  
  > **Comments:**
  > **username_1**: 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.  
  > **username_0**: 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()`.  
  > **username_1**: 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](https://github.com/lh3/minimap2/issues) on minimap2 github, Heng is pretty responsive. (If you do so, be sure to link to the issue here!)  
  > **username_0**: Thanks for the idea username_1 -- 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!  
  > **username_2**: Note also that the scores are very different with and without -c. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea  
  > **username_0**: Awesome analysis username_2! Interesting to see how much of an effect `-c` has on identity scores of different length alignments. Thanks for adding it to the conversation here.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2023-12-15. Score: 6 | Accepted answer
  
  The program's author pointed me to the [Minimap2 FAQ](https://github.com/lh3/minimap2/blob/master/FAQ.md):
  
  > **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.
  
  2. ***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.
  
  > **Comments:**
  > **username_3**: Thanks for posting the answer and accepting it. Kudos to you @username_0  
  
  ---
  
  ### Answer by username_2 on 2023-12-20. Score: 0
  
  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."
  
  
  ---
  
  