# A deletion from a BAM's own fields (NA12878, insert size) We find a deletion in a genome from the reads alone. A paired-end read pair that straddles a missing stretch of chromosome maps with its two ends further apart than the library was made, so its insert size gives the deletion away before any variant caller runs; and where one chromosome copy lacks the stretch, half the reads are gone with it. We plot those two numbers, the depth and the insert size, as marks declared over the alignment file, first across one deletion and then along a whole chromosome. The track is a `LinearMarkDisplay`: a read's fields become channels of a plot with a JSON entry, the way a BED column does in the [Alu tutorial](https://jbrowse.org/jb2/docs/tutorials/alu_age). ## Prerequisites - a JBrowse to open the figures' sessions in ([Web](https://jbrowse.org/jb2/docs/quickstart_web) or [Desktop](https://jbrowse.org/jb2/docs/quickstart_desktop)); every file here is a URL, so nothing needs hosting to read along - [samtools](https://www.htslib.org/) and htslib (`bgzip`, `tabix`), for cutting the pairs out of the file and for checking a window by hand - [bcftools](https://www.htslib.org/), for reading the callset at the end - [Node.js](https://nodejs.org/) and the [JBrowse CLI](https://jbrowse.org/jb2/docs/cli), for the build script ## Where the data comes from NA12878's reads are from the 1000 Genomes Project's high-coverage release ([Byrska-Bishop et al. 2022](https://doi.org/10.1016/j.cell.2022.08.004)), aligned to GRCh38, and the deletions checked against at the end are the same release's structural-variant callset. - the reads, 30x Illumina paired-end as a CRAM with its index beside it: https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram - the pairs whose insert exceeds 1 kb on chromosome 20, cut out of that file by the command below and rehosted with its `.tbi`: https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz - the structural-variant callset over all 3,202 samples: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz - the reference the CRAM decodes against, and the RefSeq gene lane, as one hosted configuration with both tracks below already in it: https://jbrowse.org/demos/read_marks/config.json ## What a read carries Every record in a BAM or CRAM is a read with its position, and beside the position a handful of numbers the aligner wrote: its mapping quality, and for a paired read the distance from its start to its mate's end, the template length, signed by which mate is leftmost. An alignments track draws each read as a bar and can colour it by those numbers, but it cannot plot one of them on an axis. The mark display does that, over the same file, by naming the field. The window is 30 kb of an intron of _EFCAB8_ on chromosome 20, where the callset says NA12878 carries one copy of a 3.9 kb deletion. ## Depth as a coverage step The track is the CRAM as an `AlignmentsTrack`, opened with a `LinearMarkDisplay` whose one mark is a `bar`. Its `transform` runs a `coverage` step, which replaces the reads with runs of constant depth, and its `encoding` plots that depth. Nothing names a field of the read yet; the step wrote the one being plotted. ```json addtrack { "type": "AlignmentsTrack", "trackId": "na12878_read_depth", "name": "NA12878 reads (1000 Genomes, 30x)", "assemblyNames": ["hg38"], "adapter": { "type": "CramAdapter", "cramLocation": { "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram" }, "craiLocation": { "uri": "https://s3.amazonaws.com/1000genomes/1000G_2504_high_coverage/data/ERR3239334/NA12878.final.cram.crai" } }, "displays": [ { "type": "LinearMarkDisplay", "displayId": "na12878_read_depth-LinearMarkDisplay", "marks": [ { "shape": "bar", "transform": [{ "type": "coverage" }], "encoding": { "y": "coverage", "color": "#c8d8ee" } } ] } ] } ``` Open it at `chr20:32,925,000-32,955,000`.
A CRAM needs the reference to decode, and the track takes it from the assembly it is added to, so the config carries none. ## Insert size as a point per pair The second mark reads a field. Each pair appears twice in the file, once per mate, with the template length positive on the leftmost mate and negative on the other, so a `filter` step keeps the positive ones and the pair is counted once. A `point` then plots `template_length` on the y-axis, coloured by mapping quality through a ramp pinned to the 0 to 60 the aligner writes, so a grey point is a read the aligner was unsure of. The two numbers are on different scales, depth in the tens and insert size in the thousands, so the coverage mark moves to an axis of its own with `"resolve": "independent"` and takes the right-hand side of the plot. Both marks draw from one fetch of the reads. ```json "marks": [ { "shape": "bar", "transform": [{ "type": "coverage" }], "encoding": { "y": { "field": "coverage", "resolve": "independent" }, "color": "#c8d8ee" } }, { "shape": "point", "transform": [ { "type": "filter", "expr": "jexl:feature.template_length > 0 && feature.template_length < 8000" } ], "encoding": { "y": "template_length", "color": { "field": "score", "scale": "linear", "domain": [0, 60], "ramp": ["#bdbdbd", "#1f4e9a"] } } } ] ```
The upper group is the deletion measured a second way: each of those pairs was made at the library's insert and straddles the missing 3.9 kb, so it maps that much further apart, and it sits at the left edge of the dip because a pair is drawn at its leftmost read. The `score` a read answers is its mapping quality; the field is named that way on every track type the display attaches to, and a quality the aligner left unset is skipped rather than plotted at zero. Hovering a point reads out its position, its template length and its mapping quality; clicking one opens the read. ## Which reads carry the long inserts The reads themselves stack when the `span` shape reads a row a `stack` step wrote, which is the pileup an alignments track draws, said as two steps. A `formula` step first writes the unsigned insert into a field of its own, so both mates of a pair take the same colour, and a ramp over that field pinned at 5 kb paints a spanning pair red and every other pair in the pale blue of the depth bars. ```json "marks": [ { "shape": "span", "transform": [ { "type": "formula", "expr": "jexl:abs(feature.template_length)", "as": "insert" }, { "type": "stack" } ], "encoding": { "row": "row", "color": { "field": "insert", "scale": "linear", "domain": [0, 5000], "ramp": ["#c8d8ee", "#d62728"] } } } ] ``` Zoom to the left edge of the dip, `chr20:32,936,200-32,939,200`.
Pinning the `domain` keeps the key stable: an unpinned ramp spans the values on screen, and a window with no spanning pair would paint its longest ordinary insert red. ## Scanning the chromosome for the same signature A window of 30 kb is a few thousand reads; chromosome 20 is millions, and a track that has to fetch them all stops at its byte budget. The signature is carried by the few pairs whose insert is long, so those are cut out of the file once, one row per pair, into a BED whose header line names its columns. The `required_fields` option asks the CRAM decoder for the position, flags, quality and template length and no sequence, which is most of the decode: ```bash # one row per pair with an insert over 1 kb, from the leftmost mate to the # end of the insert, with the mapping quality in the score column # -q 20 drops reads the aligner could not place; -F 0x904 drops unmapped, # secondary and supplementary records # REF_PATH lets htslib fetch each reference sequence the CRAM names by MD5 export REF_PATH='https://www.ebi.ac.uk/ena/cram/md5/%s' samtools view -q 20 -F 0x904 --input-fmt-option required_fields=0x1DF NA12878.final.cram chr20 | awk 'BEGIN { OFS = "\t"; print "#chrom", "chromStart", "chromEnd", "name", "score", "strand", "tlen" } $7 == "=" && $9 > 1000 { print $3, $4 - 1, $4 - 1 + $9, $1, $5, "+", $9 }' | bgzip > NA12878.chr20.discordant_pairs.bed.gz tabix -p bed NA12878.chr20.discordant_pairs.bed.gz ``` Chromosome 20 gives 11,327 rows, small enough to fetch whole at any zoom. The track over them is a `FeatureTrack` with two marks: - **a `point` per pair**, its `tlen` on the y-axis, coloured by the mapping quality the header put in `score`. A `filter` keeps the inserts under 20 kb, since the centromere's pairs run to tens of megabases and would set the axis; `x2` is set to `start` so a pair is a point at its leftmost read rather than a bar across its whole span. - **a `bar` per bin** counting the deletion-sized pairs, 2 to 10 kb, the bin width following the zoom, on a right-hand axis whose `domain` is pinned at 60. The centromere holds thousands of such pairs and would otherwise be the axis; pinned, its bins saturate and a deletion's ten to fifty stand up on their own. ```json addtrack { "type": "FeatureTrack", "trackId": "na12878_chr20_pairs", "name": "NA12878 chr20, pairs over 1 kb", "assemblyNames": ["hg38"], "adapter": { "type": "BedTabixAdapter", "bedGzLocation": { "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz" }, "index": { "location": { "uri": "https://jbrowse.org/demos/read_marks/NA12878.chr20.discordant_pairs.bed.gz.tbi" } } }, "displays": [ { "type": "LinearMarkDisplay", "displayId": "na12878_chr20_pairs-LinearMarkDisplay", "marks": [ { "shape": "point", "transform": [ { "type": "filter", "expr": "jexl:feature.tlen < 20000" } ], "encoding": { "x2": "start", "y": "tlen", "color": { "field": "score", "scale": "linear", "domain": [0, 60], "ramp": ["#bdbdbd", "#1f4e9a"] } } }, { "shape": "bar", "transform": [ { "type": "filter", "expr": "jexl:feature.tlen > 2000 && feature.tlen < 10000" }, { "type": "bin", "step": "auto" }, { "type": "aggregate", "groupby": ["start", "end"], "ops": [{ "op": "count" }] } ], "encoding": { "y": { "field": "count", "resolve": "independent", "domain": [0, 60] }, "color": "#d62728" } } ] } ] } ```
The tallest bar outside the centromere, at 34.2 Mb, is a deletion NA12878 carries on both chromosomes; the one at 32.9 Mb is the intron above, on one. ## Checking the bars against the callset The same release's callset lists what its callers found in NA12878 on the chromosome. Every deletion over 2 kb on a non-reference genotype: ```bash # -s keeps one sample's genotypes; -i then keeps the rows where that sample # carries the allele bcftools view -s NA12878 1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz chr20 | bcftools query -i 'GT="alt" && INFO/SVTYPE="DEL" && INFO/SVLEN<-2000' \ -f '%CHROM\t%POS\t%END\t%INFO/SVLEN\t[%GT]\t%INFO/AF\t%INFO/EVIDENCE\n' ``` Against it, the count the bar mark makes: the pairs between 2 and 10 kb in the 100 kb window around each call. | position, chr20 | size | genotype | pairs 2 to 10 kb in the window | | --------------- | ------: | -------- | -----------------------------: | | 1.58 Mb | 33.1 kb | 0/1 | not in range | | 32.94 Mb | 3.9 kb | 0/1 | 18 | | 34.23 Mb | 3.3 kb | 1/1 | 46 | | 43.64 Mb | 2.7 kb | 0/1 | 10 | | 43.85 Mb | 2.6 kb | 0/1 | 25 | | 52.14 Mb | 2.1 kb | 1/1 | 32 | | 54.03 Mb | 10.9 kb | 0/1 | not in range | | 55.86 Mb | 6.0 kb | 0/1 | 16 | Every deletion the callset holds between 2 and 10 kb is a bar, and the two homozygous ones are the tallest, both chromosomes contributing pairs. Five other windows outside the centromere hold ten or more such pairs, at the start of the chromosome and at 1.4, 2.8, 32.7 and 48.5 Mb, and the callset has no deletion of this size at any of them. The first figure's claim was a dip and a group of long inserts over one window, so read the same window out of the file: the depth inside the call against the flank, and the pairs whose insert exceeds 2 kb. ```bash samtools coverage -r chr20:32937680-32941583 NA12878.final.cram | cut -f 1-3,7 samtools coverage -r chr20:32930000-32937000 NA12878.final.cram | cut -f 1-3,7 samtools view -q 20 NA12878.final.cram chr20:32935000-32944000 | awk '{ t = $9 < 0 ? -$9 : $9; if (t > 2000) big++; else if (t > 0) norm++ } END { print norm " pairs at the library insert, " big " over 2 kb" }' ``` | window | mean depth | | --------------------------- | ---------: | | chr20:32,937,680-32,941,583 | 15.6x | | chr20:32,930,000-32,937,000 | 34.1x | Across the 9 kb around the call, 1,688 reads sit at the library's insert and 41 carry one over 2 kb. ## Reproduce it end to end Every step above is wrapped in one script, [`build_read_marks.sh`](https://github.com/GMOD/jbrowse-components/blob/main/scripts/build_read_marks.sh): ```bash curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_read_marks.sh bash build_read_marks.sh # builds ./read_marks_build/jbrowse2 npx --yes serve read_marks_build/jbrowse2 # then open the printed URL ``` With no arguments it reads NA12878's CRAM over HTTP, cuts chromosome 20's long pairs out of it, writes a JBrowse with the two tracks above and prints the callset's deletions for the chromosome. Given your own BAM or CRAM and the FASTA it was aligned to, `bash build_read_marks.sh reads.cram genome.fa`, it builds the same two tracks over your file, and `CHROM` in the environment picks the chromosome to scan. The reads track needs only what every aligner writes: a mapping quality and, for paired reads, a template length. The tools it needs are the ones under [Prerequisites](#prerequisites). ## See also - [](https://jbrowse.org/jb2/docs/config_guides/mark_display) - [](https://jbrowse.org/jb2/docs/tutorials/alu_age) - [](https://jbrowse.org/jb2/docs/tutorials/mappability_qc) - [](https://jbrowse.org/jb2/docs/tutorials/sv_multisamples) - [](https://jbrowse.org/jb2/docs/quickstart_web) ## References - Byrska-Bishop M, et al. [High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios](https://doi.org/10.1016/j.cell.2022.08.004). _Cell_ 185:3426-3440 (2022), the reads and the structural-variant callset. - Li H, et al. [The Sequence Alignment/Map format and SAMtools](https://doi.org/10.1093/bioinformatics/btp352). _Bioinformatics_ 25:2078-2079 (2009), where the template length and mapping quality fields are defined.