Answers.Body: <p><a href="http://last.cbrc.jp/" rel="nofollow noreferrer">LAST</a> has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see <a href="https://twitter.com/gringene_bio/status/1037561956066975744" rel="nofollow noreferrer">here</a>), as well as abundant base calling differences present in single-molecule sequencing.</p>
  
  <p>Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.</p>
  
  <p>I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly <em>is</em> possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-10T03:12:35.613Z
Answers.Text: Thanks for your comments, @gringer. The limitations in the problem are the limitations. The best I can do in this case is minimize false positive hits. I'll see what I can do with LAST!
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-09T11:56:58.633Z
Answers.LastActivityDate: 2018-09-09T11:56:58.633Z
Answers.UserDisplayName: gringer
Body: <h3>The dilemma:</h3>
  <p>I have around one million short sequences (21 bp to several 100s of basepairs) for which I need to identify all occurrences of in 20-30x coverage noisy long reads (both pacbio and ONT). All of the short sequences and long reads are derived from the same individual, and there are enough short sequences given the organism's genome size that each long read should have at least one or multiple short sequences.</p>
  <h3>An example:</h3>
  <p>Given three short sequences, <code>AA</code> &amp; <code>BBBB</code> &amp; <code>CCCCCC</code> , how to identify all occurrences of these sequences in long, uncorrected reads.</p>
  <pre><code>Sequences to identify in long reads:
    - AA
    - BBBB
    - CCCCCC
  
  Long Reads:
    - 1: ------AA-----------------------CCCCCC------
    - 2: ------AA--------BBBB-----------------------
    - 3: ---BBBB-----------CCCC--
    - 4: ------------------AA-----
  </code></pre>
  <p>In this example, the output sam file will contain high MAPQ hits for short sequence <code>AA</code> to reads (1, 2, 3), for short sequence <code>BBBB</code> to reads (2, 3). For short sequence <code>CCCCCC</code> there will be a high MAPQ hit to read 1 and a lower MAPQ hit to read 3</p>
  <h3>Limitations:</h3>
  <ul>
  <li>Due to the requirements of this project, I cannot perform a de novo assembly then map the short sequences to the assembly output. It is necessary to directly map the short sequences to the long reads.</li>
  <li>I cannot correct the long reads before I map/align the short sequences to them.</li>
  <li>The alignment/mapping should sensitively detect the short sequences in the long reads.</li>
  <li>The short sequences may be anywhere from 1Mbp to 30Mbp in total bases.</li>
  <li>The long reads may be anywhere from 3Gbp to 90Gbp in total bases. So the mapper/alignment technique is hopefully fast. Sensitivity is more important though.</li>
  </ul>
  <h3>Research so far:</h3>
  <ul>
  <li><code>bwa mem</code> does not output multiply-mapped in a predictable manner that fits this use case <a href="https://www.biostars.org/p/304614/" rel="nofollow noreferrer">https://www.biostars.org/p/304614/</a></li>
  <li><code>blat</code> seems to output many hits for a single short sequence, but there are no out-of-the-box parameters for mapping short reads to long reads.</li>
  <li>Maybe multiple sequence alignment is the most sensitive, but wouldn't that entail running (no. of short seqs) * (no. of long reads) alignments?</li>
  </ul>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-06T14:08:13.000Z
