Get number of reads with a single, (almost) exact match to the full length of a reference sequence

Question

I can't find an answer to this in previous questions, so hoping someone can help me now.

We have nanopore sequenced a PCR product, and I've filtered our reads to +/- 10 bp of the expected product length (~4,400bp), leaving ~55,000 reads.

We now want to work out how many of the reads correspond to each of 6 different potential arrangements. The section we're interested in is in the middle of the 4,400 bp PCR product, and each different potential arrangement is ~1,500 bp. We have reference sequences for each of the 1,500 bp arrangements. Basically, we want to know what portion of our starting cells had each of the different arrangements of genes in this operon (e.g arrangement 1=30%, arrangement 2=15%, etc.).

This seems like it should be a really easy thing to do using read mapping to count alignments for each reference separately, but I just can't get my head around it. No matter what options I've tried with minimap2 and samtools, every read maps every time on every different arrangement. I think the problem is that the genes in our reference sequences are all present in every read, just in different arrangements. So we need to only count alignments which align to the whole 1,500 bp reference sequence, not partial hits or hits with huge gaps in them.

I have tried the minimap2 options -r, -G and --secondary=no to restrict the output to single alignments with no large gaps, and combined that with samtools view -m 1000 to output only alignments longer than 1,000 bp (this is far stricter than it needs to be to only output alignments to any of the different potential arrangements). But still, every time, every read is mapping to every reference, usually multiple times.

What am I doing wrong?! Am I totally over-complicating this, or doing something obviously stupid? Almost feel like using blastn might be better at this point...

TL;DR: how do I use minimap2 and/or samtools (or any other tool better suited to this) to count how many of our ~4,400 bp reads contain only full-length, almost exact matches to our 1,500 bp reference sequence, excluding any partial hits?

My latest attempt (broken into individual steps, though I'd usually pipe the samtools steps together, also not bothered with samtools sort at this point):

minimap2 -ax map-ont -r10 -N5 --secondary=no --no-long-join reference.fasta reads.fastq > reads.minimap.sam

samtools view -m 1000 -b -o reads.minimap.bam reads.minimap.sam

samtools -c reads.minimap.bam

Answer 1

Because minimap2 is a minimiser-based mapper, it's less useful for determining mapping accuracy to single-base precision.

Where precision matters, you can use LAST, bearing in mind that it calculates local alignments, and can report multiple hits per query/reference pair. The usual approach I use for LAST mapping is as follows:

1. lastdb -uNEAR -R01 reference.fa reference.fa
2. last-train -Q 1 reference.fa <(cat reads.fastq | head -n 40000) > train.mat
3. lastal -p train.mat reference.fa reads.fastq > reads_vs_ref.maf

That final MAF file can be inspected for mapping. The default output produced from LAST is not very human-readable, so I've written a perl script that converts MAF results into a CSV format with reference-relative statistics that I have found useful:

maf2csv.pl reads_vs_ref.maf > reads_vs_ref.csv

This can be inspected in a spreadsheet viewer (e.g. LibreOffice Calc, or Excel) to look at the mapped locations, mapped fraction, identities, etc..

Alternatively, the MAF results can be converted to SAM/BAM format using maf-convert (bearing in mind the previously-mentioned caveats about multiple hits):

maf-convert sam reads_vs_ref.maf | samtools view -b -T reference.fa | samtools sort > reads_vs_ref.bam