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.
Given three short sequences,
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
- 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 memdoes not output multiply-mapped in a predictable manner that fits this use case https://www.biostars.org/p/304614/
blatseems 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?