Map reads

Put sequencing reads (DNA or RNA, short or long, single or paired) or assembled contigs where they came from on a reference, and get the alignment back as SAM, with every read's position, strand, mapping quality and alignment.

What it doesCopy Anchor Link

Takes a file of reads and a reference, and finds where each read belongs. The method is minimap2: short exact matches (minimizer seeds) between read and reference are found first, colinear seeds are chained into candidate placements, each candidate is aligned base by base, and every read is given a mapping quality that says how sure its placement is. The output is checked record for record against minimap2 2.31, so the SAM it writes is the one minimap2 would.

Pick what the reads are, and every setting takes minimap2's values for that data:

RNA is aligned across introns: a gap in the read's alignment long enough to be an intron is written as one (N in the CIGAR), and its two ends are scored by the splice signal they carry, GT..AG best, then GC..AG and AT..AC, anything else at a cost you can set. Each read is aligned twice, once assuming its transcript lies on each strand of the genome, and the better alignment decides the transcript strand (the ts:A tag). The result lists every splice junction with the number of reads that support it; clicking one selects that intron on the map.

Long reads and contigs are chained differently from short ones: when the best chain leaves much of a read unexplained, its seeds are chained again in a much wider band, so a read spanning a large deletion stays one alignment. Where an alignment stalls, the stretch is tested for an inversion (the read running backwards against the reference), and an inverted segment is written as its own record.

Give a second file of mates and the reads are mapped as pairs: the two files pair record for record, the mates are expected to face each other, and a pair whose mates land close enough to have come from one DNA fragment is marked proper. Pairing does more than label: when one mate sits in a repeat and its partner maps uniquely, the partner's placement lifts the repeated mate's mapping quality.

