Long-Read Mapping, Depth Analysis, and Read-Spanning Validation

Welcome. This lesson sits within the genome-reconstruction module: after generating a candidate hybrid assembly, you need to ask whether the long-read evidence actually supports it. For CPE genomics, this matters especially at plasmid junctions, repeat-rich resistance regions, and circular chromosome or plasmid boundaries, where an apparently complete sequence can still be incorrectly connected.

By the end of this lesson, you will be able to create a coordinate-sorted, indexed long-read alignment with minimap2 and samtools; calculate and interpret whole-replicon and local depth; and extract named long reads that genuinely span an assembly join or a region such as a carbapenemase cassette.

---

Mapping as an assembly-validation test

A long-read assembly is a hypothesis about the order and adjacency of bases. Mapping reads back to that assembly tests whether the original molecules are consistent with the hypothesis.

A read that aligns continuously across a proposed join provides evidence that the two flanking sequences occur on the same physical DNA molecule. This is particularly informative when evaluating:

• a join through a repeat or mobile element;
• the connection between a carbapenemase gene and a plasmid backbone;
• a suspected plasmid cointegrate;
• the junction between the last and first bases of a contig reported as circular;
• a chromosomal integration boundary.

It is important, however, to distinguish read-to-assembly consistency from fully independent validation. If the same reads contributed to the assembly, remapping them is not independent evidence that the sequence is biologically correct. It is still highly valuable for finding unsupported joins, coverage anomalies, competing placements, and possible collapsed repeats. Where possible, a held-out set of long reads is stronger validation evidence, although withholding reads reduces assembly coverage.

Before working with commands, spend a few minutes consolidating the distinction between depth and breadth.

For a reference sequence of length $L$, where $d_i$ is the number of qualifying aligned reads at position $i$, mean observed depth is:

$$\bar{D}=\frac{{\sum }_{i=1}^L{d}_i}{L}$$

Breadth is the fraction of positions with at least one qualifying alignment:

