How to Check a Whole-Plasmid Sequencing Result

What arrives, and which file is evidence

Plasmidsaurus's own technical documentation for whole-plasmid sequencing describes the pipeline: Oxford Nanopore sequencing on R10.4.1 flow cells with v14 library chemistry and the latest super-accurate basecalling model, after a "linearization of the circular input DNA in a sequence-independent manner". The delivered set is a consensus sequence as .fasta and .gbk, an interactive plasmid map as .html, a read-length histogram, a virtual gel, a chromatogram as .ab1, a coverage plot, per-base data as .txt and .tsv, and the raw reads as .fastq.gz, available separately.

That list has one observation in it and several derivations. The .fastq.gz is what the sequencer read. Everything else — including the consensus, the map and the chromatogram — is computed from it. Other providers deliver a similar mix, and the same distinction applies to all of them.

The consensus is a very good derivation. The same documentation states that "Coverage over ~20x indicates a very accurate consensus" and that "Consensus accuracy is often above Q60, which corresponds to 99.9999%, or one error per 1,000,000 bases", off runs that "typically yield in the high dozens to hundreds (or thousands!) of raw sequencing reads". It also says plainly: "We do not guarantee any specific level of coverage, as the number of raw reads generated can vary substantially depending on sample quality."

Q60 on a 5 kb plasmid means the expected number of consensus errors is well under one. It does not mean the consensus answers your question, for two reasons covered below.

Start by comparing the consensus with your design

This is the cheapest useful check and it catches most real problems: assembly junctions that closed wrong, an insert in the wrong orientation, a fragment that dropped out, a mutation carried in from a PCR step.

Align the consensus FASTA against your intended sequence and read the difference list. SeqBench's Sanger vs Reference tool takes a reference plus a read as raw text or FASTA — despite the name it does not need an .ab1 — and reports identity plus every mismatch, insertion and deletion, so it works directly on a consensus. For a pure two-sequence comparison the pairwise alignment tool does the same job with the alignment shown.

One parameter to watch: that tool's minCoverage defaults to 0.9, meaning a read spanning less than 90% of the reference returns "ambiguous_low_coverage" regardless of identity. That default is right for a Sanger read against a whole plasmid and wrong when your reference is deliberately just the junction you care about — lower it in that case rather than reading the verdict as a failure.

If the consensus matches your design end to end, you are nearly done. If it does not, the difference list tells you where to look, and the next question is whether each difference is real.

Why the consensus is not the end of it

Two failure modes survive a perfect-looking consensus.

The first is a mixed population. A consensus is a majority vote. If your tube contains two plasmid species — a correct clone and a deletion variant, or a colony that was not clonal — the consensus reports the winner and says nothing about the loser. A 70/30 mixture and a pure prep produce consensus files that look identical.

The second is that the consensus inherits the systematic errors of the reads it was built from. Nanopore's residual error is not uniformly distributed; it concentrates in homopolymers. R10.4.1 chemistry reads at roughly 99% modal single-read accuracy — the Q20+ generation — and duplex basecalling reaches about Q30, but reported basecalling accuracy still falls off as homopolymer length grows, with a marked drop past about 10 bases. An error that every read makes in the same place is not fixed by more reads.

Both failure modes are visible in the raw reads and invisible in the consensus. That is the whole reason the .fastq.gz is worth downloading.

A worked example on pUC19

To show what read-level checking actually returns, here is a run against pUC19 (GenBank L09137, 2,686 bp) with eight constructed reads of 687 to 1,601 nt. The reads are constructed rather than real, and deliberately much cleaner than nanopore reads — the point is to show what the output looks like when you know the right answer, not to simulate an instrument. One of the eight straddles the reference's arbitrary start, which is its own problem and gets a section below.

The read set carries three things: a G to A substitution at position 2,124 present in six of the eight reads; one unique substitution in a single read; and a one-base deletion inside a seven-T homopolymer at 1,407-1,413, present in two reads.

SeqBench's read-back verification tool aligns every read to the reference with the real minimap2 binary (Li, Bioinformatics 2018;34:3094-3100, doi:10.1093/bioinformatics/bty191) and separates positions by how many reads agree. All eight reads mapped, at 99.86% to 100% identity. The output:

• Position 2,124, G>A — 6 supporting reads. This is the one to act on. Deep annotation puts it inside the AmpR (bla) coding sequence, which runs 1,626-2,486 on the minus strand; translating that codon makes it TGC to TGT at codon 121, cysteine to cysteine. A silent change, which is consistent with the clone having grown on ampicillin at all.
• Position 1,413, a one-base deletion — 2 supporting reads. Inside the seven-T homopolymer. Two reads is enough to clear the default threshold and be reported as consensus-supported, and it is exactly the artefact class nanopore is known for.
• One further position reported as single-read-only, and set aside: read 3's g.1780C>T, seen once and by nothing else.

The support threshold is a real decision, not a default

The tool's minSupportingReads defaults to 2, and this example shows why you should not leave it there on real data. The corrected consensus it returns — the reference with every consensus-supported edit applied — came back at 2,685 bp rather than 2,686. The two-read homopolymer deletion was promoted into it, and the plasmid lost a base it never lost.