The output is SAM (each read with its alignment) or PAF (one line per hit, minimap2's default, and the usual input to assemblers and overlap tools); PAF can skip base-by-base alignment, which is much faster and leaves each hit's ends approximate.

The alignment also downloads as a sorted BAM with its index (.bam + .bai), the form IGV and the samtools family read directly.

Open alignment writes the SAM into the project and opens it in the pileup against the same reference, so the reads can be read in place.

When you'd use thisCopy Anchor Link

• You have raw reads and a reference. Mapping is the step between a FASTQ and everything downstream: the pileup, a variant call, a consensus.
• "Is my plasmid what I think it is?" Map whole-plasmid or amplicon reads onto the design and read the pileup at the junctions.
• "How much of this library is my target?" The share of reads that map, and where the unmapped ones come from, is a direct read-out of contamination or off-target sequence.

Either side of it: Alignment to look at the result, and Vector contamination when unmapped reads look like vector.

What this computesCopy Anchor Link

Before any of that — read QC. Adapters and low-quality ends are trimmed first, and each is reported separately:

• Bases removed by quality trimming, and bases removed by adapter removal, counted apart rather than as one total. Losing a large share of your bases to quality points at the sequencing run; losing it to adapter points at a DNA insert shorter than your read length. They are different problems and one number cannot tell them apart.
• How many reads carried the adapter, and what share of the run that is.
• Mean read length before and after, and how many reads were trimmed to nothing at all or left shorter than a length you name. Reads are counted, never dropped: nothing disappears silently.

Quality trimming is not a per-base threshold, which is the most common misunderstanding of it. Cutting at the first base below a cutoff would truncate a read on a single bad base in an otherwise good tail. Instead the cutoff is subtracted from each quality and the running total accumulated from the end inward, and the read is cut where that total peaks — so a lone dip is absorbed by the good bases after it, while a sustained decline is not. This is the rule BWA uses.

The adapter is found by alignment, not by search. It is sequenced with the same error rate as everything else, and a read usually ends part way through it, so an exact-match search misses the commonest case. A partial adapter at the very end of a read, with a sequencing error in it, is still found — and how many errors were allowed is a number you set.

Quality trimming runs first and adapter removal second. That order changes the answer whenever a low-quality tail is also the start of an adapter, and it is the order the method below uses.

• Placement of every read: reference, position, strand, and the alignment itself (CIGAR), including soft-clipped ends where the read runs past what matches.
• Mapping quality (MAPQ) for every placement, 0 to 60. MAPQ 0 means the read aligns equally well elsewhere, so its position is a guess.
• Edit distance (NM), the alignment score, and the peak local score of each alignment.
• Proper pairs: both mates mapped, facing each other, within the fragment-length limit.
• Secondary placements: other places a read aligns nearly as well (repeats). Counted always, written to the SAM on request.
• Supplementary alignments: a read split into pieces that map to different places, which marks a junction the reference does not have.
• Unmapped reads, still written to the SAM so no read is silently dropped.
• Repetitive seed length: how much of each read was covered by seeds too frequent in the reference to use; it lowers MAPQ.
• Splice junctions (RNA presets): each intron's start, end, length, the transcript strand, a motif score for its two ends (6 = canonical GT..AG) and the number of reads supporting it, counted over primary alignments of MAPQ 10 or more; how many reads were spliced; how many junctions are canonical and how many rest on a single read. Downloadable as the BED-like junction file minimap2 writes.
• Transcript strand (RNA presets): which strand of the genome the read's transcript lies on, decided by which orientation of its splice signals aligns better.
• Overlaps (all-vs-all presets): every pair of reads that overlap, each pair once, with the overlap's coordinates on both reads, its strand, its matching bases and length; how many reads overlap at least one other and how many overlap none.
• Inversions (long reads and contigs): an inverted segment inside a read, written as its own record on the opposite strand, with the lower MAPQ of the two alignments flanking it.
• Divergence estimate (long reads and contigs): per alignment, estimated from the share of the read's seeds the alignment kept; a placement kept only because it is on the other strand is dropped when it is far more divergent than the best one.
• samtools flagstat: the standard QC table of the alignment: records in total, primary, secondary and supplementary, duplicates, mapped, paired, properly paired, singletons and mates on another reference, each split into QC-passed and QC-failed, counted exactly as samtools counts them.
• samtools coverage: per reference, the reads that pass samtools' filters, the bases covered, the percent covered, the mean depth, the mean base quality and the mean mapping quality.
• samtools depth (when asked for): the read depth at every position, either every position a read spans or every position of the reference (-a), with samtools' own filters: count deletions (-J), count the overlap of a read pair once (-s), and minimum base quality, mapping quality and aligned read length. The run reports how many positions sit at depth 0, the deepest position and the mean depth, and the full table downloads as samtools prints it.
• samtools markdup (when asked for): duplicate marking after samtools fixmate -m and a coordinate sort. Reads whose outer ends, strands and mate positions match are copies of one fragment; the one with the highest base-quality sum is kept and the rest are flagged 0x400 (or removed). Duplicates can be keyed on the template's two ends (samtools' default, as Picard) or on each read's start, and with an optical distance the duplicates whose Illumina tile x/y lies within that many pixels of the original are counted as optical. The run reports duplicate pairs and single reads, optical duplicates and the estimated library size (Picard's estimator), and flagstat, coverage and depth then count the flags.
• MD tags (when asked for): per alignment, the reference base at every mismatch and the bases a deletion removed, so a later step can read each difference without the reference.
• samtools consensus (when asked for): the sample's own sequence, called column by column from the reads over the reference. By default this is samtools' Bayesian model from the Gap5 assembler. At each column it weighs the likelihood of every base, and of every two-base mixture at a heterozygous site, given each read's base quality and mapping quality. A read's mapping quality is lowered near its own mismatches, indels and soft clips (read from its MD tags), and a base's quality is lowered beside low-quality neighbours. A call below the quality cutoff is written as N. The other callers are precise and blended variants of the model, samtools 1.16's model, and simple counting, where a base must carry a set fraction of the column and a second base makes the call heterozygous once it reaches a set fraction of the first. Platform profiles (Illumina HiSeq, PacBio HiFi, Nanopore R10.4 simplex and duplex, Ultima) bring their own quality calibration and priors. The output is FASTA, FASTQ with a Phred quality for every base, or the pileup behind every call. Heterozygous calls can be written as IUPAC ambiguity codes. Deletions can be shown as * or left out, and inserted bases kept or dropped. The consensus covers the span the reads cover or whole references, with uncovered bases written as N or taken from the reference. The run reports the sequences and bases written, the N count, the ambiguity codes, the lower-case calls (a base or a deletion), the deletions shown and the mean quality of the called bases. Every samtools consensus option can be set, and the result is checked byte for byte against samtools 1.24. Open consensus adds the sequence to the project and opens it.
• Variants (when asked for): SNVs and indels called from the alignments, as bcftools mpileup | bcftools call -mv calls them, written as VCF. At every covered position the reads' bases are piled up (with a depth cap, the overlap of two mates counted once, orphan and flagged reads skipped) and each base's quality is capped by its neighbours', by the read's mapping quality and, by default, by its base alignment quality (BAQ): the probability that the base is aligned to the wrong reference position, from a banded pair-HMM realignment. The error model turns those qualities into genotype likelihoods (the PL field), and the multiallelic caller picks the allele set, the genotype (GT) and the site's QUAL under a mutation-rate prior. Indels are scored by realigning each read against a consensus of the sample carrying each candidate indel; an indel's quality is capped by its length through the gap-open and gap-extension error probabilities and, in a homopolymer run, by the homopolymer error coefficient, and a candidate needs a minimum number and fraction of reads carrying it. Each call carries its allelic depths (AD), the high-quality strand counts for REF and ALT (DP4), mean mapping quality, the fraction of reads at mapping quality 0, the Mann–Whitney bias z-scores for read position, mapping quality, base quality, strand versus mapping quality and soft-clip length, the variant distance bias and the segregation-based metric; indels also carry the largest number and fraction of reads supporting one. The run reports SNVs and indels, homozygous and heterozygous calls, and transitions to transversions among the SNVs; each of the strongest calls selects its bases on the map. Every call is checked record for record against bcftools 1.24. Open variants opens them in the variant viewer.
• Summary: reads mapped, proper pairs, MAPQ distribution, mean identity of primary alignments, and per reference the reads placed and the share of its length covered.

