4
$\begingroup$

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?
Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ 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. $\endgroup$
    – user172818
    Sep 6, 2018 at 14:08
  • $\begingroup$ 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. $\endgroup$
    – conchoecia
    Sep 7, 2018 at 6:13
  • 1
    $\begingroup$ 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. $\endgroup$
    – user172818
    Sep 7, 2018 at 11:52
  • $\begingroup$ 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. $\endgroup$
    – conchoecia
    Sep 8, 2018 at 2:31
  • 1
    $\begingroup$ "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. $\endgroup$
    – user172818
    Sep 9, 2018 at 0:24

1 Answer 1

Reset to default
3
$\begingroup$

LAST 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 here), 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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – conchoecia
    Sep 10, 2018 at 3:12

Your Answer

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

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