$$B_{>0}=\frac{\#\{i:d_i>0\}}{L}$$

These differ from the nominal long-read coverage calculated during read QC:

$$C_{\text{nominal}}=\frac{\text{total retained long read bases}}{\text{estimated genome size}}$$

Nominal coverage is an input-data estimate. Observed alignment depth tells you what portion of those data supports this particular assembly, after mapping and filtering.

For complete bacterial assemblies, assess depth per replicon, not just as one genome-wide number. A small multicopy plasmid can have substantially greater depth than the chromosome. Conversely, unusually high depth in a chromosomal repeat may signal repeat collapse or ambiguous mapping rather than a genuine increase in copy number.

---

Build an indexed long-read alignment

The essential inputs are:

1. A final or candidate assembly FASTA, ideally your Autocycler consensus after Polypolish.
2. The quality-filtered long-read FASTQ used for assembly, or preferably an explicitly held-out validation FASTQ if one is available.
3. A clear record of which read set was mapped and which assembly version was tested.

The mapping output should be a coordinate-sorted BAM plus index. BAM is a compressed binary alignment format; the index lets samtools and IGV retrieve reads from a narrow genomic interval without scanning the entire file.

Choose an alignment preset deliberately

minimap2 presets bundle alignment parameters for particular read technologies and expected error profiles. For many ONT read sets, -x map-ont remains an appropriate general-purpose choice. For highly accurate recent ONT reads, -x lr:hq may be preferable. PacBio HiFi reads require -x map-hifi.

For a typical ONT mapping-validation run, use a streamed command: minimap2 writes SAM to standard output, and samtools immediately sorts and compresses it to BAM. This avoids retaining a very large intermediate SAM file.

sample="ISO001"
asm="assemblies/${sample}.autocycler.polypolish.fasta"
reads="reads/${sample}.ont.filtered.fastq.gz"

mkdir -p validation logs tmp

samtools faidx "$asm"

set -o pipefail

minimap2 \
  -t 16 \
  -a \
  -x map-ont \
  --MD \
  "$asm" "$reads" \
  2> "logs/${sample}.minimap2.log" \
| samtools sort \
    -@ 8 \
    -T "tmp/${sample}.longread_sort" \
    -o "validation/${sample}.longreads_to_assembly.bam" \
    - \
    2> "logs/${sample}.samtools_sort.log"

samtools index -@ 8 "validation/${sample}.longreads_to_assembly.bam"

samtools quickcheck -v "validation/${sample}.longreads_to_assembly.bam"

minimap2 --version > "validation/${sample}.minimap2.version.txt"
samtools --version > "validation/${sample}.samtools.version.txt"

The principal components are:

Component	Purpose
-a	Requests SAM alignments. Without it, minimap2 produces PAF, which is useful for some tasks but not for standard BAM-based depth and IGV workflows.
-x map-ont	Uses parameters designed for general ONT genomic reads. Substitute lr:hq only when the read accuracy and chemistry justify it.
--MD	Adds an MD tag that describes reference mismatches. It is not essential for depth but can assist later alignment inspection.
-t 16	Allocates 16 minimap2 threads. Adjust to the compute allocation actually available.
samtools sort	Produces a coordinate-sorted BAM suitable for indexing and region-based queries.
-T	Sets a temporary-file prefix. This is useful on shared systems, where default temporary locations can fill unexpectedly.
set -o pipefail	Makes the shell report failure if minimap2 fails even when samtools receives partial input and exits normally.
samtools faidx	Creates the assembly FASTA index required by many viewers and useful for extracting assembly intervals.
samtools quickcheck -v	Checks basic BAM integrity. No output indicates that no basic structural problem was found.

For high-accuracy ONT reads, the minimap2 line becomes:

minimap2 -t 16 -a -x lr:hq --MD "$asm" "$reads"

Do not select a preset solely because an assembly looks polished. The preset concerns the error profile of the reads being aligned, not the apparent quality of the reference.

First-pass checks on the alignment

bam="validation/${sample}.longreads_to_assembly.bam"

samtools flagstat "$bam" > "validation/${sample}.longreads.flagstat.txt"

samtools idxstats "$bam" > "validation/${sample}.longreads.idxstats.txt"

samtools view -H "$bam" | less

Inspect flagstat for a broadly plausible mapped fraction. A low mapping fraction may reflect contamination, an incorrect read file, a mixed isolate, a highly fragmented or incomplete assembly, or a substantial amount of host/background DNA. It is a diagnostic, not a single pass/fail threshold.

idxstats reports alignments per assembly sequence. It gives a useful first look at whether reads map across the chromosome and each candidate plasmid, but do not treat raw alignment counts as copy-number estimates: reads differ greatly in length and may have secondary or supplementary mappings.

The SAM/BAM records contain the evidence you will interrogate. Four fields are especially important:

Field	Meaning for validation
RNAME and POS	Assembly contig and leftmost mapped position.
MAPQ	The aligner’s confidence that this placement is correct relative to alternatives. Low MAPQ is common in repeats and conserved mobile elements.
CIGAR	How the read aligns: aligned bases, insertions, deletions, clipping, and sometimes split alignment structure.
FLAG	Whether an alignment is primary, secondary, supplementary, unmapped, and its strand orientation.

A supplementary alignment has FLAG 0x800; it represents another segment of the same read. Supplementary records are often biologically informative near rearrangements, integrations, circular-contig origins, and repeated elements. Do not discard them prematurely during diagnosis, even though you will usually exclude them when calculating a conservative primary-alignment depth summary.

---

Calculate depth at whole-replicon and local scales

Start with samtools coverage, which provides a compact per-reference summary including covered bases, breadth, mean depth, and mean mapping quality.

samtools coverage --no-header "$bam" \
  > "validation/${sample}.longreads.coverage.tsv"

column -t "validation/${sample}.longreads.coverage.tsv" | less -S

Interpret the chromosome and each plasmid separately. For example, a pattern such as a chromosome at approximately $80$-fold depth and a small $100$-kb plasmid at approximately $300$-fold depth can be consistent with a plasmid present at several copies per cell. Yet copy-number interpretation is only credible if the plasmid is well reconstructed, its sequence is not shared extensively with another replicon, and the reads map uniquely over most of its length.

For explicit per-position depth, use samtools depth.

samtools depth \
  -aa \
  -d 0 \
  -G 0xF04 \
  "$bam" \
  > "validation/${sample}.longreads.primary_depth.tsv"

Here:

• -aa emits every assembly position, including positions with zero qualifying coverage;
• -d 0 removes the default maximum-depth cap;
• -G 0xF04 excludes unmapped, secondary, quality-fail, duplicate, and supplementary alignments. This yields a conservative primary-alignment depth profile.

Calculate a mean depth and zero-depth proportion for every replicon:

awk '{
    sum[$1] += $3
    n[$1]++
    if ($3 == 0) zero[$1]++
}
END {
    print "contig\tlength_bp\tmean_primary_depth\tzero_depth_fraction"
    for (contig in sum) {
        printf "%s\t%d\t%.2f\t%.6f\n",
               contig, n[contig], sum[contig] / n[contig], zero[contig] / n[contig]
    }
}' "validation/${sample}.longreads.primary_depth.tsv" \
| sort -k3,3nr \
> "validation/${sample}.longreads.primary_depth_summary.tsv"

For a local region, use the same filtering logic. Suppose an AMR screen or annotation later identifies a carbapenemase region on contig_4, from 42,000 to 52,000 bp:

