We sequenced human mtDNA for multiple individuals. Our method involved making 2 partially overlapping PCR amplicons of mtDNA and then sequencing them on Illumina MiSeq (paired-end, ~350 average insert size).

I removed adapters, mapped the reads to hg38 reference using bwa and I am observing a strange coverage drop at position 3,234. It is present at this exact position in all samples (multiple separate batches) but in some samples the coverage drops from several hundred fold to 10-15x while in others it only drops to about half the coverage of the neighboring region:

This position is in the middle of one of the PCR amplicons and not near any of the primers used to make the amplicons.

I am trying to figure out what is happening.

1) Does the shape of the coverage track indicate anything? I am confused why the coverage drops so suddenly at position 3,234 (in a straight vertical line) but then gradually increases from almost 0 starting at position 3,235. How do I interpret that?

2) What should I try to figure this out (things I have already done are below)?. I'd like to solve this in silico but if I can't, will PCR amplify the region.

Here is what I have done so far:

1) checked primer specificity using BLAST - none of the primer sequences match anywhere besides the intended positions on human mtDNA

2) de novo assembled several samples - in all cases the region downstream of position 3,234 didn't assemble correctly (when I map reads back to the assembly, that region looks messy and is not uniformly covered; when I BLAST the part of the assembly downstream of 3,234 it aligns to various parts of mtDNA that are not supposed to be near position 3,234 but there is no clear pattern between different samples

3) looked at soft-clipped reads on both sides of position 3,234. Some soft-clipped parts of reads cannot be found anywhere at all, others map to various parts of mtDNA - again no clear pattern.

4) for reads that map near position 3,234, I looked where the other read in the pair maps - again, there is no clear pattern

5) looked at incorrectly oriented reads - there are some throughout the alignment and for the ones where one read is near position 3,234 again there is no clear pattern to where the second read in the pair is.

Thank you very much in advance!

First thing to do is to investigate the region in a genome browser. Here is how your region looks in UCSC. Your coordinate overlaps an annotated repeat: tRNA-Leu-TTA(m) seemingly a tandem repeat of variable length. Your read pairs likely do not span the real repeat region which is actually longer than the sequence depicted in the assembly, or if they do span it, the resulting read sequence is not identical to the reference. This is also the reason why you couldn't assemble this region.

So, when you look at the genome at this spot, for which the actual assembly gives:

>hg38_rmsk_tRNA-Leu-TTA(m) range=chrM:3230-3307 5'pad=0 3'pad=0 strand=+ repeatMasking=lower
gttaagatggcagagcccggtaatcgcataaaacttaaaactttacagtc
agaggttcaattcctcttcttaacaaca

it could be in fact 10 or 100 times that sequence, and with variability between individuals as well. No read pair could span the beginning of the repeat and still align to the reference.

If you need the exact sequence length at this site, you could do nanopore sequencing. I estimate that the actual repeat rate is not much more than 10x - 20x, limited by the fragment size of the specific polymerase used during PCR (pcr_length = 1400bp+78*x if there is only this repeat). You could maybe also measure your real fragment length using a southern blot.

You've got some kind of indel going on there. That's what causes that coverage pattern.

de novo assembled several samples - in all cases the region downstream of position 3,234 didn't assemble correctly

Not 'correctly'? Or not what you expect?