Every setting is one of minimap2's own options under its own flag (seed size and window, homopolymer compression, the seed frequency cut-off, the chaining band and the wider long-join band, the longest gap a chain may cross, chaining score, match, mismatch and transition scores, the two gap costs, Z-drop and the inversion Z-drop, the non-canonical splice-site cost and the longest intron, end bonus, secondary-alignment ratio and count, and the fragment length for pairs), starting at the chosen preset's values. A changed setting is named in the result.

Read trimmingCopy Anchor Link

Trimming is cutadapt's algorithm (Martin 2011), with quality trimming following the rule from BWA. Correctness here is provenanced by authority rather than by our own testing: the trim positions, the adapter alignments and the trimmed reads are each checked against cutadapt's own output over every read of the sequencing files bundled with this app — 94,200 trim positions, 376,800 adapter alignments, and the trimmed FASTQ compared read for read against what cutadapt writes from the command line.

It trims, and does not filter. A read trimmed to nothing is returned empty and counted, never removed, so the number of reads going into mapping is the number that came out of the sequencer.

ControlsCopy Anchor Link

Reference is what the reads are placed on, Reads the FASTQ or FASTA to place (FASTQ keeps base qualities in the alignment), and Mate 2 the second file of a pair, which must pair record for record. Leaving it on None (single-end) maps single-end.

What it does not doCopy Anchor Link

Splice junctions are found from the reads alone: an annotation of known junctions cannot be given to bias them. Read pairs are mapped by the two short-read presets (DNA and RNA); the long-read and assembly presets map every read on its own. Coverage counts the aligned bases of each read, so an intron a read spans is not counted as covered. The consensus takes one region at a time (no region file), and its quality calibrations are the built-in platform tables rather than a table read from a file. Variant calling is bcftools' default model for one diploid sample: its experimental indel models (--indels-2.0, --indels-cns and the platform profiles that switch them on) are not offered, and neither is a depth cap above 255, where bcftools downsamples the reads at random.