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.

Prerequisites

a JBrowse to open the figures' sessions in (Web or Desktop); every file here is a URL, so nothing needs hosting to read along
samtools and htslib (bgzip, tabix), for cutting the pairs out of the file and for checking a window by hand
bcftools, for reading the callset at the end
Node.js and the JBrowse 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), 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; what it cannot do is plot one of them on an axis. The mark display does exactly 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.

Config file CLI (add-track-json) JBrowse Desktop

{
  "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" }
        }
      ]
    }
  ]
}

jbrowse add-track-json '{
  "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" }
        }
      ]
    }
  ]
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "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.

Thirty kilobases of an EFCAB8 intron in NA12878, the read depth as bars. Between 32,937,500 and 32,941,500 the depth runs at about half of what it is on either side. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

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.

"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 same window with each pair's insert size as a point on the left axis and the depth on its own axis on the right. The pairs sit in a band under 1,000 bases, and over the left edge of the dip a second group appears between 4,300 and 4,700 bases, in the full blue of a mapping quality of 60. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

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.

"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.

The left breakpoint at 3 kb, the reads stacked and coloured by their pair's insert. The red reads end together at 32,937,680, where their mates lie 4 kb to the right; the pale reads run across it, and thin out on the far side. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The pinned domain is what makes 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:

# 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.

Config file CLI (add-track-json) JBrowse Desktop

{
  "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"
          }
        }
      ]
    }
  ]
}

jbrowse add-track-json '{
  "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"
          }
        }
      ]
    }
  ]
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "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"
          }
        }
      ]
    }
  ]
}

Chromosome 20 end to end. Every pair with an insert under 20 kb is a point at its insert size, and the red bars on the right axis count the pairs between 2 and 10 kb per bin. The centromere, from 26 to 32 Mb, saturates both; outside it the bars stand up in a handful of places, each under a short stack of dark points. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

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:

# -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.

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:

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. What the reads track needs from any file is what every aligner writes: a mapping quality and, for paired reads, a template length. The tools it needs are the ones under Prerequisites.

References

Byrska-Bishop M, et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell 185:3426-3440 (2022), the reads and the structural-variant callset.
Li H, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078-2079 (2009), where the template length and mapping quality fields are defined.

Feedback on this tutorial is welcome: contact us.

A deletion from a BAM's own fields (NA12878, insert size)

Prerequisites#

Where the data comes from#

What a read carries#

Depth as a coverage step#

Insert size as a point per pair#

Which reads carry the long inserts#

Scanning the chromosome for the same signature#

Checking the bars against the callset#

Reproduce it end to end#

See also#

References#