contig="contig_4"
start=42000
end=52000

samtools depth \
  -aa \
  -d 0 \
  -G 0xF04 \
  -r "${contig}:${start}-${end}" \
  "$bam" \
  > "validation/${sample}.${contig}_${start}_${end}.depth.tsv"

A useful summary is not merely the mean. Inspect the full depth trace and ask whether there are abrupt changes at the feature boundary or a proposed assembly join.

awk '{
    sum += $3
    n++
    if ($3 == 0) zero++
}
END {
    printf "mean_depth=%.2f\nzero_depth_fraction=%.6f\n",
           sum / n, zero / n
}' "validation/${sample}.${contig}_${start}_${end}.depth.tsv"

What depth anomalies can mean

A local anomaly should trigger investigation, not an automatic conclusion.

Pattern	Possible explanations	Next diagnostic
Sharp low-depth trough	Misassembly, poor read mapping through a divergent region, basecalling issue, or chimeric contig connection	Inspect individual reads and soft clipping in IGV.
Sustained high depth across a repeat	Collapsed repeat, multi-mapping reads, or genuine higher copy number	Inspect MAPQ, secondary alignments, and repeat annotation.
Whole plasmid at higher depth than chromosome	Multicopy plasmid	Confirm broad unique support across the plasmid rather than only shared backbone segments.
Plasmid end regions with reduced depth	Incorrect circularisation, incomplete plasmid, or repeat-linked unresolved structure	Test the end-to-start circular junction.
Carbapenemase locus with irregular depth	Transposon repeat, tandem duplication, multi-replicon sequence, or an assembly error	Examine read placements and read-level junction support.

By default, samtools depth does not count deletions in a read’s CIGAR string as coverage of deleted reference bases. The -J option includes them. For an ordinary base-support profile, keeping deletions excluded is usually sensible. For a physical question such as “does the molecule bridge this coordinate?”, a read with a deletion still spans the locus; that distinction is why depth summaries cannot replace direct inspection of alignments.

---

Identify reads that span a join

A read merely overlapping a breakpoint is not sufficient. To validate a join, define two flanks around it and require a single primary alignment to extend from the left flank, across the join, into the right flank.

Suppose the proposed join is at coordinate 500,000 on contig_4. Requiring at least 1 kb of aligned reference sequence on either side provides a more meaningful test than querying the exact base alone.

contig="contig_4"
join=500000
flank=1000

left=$((join - flank))
right=$((join + flank))

The following command retrieves alignments overlapping the window, calculates each alignment’s reference end from its CIGAR string, and retains read names only if the alignment covers both flanks.

samtools view \
  -F 0xF04 \
  -q 20 \
  "$bam" "${contig}:${left}-${right}" \
| awk -v left="$left" -v right="$right" '
function reference_length(cigar,   token,op,len,total) {
    total = 0
    while (match(cigar, /[0-9]+[MIDNSHP=X]/)) {
        token = substr(cigar, RSTART, RLENGTH)
        op = substr(token, length(token), 1)
        len = substr(token, 1, length(token) - 1) + 0

        if (op ~ /[MDN=X]/) total += len
        cigar = substr(cigar, RSTART + RLENGTH)
    }
    return total
}
{
    alignment_end = $4 + reference_length($6) - 1

    if ($4 <= left && alignment_end >= right) {
        print $1 "\t" $3 "\t" $4 "\t" alignment_end "\t" $5 "\t" $6
    }
}' \
> "validation/${sample}.${contig}.join_${join}.spanning_reads.tsv"

The output columns are:

1. read name;
2. contig;
3. alignment start;
4. alignment end;
5. mapping quality;
6. CIGAR string.

The filter has several deliberate choices:

• -q 20 retains reads whose reported probability of an incorrect placement is relatively low. It is a practical starting point, not a universal threshold.
• -F 0xF04 keeps primary, mapped, non-duplicate, non-QC-fail alignments. This prevents a read’s secondary or supplementary alignment from being mistaken for independent spanning evidence.
• The CIGAR-aware calculation matters because read length alone does not tell you how much reference sequence an alignment spans. Insertions and clipping consume read bases but not reference bases; deletions consume reference bases.

Now produce a simple list of read names:

cut -f1 "validation/${sample}.${contig}.join_${join}.spanning_reads.tsv" \
| sort -u \
> "validation/${sample}.${contig}.join_${join}.spanning_read_names.txt"

wc -l "validation/${sample}.${contig}.join_${join}.spanning_read_names.txt"

A join supported by multiple long reads that extend substantially into unique sequence on both sides is reassuring. A single read is weaker evidence, especially if it has extensive soft clipping, low MAPQ, or a highly repetitive flank. No reads spanning a join despite otherwise sufficient local coverage is a warning that requires visual inspection.

---

