Low-mappability regions (SMN)
TL;DR: a pileup looks the same whether its reads belong at a locus or merely landed there. Four tracks tell the difference, and genomes.jbrowse.org publishes all of them for hg38, so this page is a click-path.
Prerequisites
- nothing to install to read along: every track comes from the hosted hg38 config, and the one read track is a public CRAM added as a session track
- to re-measure the numbers on this page, the tools listed under Reproduce it end to end
Where the data comes from
Every hosted lane here is a file already wired into the genomes.jbrowse.org hg38 config, plus the 1000 Genomes high-coverage short-read CRAM and its ONT long-read release (Gustafson et al. 2024).
- Umap k100 multi-read mappability: https://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k100.Umap.MultiTrackMappability.bw
- gnomAD v3 mean genome coverage: https://hgdownload.soe.ucsc.edu/gbdb/hg38/gnomAD/coverage/v3-genome/gnomad.coverage.mean.bw
- GIAB's low-mappability and segmental-duplication regions: https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/GIAB/alllowmapandsegdupregions.bb
- ENCODE's blacklist: https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/encBlacklist.bb
- the GRC's exclusion list: https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/grcExclusions.bb
- UCSC's own problematic-regions comments: https://hgdownload.soe.ucsc.edu/gbdb/hg38/problematic/comments.bb
- the DGV merged CNV catalogue: https://hgdownload.soe.ucsc.edu/gbdb/hg38/dgv/dgvMerged.bb
- the 1000 Genomes ONT long-read SV callset: https://hgdownload.soe.ucsc.edu/gbdb/hg38/lrSv/1kgOnt.bb
- NA12878 at 30x, GRCh38, the read track added to the session: https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram
- UCSC's hg38-to-CHM13 liftOver chain set: https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz
- GM18501 ONT long reads aligned to GRCh38, counted rather than shown since the bucket serves no CORS headers: https://s3.amazonaws.com/1000g-ont/PROCESSED_DATA/ALIGNED_TO_HG38/MINIMAP2_ALIGNED_BAMS/GM18501-ONT-hg38-R9-LSK110-guppy-sup-5mC.phased.bam
- the same sample aligned to T2T-CHM13: https://s3.amazonaws.com/1000g-ont/PROCESSED_DATA/ALIGNED_TO_CHM13/MINIMAP2_ALIGNED_BAMS/GM18501-ONT-chm13-R9-LSK110-guppy-sup-5mC.phased.bam
The SMN1 and SMN2 duplication
SMN1 and SMN2 sit about 900 kb apart on chromosome 5 and are roughly 99.9% identical across their ~28 kb. Which of the two a read came from is the clinically interesting question, since it is the copy number of SMN1 that spinal muscular atrophy turns on, and it is also the question a 150 bp read cannot answer: the same sequence exists twice, so an aligner given a read from either copy has two equally good places to put it.
An aligner reports that as MAPQ 0. The read is still aligned and still drawn
where it aligned; MAPQ is -10 log10 Pr{mapping position is wrong}
(SAM specification), so 0 is
the aligner saying the position it chose is about as likely wrong as right.
The block, and the reads inside it
The affected sequence is a much larger block than the gene. The two published
annotations disagree about where it ends: GIAB's interval stops well short of
where ENCODE's blacklist continues to. scan_mappability_qc.sh bins the
coverage lane so a locus between the two edges can be settled by measurement,
and the lane stays low across the span GIAB has let go of.
The block is one interval. GIAB's annotation over this arm is a single megabase-and-a-half region, a second one a few kilobases past it, and then nothing larger than a few kilobases for megabases in either direction. Short reads fail across a whole gene neighbourhood here.
The lower panel is the same block at the scale a read lives at, where reads do not recover until well past the end of SMN1.
The same block in T2T-CHM13
T2T-CHM13 is a finished assembly of this chromosome. UCSC publishes an hg38 to CHM13 liftOver chain set, and over this block that chain set does not resolve to one correspondence:
tabix https://jbrowse.org/ucsc/hg38/liftOver/hg38ToHs1.over.pif.gz \
tchr5:69200000-71700000Several of the chains it returns are long, they overlap each other on both sides, and some of them run backwards.
The gene order is the same in both assemblies (SMN2 first, then SMN1), so this is two copies similar enough that a whole-genome chainer can join either one to either one, which is what the Umap and MAPQ lanes below say per base. The array is a different length in the two assemblies, the genes sitting closer together in CHM13.
The 1000 Genomes ONT release, the same project as the long-read SV callset in
the wide figure, aligned some of its samples to both references with the same
minimap2 pipeline, so one sample answers the question twice.
scan_mappability_qc.sh counts GM18501's records over SMN1 in each assembly's
own coordinates:
| reference | records | MAPQ 0 | MAPQ 60 |
|---|---|---|---|
| GRCh38 | 290 | 46.6% | 6.9% |
| T2T-CHM13 | 290 | 46.6% | 9.7% |
The long reads place better than the short-read lane does at the same gene, and a large share of them still fit somewhere else as well. Between the two references the columns barely move.
What the lanes are
The lanes are independent of each other:
- Umap k100 multi-read mappability is computed from the reference alone. For each position it gives the fraction of overlapping 100-mers that are unique in the genome. Positions where no 100-mer is unique are absent from the file rather than stored as zero, so the lane goes blank rather than to the floor. Most of the genome scores near the top of its range, so a blank stretch is unusual. How it summarizes decides whether it survives a wide window: the default Score → Summary score mode → Whiskers draws each pixel's min and max, so a bin touching one unique position paints full height, while Minimum takes the worst position in the bin and sits on the floor across the block, stepping up at the same coordinate as the MAPQ 0 to MAPQ 60 transition and the gnomAD coverage step. Past about a kilobase per pixel even Minimum saturates low, so this lane belongs in the narrower frame.
- gnomAD v3 mean genome coverage is the outcome of that annotation on real data, averaged over tens of thousands of sequenced genomes. gnomAD drops non-uniquely-placed reads before computing it, so wherever the lane above is blank this one falls.
- Mapping quality on the reads is the aligner's own account, one read at a time, in the sample on screen. Red is MAPQ 0, meaning the aligner found another place the read fits equally well; yellow is MAPQ 60 and above.
- The GIAB low-mappability + segdup lane in the wide panel is a published opinion of the same sequence, drawn by a project that had to decide where its benchmark regions stop.
Everything on this page except the read track comes out of the hosted hg38 config at genomes.jbrowse.org: find them in the track selector under Multi-read mappability, gnomAD v3 Genome Coverage and Problematic Regions. The reads are the public 1000 Genomes NA12878 high-coverage CRAM, added to the session; the figure's link opens both together.
Depth at the locus and at a control
scan_mappability_qc.sh counts the reads in equal windows over SMN1 and over
the right-hand end of the frame, from the same library, and they come back at
the same depth: without a MAPQ filter a coverage track draws flat across both.
What separates them is the share of those reads sitting at MAPQ 0, which is most
of them at SMN1 and almost none at the control.
The gnomAD lane is what a MAPQ filter does to a depth track, in a different set of samples. It drops to a fraction of the control's depth over SMN1, because MAPQ 0 reads were dropped before the average was taken.
SV calls over the block
The long-read lane in the wide figure is empty across the block. Counting over
the flagged block and an equal-width window on either side of it,
scan_mappability_qc.sh finds the callset nearly silent inside and populated on
both sides, where the older DGV catalogue carries records throughout.
Widen the same count to the whole chromosome and both catalogues put a larger share of their calls inside the flagged regions than those regions' share of chr5 would predict. Segmental duplications are copy-number variable, so this is where real variation lives as well as where artifacts do.
At this locus, in this sample, the reads carry no information about which copy they came from, so a short-read call over it cannot be checked against them.
Checking your own locus
The same three tracks and a control work anywhere in hg38:
- Open the hosted hg38 config and turn on Umap M100, gnomAD v3 Genome Coverage - Mean Coverage, and the GIAB Problematic Regions and Problematic Regions annotation tracks.
- Add your reads and set Color by... → Mapping quality from the track menu. Turn on Show legend in the same menu.
- Take a second window of the same width, from the same sample, outside every flagged interval, and put the two side by side.
The same comparison is three counts per window, -q being a minimum MAPQ:
samtools view -c "$CRAM" chr5:70,900,000-71,000,000 # every read
samtools view -c -q 1 "$CRAM" chr5:70,900,000-71,000,000 # placed at all
samtools view -c -q 60 "$CRAM" chr5:70,900,000-71,000,000 # placed uniquelyRun it on the control window too, at the same width and from the same sample.
The eye reads the Umap lane as high or low, and the number behind it is one command over the same bigWig the lane draws:
# A position where NO 100-mer maps uniquely is ABSENT from the file rather than
# stored as zero, so bigWigToBedGraph emits no interval there at all. That makes
# the fraction of the span carrying any value the number to read, and it is
# stronger than the mean: a mean is taken only over what is present.
bigWigToBedGraph -chrom=chr5 -start=70049000 -end=70077000 \
k100.Umap.MultiTrackMappability.bw stdout |
awk -v s=70049000 -v e=70077000 \
'{cov += $3 - $2} END {printf "%.1f%% has a value\n", 100 * cov / (e - s)}'Run it over the control window too. A percentage means nothing until an ordinary gene has produced one.
Reproduce it end to end
Every number on this page comes from
scripts/scan_mappability_qc.sh,
run against the same files the figures draw. It needs kent tools (bigWigInfo,
bigWigToBedGraph, bigBedToBed), bedtools, samtools, curl and awk,
downloads the four small annotation files it reads twice, and streams the rest.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/scan_mappability_qc.sh
bash scan_mappability_qc.shIt prints the mappability, coverage, region-annotation, MAPQ and callset
sections in the order this page uses them, so a locus swapped into its LOCI
list is measured the same way.
See also
References
- Li H, Handsaker B, Wysoker A, et al.
The Sequence Alignment/Map format and SAMtools.
Bioinformatics 25:2078-2079 (2009), and the current
SAM specification, which
defines MAPQ as
-10 log10 Pr{mapping position is wrong}. - Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Research 18:1851-1858 (2008), which introduced that estimator.
- Karimzadeh M, Ernst C, Kundaje A, Hoffman MM. Umap and Bismap: quantifying genome and methylome mappability. Nucleic Acids Research 46:e120 (2018), the source of the k100 mappability track.
- Gustafson JA, Gibson SB, Damaraju N, et al. High-coverage nanopore sequencing of samples from the 1000 Genomes Project to build a comprehensive catalog of human genetic variation. Genome Research 34:2061-2073 (2024), the source of both the long-read SV callset in the wide figure and the GRCh38 / T2T-CHM13 read counts above.
Feedback on this tutorial is welcome: contact us.