The right threshold scales with your depth. On a run with hundreds of reads, a variant supported by two of them is present in around 1% of the population, which is a basecalling artefact far more often than it is a subclone. There is no published rule for where to put the cut, but something on the order of a fifth of your read count separates "the clone carries this" from "a few reads made the same mistake", which at these error rates they will. Set it deliberately and say what you used, rather than inheriting a default that was chosen to be safe on three Sanger reads.

And treat indels inside homopolymers as suspect regardless of what the threshold says. A run of seven or more identical bases is the one place where a systematic error can accumulate consistent support across many reads. If a length change in a homopolymer changes your answer — a frameshift call, a promoter spacing — confirm it another way rather than trusting depth alone.

The inverse case matters too. Single-read differences are the null hypothesis, not a warning. At 99% single-read accuracy a 2,686-base read is expected to carry about 27 differences from a perfect reference by chance alone. On a real run you should see a long list of single-read positions, and a clean report is the one where none of them recur.

The seam: the failure mode that looks exactly like a pass

A FASTA reference is linear and its start position is arbitrary — whatever the depositor or your cloning software picked. Your plasmid is not linear, and whole-plasmid protocols cut it at a sequence-independent position, so the reads begin all over the molecule and most of them cross that arbitrary join. An aligner working against a linear string cannot place a read that runs off the end and back to the beginning: it keeps the better-scoring half and soft-clips the rest.

That is not a cosmetic problem, and it is worth seeing what it costs. Here is a second run: a single G>A substitution at position 150 of pUC19, carried by four reads that each straddle the join — they start between positions 2,280 and 2,350, run off the end, and continue to somewhere between 260 and 340. Two ordinary reads elsewhere carry nothing. The plasmid genuinely has a mutation, in four of six reads.

Aligned against the linear reference, all four seam reads come back mapped at 100% identity, each running from wherever it starts to the last base of the reference, with 54-57% query coverage and 260-340 bases clipped. The list of consensus variants is empty. The corrected consensus is 2,686 bp — unchanged. The mutation is in the clipped remainder of every read that covers it, so it was never called, and the hard "no consensus variants" check passes with the message "No position had multiple reads disagreeing with the reference in the same way." The overall verdict is a pass on a plasmid that carries a mutation.

Rerun the identical input with the circular option on and all four align end to end: 100% query coverage, nothing clipped, spans now reported with the end coordinate before the start — 2,300-300 for the first of them — plus a wrapsOrigin flag. And g.150G>A appears with four supporting reads. The hard check now fails, correctly.

So the practical rule is short: turn circular on for plasmid data. It is off by default because a reference is not necessarily circular and the tool will not assume; the mechanism is that the reference is laid end to end with itself so a straddling read has one contiguous place to sit, with coordinates folded back into the real range afterwards.

If you forget, the tool now says so rather than leaving you to notice. In the linear run above the end-to-end check fires with: "4 of 6 mapped read(s) had a meaningful share of the READ clipped out of the alignment — their identityPct describes only the aligned part, and any variant in the clipped remainder was never called ... re-run with circular: true to align through the join". That is a soft check, so the overall verdict still reads pass — which is exactly why it is worth reading the checks and not just the verdict.

The same reasoning applies to any tool you use for this, most of which do not offer a circular mode at all. Compare each read's aligned length against its own length; a pile of reads clipped against the two ends of the reference is a seam, not a deletion, and a region that looks uncovered at position 1 usually is not.

What whole-plasmid sequencing catches that Sanger does not

The variant in the example sits at position 2,124, inside the AmpR coding sequence — backbone, not insert. No primer pointed at your cloning junctions would have covered it, and the classic verification protocol of two reads across the insert would have returned a clean result on a plasmid that carries a change.

That is the actual argument for the technique. Sanger tells you the region you asked about is right. Whole-plasmid sequencing tells you about the kilobases you never asked about — the resistance marker, the origin, the region a previous PCR step passed through — and those are where the surprises accumulate in constructs that have been through several rounds of cloning.

It also changes what a negative result means. "The insert sequence is correct" and "the plasmid is correct" are different claims, and only one of them survives a backbone mutation.

A checklist

1. Check the coverage plot and read count first. A consensus built from a handful of reads is a different object from one built from hundreds, and the provider does not guarantee depth.
2. Align the consensus against your intended sequence and read the difference list. If it is empty and the coverage is good, stop here.
3. For every difference, go back to the raw reads and count support. A difference in the consensus that only a minority of reads carry is a different finding from one that all of them carry.
4. Raise the support threshold to match your depth instead of accepting a default of 2.
5. Treat indels in homopolymers of seven or more bases as unproven until confirmed another way.
6. Turn on circular mode for a plasmid, and compare each read's aligned length against its own length. Reads clipped against the two ends of the reference are crossing the seam, and a variant in the clipped part is not reported at all.
7. Translate any coding change before panicking or relaxing. Silent, missense and nonsense are three different results, and a construct linter will also flag a premature stop, an unwanted restriction site or a repeat introduced along with it.
8. Keep the .fastq.gz. The consensus can be rebuilt from the reads; the reads cannot be rebuilt from the consensus, and they are the only file that can settle an ambiguous call later.

Related references

Related tools

Related guides

• How Long Should a Gibson Assembly Overlap Be? →
• How to Find an Enzyme That Cuts Your Plasmid Exactly Once →
• How Allele-Specific Primers Actually Discriminate (KASP and ARMS) →