Inspect the evidence visually in IGV

IGV is useful because a numerical depth value can hide the structure of its underlying alignments. For each assembly you inspect, load:

• the assembly FASTA as the reference genome;
• its .fai FASTA index, created by samtools faidx;
• the coordinate-sorted BAM;
• the BAM index, normally a .bai file.

Navigate to the join window, for example contig_4:499000-501000. The named reads in your spanning-read table can be searched individually.

A credible long-read bridge should show:

1. One alignment extending well across both flanks, not two disconnected fragments that merely approach the join.
2. Adequate MAPQ, unless you have established that the region is inherently repetitive.
3. Limited soft clipping at the join. A modest amount is normal for ONT reads; a cluster of reads clipped exactly at the same coordinate is suspicious.
4. Several distinct reads, ideally from both strand orientations where coverage permits.
5. No obvious competing placements that would make the observed bridge ambiguous.
6. Reasonably stable local depth on both flanks.

Examine supplementary alignments separately when the assembly is structurally complex:

samtools view \
  -f 0x800 \
  "$bam" "${contig}:${left}-${right}" \
| cut -f1 \
| sort -u \
> "validation/${sample}.${contig}.join_${join}.supplementary_read_names.txt"

Supplementary records do not automatically invalidate a join. They may reflect a genuine plasmid rearrangement, a chromosomal integration, a repeated transposon, or a read spanning the linearised origin of a circular contig. They do mean that the simple “one primary alignment bridges the join” rule is insufficient by itself.

---

Test the circular contig origin explicitly

Circular bacterial replicons are represented in FASTA as linear strings. The biological adjacency between the last base and the first base is therefore invisible as an internal coordinate in the standard contig. Reads crossing that origin often appear as split or supplementary alignments near opposite contig ends.

Create an artificial reference containing the final 10 kb followed directly by the first 10 kb. Any genuine origin-spanning long read should then align continuously across the central artificial junction.

contig="contig_4"
window=10000

length=$(awk -v c="$contig" '$1 == c {print $2}' "${asm}.fai")

{
    echo ">${contig}_origin_junction"
    samtools faidx "$asm" \
        "${contig}:$((length - window + 1))-${length}" \
        "${contig}:1-${window}" \
    | awk '!/^>/ {printf "%s", $0} END {print ""}'
} > "validation/${sample}.${contig}.origin_junction.fa"

Map the long reads to this short junction reference:

minimap2 \
  -t 16 \
  -a \
  -x map-ont \
  "validation/${sample}.${contig}.origin_junction.fa" \
  "$reads" \
| samtools sort \
    -@ 8 \
    -o "validation/${sample}.${contig}.origin_junction.bam" \
    -

samtools index "validation/${sample}.${contig}.origin_junction.bam"

The synthetic junction lies at coordinate 10,000. Apply the same spanning-read method with flanks around this central point. When interpreting the result, require reads to extend into sufficiently unique sequence on both sides. A repeat located at the natural contig origin can make apparent support ambiguous.

---

A defensible interpretation record

For every join or resistance region examined, record the evidence in a compact validation table. This avoids later overstatement when comparing environmental and patient plasmids.

Field	Example entry
Assembly and contig	ISO001, contig_4
Feature or join	blaKPC region or Autocycler join at 500,000
Long-read input	Filtered ONT reads, same-run or held-out
Mapping preset	map-ont
Local primary mean depth	74.2-fold
Breadth and zero-depth positions	100% breadth; 0 zero-depth bases
Number of stringent spanning reads	11 reads, MAPQ at least 20, 1-kb flanks
Competing/supplementary evidence	Two supplementary reads in adjacent repeat
Interpretation	Join supported, but repeat-associated ambiguity retained as a caution

The calibrated conclusion is not “the plasmid is proven correct.” It is: the proposed adjacency is supported by a specified number of long-read molecules under documented alignment and filtering criteria. That wording remains valid when you later integrate plasmid reconstruction, mobile-element annotation, and cross-isolate comparisons.

---

Key takeaways

Long-read remapping turns an assembly into a testable model. Build a sorted and indexed BAM with minimap2 and samtools, preserve logs and versions, and inspect basic mapping summaries before interpreting depth.

Use samtools coverage for compact per-replicon summaries and samtools depth for explicit base-level or local depth. Interpret anomalies in the context of repeats, multi-mapping, plasmid copy number, and assembly uncertainty rather than treating depth as a simple quality score.

Finally, validate a proposed join with reads that span substantial flanks on both sides, then inspect those reads visually. For circular replicons, create an artificial end-to-start junction sequence and test it directly.

The next lesson moves from read-level support to assembly-level acceptance: using QUAST and CheckM2 to assess contiguity, completeness, contamination, and whether a reconstructed genome meets a documented quality rubric.