Comments.Text: Don't map short reads to noisy long reads. This won't work. As you already have decent coverage, assemble first (or correct long reads by themselves with canu) and then map. It will be faster and give you *much* better result.
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-07T06:13:21.200Z
Comments.Text: Hi user17818 - unfortunately that isn't an option in this scenario based on the requirements for this project. I'll update my question to better clarify that de novo assembly, then mapping the short sequences is not an option.
Comments.UserDisplayName: conchoecia
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-07T11:52:20.750Z
Comments.Text: This is probably an intermediate problem of a larger project. If so, you need to change the strategy. Mapping a 21bp oligo to 90Gb noisy reads will lead to many false positives and negatives. This is a theoretical limitation that no direct mapping algorithms can get around.
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-08T02:31:19.080Z
Comments.Text: You're right - there will be false negatives and positives. My goal is to minimize them by finding the best mapping strategy to get multiple mappings of each query. The size range of the illumina-derived "short" query sequences vary from 21 to hundreds/thousands of bp. The chances of correctly mapping the short queries to a long read will increase as the lengths of the illumina-derived queries increase. The kmer space of even a 31-mer is 4.61x10^18, or around 9 orders of magnitude larger than the kmer space of a human genome. A few indels in the long reads should still be uniquely mappable.
Comments.UserDisplayName: conchoecia
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-09T00:24:36.613Z
Comments.Text: "Uniquely mapped" doesn't mean the mapping is correct. This is particularly true when you look for all partial mappings. For >100bp sequences, you may use minimap2, taking Illumina sequences as the reference. It won't work well with shorter sequences, though.
Comments.UserDisplayName: user172818
Tags: <phylogenetics><software-recommendation><long-reads><homework>
Title: How to map short sequences to long reads, recovering all multiply-mapped high-quality matches
ThreadText: # How to map short sequences to long reads, recovering all multiply-mapped high-quality matches
  
  (Asked by: username_0 on 2018-09-06)
  
  ### The dilemma:
  
  I have around one million short sequences (21 bp to several 100s of basepairs) for which I need to identify all occurrences of in 20-30x coverage noisy long reads (both pacbio and ONT). All of the short sequences and long reads are derived from the same individual, and there are enough short sequences given the organism's genome size that each long read should have at least one or multiple short sequences.
  
  ### An example:
  
  Given three short sequences, `AA` & `BBBB` & `CCCCCC` , how to identify all occurrences of these sequences in long, uncorrected reads.
  
  ```
  Sequences to identify in long reads:
    - AA
    - BBBB
    - CCCCCC
  
  Long Reads:
    - 1: ------AA-----------------------CCCCCC------
    - 2: ------AA--------BBBB-----------------------
    - 3: ---BBBB-----------CCCC--
    - 4: ------------------AA-----
  
  ```
  
  In this example, the output sam file will contain high MAPQ hits for short sequence `AA` to reads (1, 2, 3), for short sequence `BBBB` to reads (2, 3). For short sequence `CCCCCC` there will be a high MAPQ hit to read 1 and a lower MAPQ hit to read 3
  
  ### Limitations:
  
  * Due to the requirements of this project, I cannot perform a de novo assembly then map the short sequences to the assembly output. It is necessary to directly map the short sequences to the long reads.
  * I cannot correct the long reads before I map/align the short sequences to them.
  * The alignment/mapping should sensitively detect the short sequences in the long reads.
  * The short sequences may be anywhere from 1Mbp to 30Mbp in total bases.
  * The long reads may be anywhere from 3Gbp to 90Gbp in total bases. So the mapper/alignment technique is hopefully fast. Sensitivity is more important though.
  
  ### Research so far:
  
  * `bwa mem` does not output multiply-mapped in a predictable manner that fits this use case <https://www.biostars.org/p/304614/>
  * `blat` seems to output many hits for a single short sequence, but there are no out-of-the-box parameters for mapping short reads to long reads.
  * Maybe multiple sequence alignment is the most sensitive, but wouldn't that entail running (no. of short seqs) \* (no. of long reads) alignments?
  
  **Tags:** phylogenetics, software-recommendation, long-reads, homework
  
  > **Comments:**
  > **username_1**: Don't map short reads to noisy long reads. This won't work. As you already have decent coverage, assemble first (or correct long reads by themselves with canu) and then map. It will be faster and give you *much* better result.  
  > **username_0**: Hi user17818 - unfortunately that isn't an option in this scenario based on the requirements for this project. I'll update my question to better clarify that de novo assembly, then mapping the short sequences is not an option.  
  > **username_1**: This is probably an intermediate problem of a larger project. If so, you need to change the strategy. Mapping a 21bp oligo to 90Gb noisy reads will lead to many false positives and negatives. This is a theoretical limitation that no direct mapping algorithms can get around.  
  > **username_0**: You're right - there will be false negatives and positives. My goal is to minimize them by finding the best mapping strategy to get multiple mappings of each query. The size range of the illumina-derived "short" query sequences vary from 21 to hundreds/thousands of bp. The chances of correctly mapping the short queries to a long read will increase as the lengths of the illumina-derived queries increase. The kmer space of even a 31-mer is 4.61x10^18, or around 9 orders of magnitude larger than the kmer space of a human genome. A few indels in the long reads should still be uniquely mappable.  
  > **username_1**: "Uniquely mapped" doesn't mean the mapping is correct. This is particularly true when you look for all partial mappings. For >100bp sequences, you may use minimap2, taking Illumina sequences as the reference. It won't work well with shorter sequences, though.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2018-09-09. Score: 3 | Accepted answer
  
  [LAST](http://last.cbrc.jp/) has given the best results for me when I've tried to do this, although I agree with @username_1 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see [here](https://twitter.com/gringene_bio/status/1037561956066975744)), as well as abundant base calling differences present in single-molecule sequencing.
  
  Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.
  
  I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly *is* possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.
  
  > **Comments:**
  > **username_0**: Thanks for your comments, @username_2. The limitations in the problem are the limitations. The best I can do in this case is minimize false positive hits. I'll see what I can do with LAST!  
  
  ---
  
  