A loss-of-function allele across breeds (Dog10K)
The tutorials target the JBrowse v5 beta, and the v4.3.0 release
on the download page lacks some of what they show. To install the
beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then
jbrowse create jbrowse2 --branch v5.0.0-beta.13. Desktop beta builds are coming soon.
Some dogs are poor metabolizers of the drugs the liver enzyme CYP1A2 clears, because they have a premature stop codon in CYP1A2. The literature names the variant only by its protein change, so we derive its genome coordinate by translating the reference coding sequence, slice the gene out of the 397 GB Dog10K SNV callset over HTTP, and read the genotypes across breeds with the wild canids as the control. A copy-number track then shows the gene's expansion across every canid in the collection.
Prerequisites
- the
UU_Cfam_GSD_1.0dog assembly (UCSC's canFam4) set up in JBrowse. The canFam4 hub on genomes.jbrowse.org is a config that loads it with its gene and repeat tracks, and the assemblies guide builds one by hand. bcftoolsbuilt with libcurlcurlpython3- htslib (
tabix) samtoolsbuilt with libcurl, for the build script's CRAM cross-check on the copy-number track, which is not a step on this page
On Debian/Ubuntu, apt install bcftools samtools tabix curl python3 covers it;
the packaged builds are linked against libcurl, so both can read the remote
callset and CRAMs. The scripts write local files, which
JBrowse Desktop opens by path and JBrowse Web takes
through Add track.
Where the data comes from
The Dog10K consortium's public share (Meadows et al. 2023), read directly over HTTP with no local copy of the 397 GB callset.
The build scripts take these files from their URLs, so there is nothing to download by hand.
- the SNV/indel callset the gene is sliced from, 397 GB over 1,987 canids: kiddlabshare.med.umich.edu/…/AutoAndXPAR.SNPs.vqsr99.vcf.gzhttps://kiddlabshare.med.umich.edu/dog10K/SNP_and_indel_calls_2021-10-17/AutoAndXPAR.SNPs.vqsr99.vcf.gz
- the sample table, breed and category per animal: kiddlabshare.med.umich.edu/…/dog10K-alignment-sample-table.2022-02-23-v7.txthttps://kiddlabshare.med.umich.edu/dog10K/sample-information/dog10K-alignment-sample-table.2022-02-23-v7.txt
- the reference sequence the stop codon is derived from, over UCSC's canFam4 REST API: api.genome.ucsc.edu/…/sequencehttps://api.genome.ucsc.edu/getData/sequence?genome=canFam4;chrom=chr30;start=38258000;end=38265000
- the RefSeq gene structure that same derivation reads exon boundaries from: api.genome.ucsc.edu/…/trackhttps://api.genome.ucsc.edu/getData/track?genome=canFam4;track=ncbiRefSeqCurated;chrom=chr30;start=38258000;end=38265000
- the 15 published CRAMs the copy-number track validates callset depth against: kiddlabshare.med.umich.edu/…/cram-sharehttps://kiddlabshare.med.umich.edu/dog10K/cram-share/
Loading the UU_Cfam_GSD_1.0 dog assembly
The tracks name UU_Cfam_GSD_1.0, the Dog10K reference that UCSC calls canFam4.
We load it from UCSC's 2bit, with the alias file that maps the chr names to
GenBank accessions:
Goes in the assemblies array of config.json. See Assemblies.
{
"name": "UU_Cfam_GSD_1.0",
"aliases": ["canFam4"],
"sequence": {
"adapter": {
"type": "TwoBitAdapter",
"uri": "https://hgdownload.soe.ucsc.edu/goldenPath/canFam4/bigZips/canFam4.2bit"
}
},
"refNameAliases": {
"uri": "https://hgdownload.soe.ucsc.edu/goldenPath/canFam4/bigZips/canFam4.chromAlias.txt"
}
}jbrowse add-assembly https://hgdownload.soe.ucsc.edu/goldenPath/canFam4/bigZips/canFam4.2bit \
--name UU_Cfam_GSD_1.0 \
--alias canFam4 \
--refNameAliases https://hgdownload.soe.ucsc.edu/goldenPath/canFam4/bigZips/canFam4.chromAlias.txtThe CYP1A2 nonsense variant
CYP1A2 is a drug-metabolizing cytochrome P450 in which dogs have a nonsense variant. The truncating allele acts recessively: liver microsomes from dogs homozygous for it have no CYP1A2 protein and those dogs are poor metabolizers of drugs the enzyme clears, while heterozygotes express the enzyme normally (Mise et al. 2004). We ask which breeds have the allele and whether wild canids do. Wolves are the control, because an allele they share predates domestication.
Deriving the CYP1A2 stop codon's coordinate from the reference
The literature names this variant by its protein consequence, p.Arg373Ter. The build script rebuilds CYP1A2's coding sequence from the reference and the RefSeq exon structure, translates it, and reports codon 373:
NM_001008720.1 CDS 1539 bp, 513 aa
codon 373 = CGA (R) at chr30:38261635
C>T at that first base makes TGA, a stopThe C>T at chr30:38,261,635 truncates the protein at residue 373 of 513. Checking the callset at exactly that position finds it, passing every filter:
SNVS=https://kiddlabshare.med.umich.edu/dog10K/SNP_and_indel_calls_2021-10-17/AutoAndXPAR.SNPs.vqsr99.vcf.gz
bcftools query -r chr30:38261635-38261636 -f '%POS\t%REF\t%ALT\t%FILTER\t%AC\t%AN\n' \
"$SNVS"
# chr30 38261635 C T PASS 174 3974Slicing CYP1A2 out of the Dog10K SNV callset
The Dog10K SNV callset is a single 397 GB VCF over 1,987 canids with a tabix
index, so bcftools fetches one gene over HTTP:
SNVS=https://kiddlabshare.med.umich.edu/dog10K/SNP_and_indel_calls_2021-10-17/AutoAndXPAR.SNPs.vqsr99.vcf.gz
# --force-samples: proceed even if a name in cyp.samples isn't in the VCF's
# sample list, instead of exiting
bcftools view -r chr30:38258000-38265000 -S cyp.samples --force-samples \
-Oz -o dog10k_cyp1a2_snvs.vcf.gz "$SNVS"
tabix -p vcf dog10k_cyp1a2_snvs.vcf.gzcyp.samples holds breeds that have the allele, two that do not, and four Greek
gray wolves.
Loading the CYP1A2 slice as a variant track
An SNV VCF loads as an ordinary VariantTrack:
Goes in the tracks array of config.json. See Tracks.
{
"trackId": "dog10k_cyp1a2_snvs",
"name": "Dog10K SNVs at CYP1A2",
"uri": "dog10k_cyp1a2_snvs.vcf.gz",
"assemblyNames": ["UU_Cfam_GSD_1.0"]
}jbrowse add-track dog10k_cyp1a2_snvs.vcf.gz \
--trackId dog10k_cyp1a2_snvs \
--name "Dog10K SNVs at CYP1A2" \
--assemblyNames UU_Cfam_GSD_1.0 \
--load copyIn JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track... and, in Add a track from file or URL, enter:
- Main file:
dog10k_cyp1a2_snvs.vcf.gz
Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:
- Track name:
Dog10K SNVs at CYP1A2 - Assembly:
UU_Cfam_GSD_1.0
Click Add.
dog10k_cyp1a2_snvs.vcf.gz is relative to a config.json. Replace it with its URL or its path on this computer.
The display draws one row per sample, labelled with the Dog10K IDs. In
grammar-of-graphics terms, the figure's rows orders the rows through domain
and names them through labels ({ "GHND000001": "German Hound 1" }), and
rowColor, a color scale on each row's label bar, tints them.
Edit colors/arrangement... in the track menu writes both as you rename and
recolor rows, and Advanced → Edit plot... shows them as text. For a larger
panel, a samplesTsvLocation on the adapter supplies the names instead.
A whole-gene view of the slice is a field of one-pixel ticks, so zoom to the codon. At base level each sample's call is a block, and the gene track shows which exon it sits in.
CYP1A2 stop-gained genotypes across breeds and wolves
The build script genotypes the stop site over every canid in the callset. Dozens of breeds have the allele and it reaches homozygosity in several: every German Hound and every Shetland Sheepdog sampled here has at least one copy, while every wolf and every coyote in the collection is homozygous reference.
Three other variants lie within about a hundred bases of the stop site, and the config below filters them out. Dropping the filter shows them: two are reference in every animal of this sample set, including the one at the codon's second base, so each draws an empty column, and the third sits 15 bp along, and every wolf has it.
Goes in the tracks array of config.json. See Tracks.
{
"type": "VariantTrack",
"trackId": "dog10k_cyp1a2_snvs",
"name": "Dog10K SNVs at CYP1A2",
"uri": "dog10k_cyp1a2_snvs.vcf.gz",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"filter": ["jexl:feature.start == 38261634"]
}
]
}jbrowse add-track-json '{
"type": "VariantTrack",
"trackId": "dog10k_cyp1a2_snvs",
"name": "Dog10K SNVs at CYP1A2",
"uri": "dog10k_cyp1a2_snvs.vcf.gz",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"filter": ["jexl:feature.start == 38261634"]
}
]
}'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": "VariantTrack",
"trackId": "dog10k_cyp1a2_snvs",
"name": "Dog10K SNVs at CYP1A2",
"uri": "dog10k_cyp1a2_snvs.vcf.gz",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"filter": ["jexl:feature.start == 38261634"]
}
]
}dog10k_cyp1a2_snvs.vcf.gz is relative to a config.json. Replace it with its URL or its path on this computer.
Copy number at CYP1A2
The Dog10K paper reports half the collection at three or more copies of
CYP1A2. The per-animal estimates are unpublished, but the SNV callset records
a per-sample read depth (DP) at every site, so one slice of it, stripped to
that field, covers every canid in the collection:
# -r reads only the locus over HTTP; -x keeps FORMAT/DP and drops the rest
bcftools view -r chr30:38205000-38400000 -Ou "$SNVS" |
bcftools annotate -x 'INFO,^FORMAT/DP' -Oz -o dp.vcf.gz
bcftools query -l dp.vcf.gz > cohort.samples
bcftools query -f '%POS[\t%DP]\n' dp.vcf.gz > cohort.dpThe build script converts depth to copy number within each dog, taking the sequence around the element in that dog as two copies:
CN = 2 * depth over the element / depth over the sequence around itThe flanking sequence should come back at two copies in every dog, which checks the ratio.
Each window is 5 kb of depth stepped by 1 kb, so a call rests on 5 kb of evidence and is painted at 1 kb resolution.
The callset records depth only where a variant was called, so the build script
checks it against the 15 CRAMs the Dog10K share publishes. Over the shared
windows the two depth sources agree closely, with no bias. The CRAM-based
painting is in the config as dog10k_cyp1a2_cn.
The output is a BED with the color in itemRgb and the rounded call in
copyNumber. The identity color scale uses the file's own colors, and its
labels pair each color with a copy number so the key reads as copy number. The
figure draws one key, from the named-animals track above this one (described
below), so this track sets showLegend to false. Drop that line to get a key
on this track:
Goes in the tracks array of config.json. See Tracks.
{
"type": "FeatureTrack",
"trackId": "dog10k_cyp1a2_cohort_cn",
"name": "CYP1A2 copy number, every canid",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "dog10k_cyp1a2_cohort_cn.bed.gz",
"disableGeneHeuristic": true,
"columnNames": [
"chrom",
"chromStart",
"chromEnd",
"name",
"score",
"strand",
"thickStart",
"thickEnd",
"itemRgb",
"sample",
"copyNumber"
]
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": "sample",
"showLegend": false,
"color": {
"scale": "identity",
"title": "Copy number",
"domain": [
"rgb(33,102,172)",
"rgb(146,197,222)",
"rgb(224,224,224)",
"rgb(244,165,130)",
"rgb(214,96,77)",
"rgb(178,24,43)",
"rgb(103,0,31)"
],
"labels": ["CN 0", "CN 1", "CN 2", "CN 3", "CN 4", "CN 5", "CN 6+"]
}
}
]
}jbrowse add-track-json '{
"type": "FeatureTrack",
"trackId": "dog10k_cyp1a2_cohort_cn",
"name": "CYP1A2 copy number, every canid",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "dog10k_cyp1a2_cohort_cn.bed.gz",
"disableGeneHeuristic": true,
"columnNames": [
"chrom",
"chromStart",
"chromEnd",
"name",
"score",
"strand",
"thickStart",
"thickEnd",
"itemRgb",
"sample",
"copyNumber"
]
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": "sample",
"showLegend": false,
"color": {
"scale": "identity",
"title": "Copy number",
"domain": [
"rgb(33,102,172)",
"rgb(146,197,222)",
"rgb(224,224,224)",
"rgb(244,165,130)",
"rgb(214,96,77)",
"rgb(178,24,43)",
"rgb(103,0,31)"
],
"labels": ["CN 0", "CN 1", "CN 2", "CN 3", "CN 4", "CN 5", "CN 6+"]
}
}
]
}'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": "dog10k_cyp1a2_cohort_cn",
"name": "CYP1A2 copy number, every canid",
"assemblyNames": ["UU_Cfam_GSD_1.0"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "dog10k_cyp1a2_cohort_cn.bed.gz",
"disableGeneHeuristic": true,
"columnNames": [
"chrom",
"chromStart",
"chromEnd",
"name",
"score",
"strand",
"thickStart",
"thickEnd",
"itemRgb",
"sample",
"copyNumber"
]
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": "sample",
"showLegend": false,
"color": {
"scale": "identity",
"title": "Copy number",
"domain": [
"rgb(33,102,172)",
"rgb(146,197,222)",
"rgb(224,224,224)",
"rgb(244,165,130)",
"rgb(214,96,77)",
"rgb(178,24,43)",
"rgb(103,0,31)"
],
"labels": ["CN 0", "CN 1", "CN 2", "CN 3", "CN 4", "CN 5", "CN 6+"]
}
}
]
}dog10k_cyp1a2_cohort_cn.bed.gz is relative to a config.json. Replace it with its URL or its path on this computer.
The build script prints this color block from the palette it painted with. The
figure's upper track is a second BED the script writes,
dog10k_cyp1a2_breed_cn.bed.gz, whose sample column already holds breed
names; its track takes rows: { "field": "sample", "domain": [...] } with those
names in drawing order.
The upper track holds every Golden Retriever, Labrador Retriever and Boxer in the collection, plus the four wolves from the genotype figure. Every Golden has the expansion, every Boxer has two copies, and the Labradors split one dog to the next.
The four wolves, the control, all have the expansion, so unlike the stop-gained allele it is shared with wild canids and predates domestication. Their calls rest on callset depth alone, since none of the dogs with published reads is a wolf.
The white stripes through both tracks are windows with no call. A window whose median across the whole collection is not two copies measures a quirk of the reference, so the build script drops it from every row.1
The lower track repeats the estimate for every canid.
Clustering → Cluster rows by similarity... in its track menu, or
runClustering in a session's display entry, orders the rows by each animal's
profile across the window. Clustering groups animals whose expansion starts and
ends in the same place, and the blocks on either side of the gene are deletion
polymorphisms.
The depth-based estimate puts far more of the collection at three or more copies than the paper reports. The two counts cover different intervals: the paper ran QuicK-mer2 over an element whose extent is unpublished, and the build script counts the windows the collection puts above two.
Reproduce it end to end
build_dog10k_cyp1a2.sh
builds the genotype track. It:
- derives the stop codon's position from the reference
- builds the sample list from the Dog10K sample table and slices the gene out of the callset
- prints the genotypes at the stop so you can check the figure against the data
- genotypes that one site over all 1,987 canids for the breed and wild-canid counts quoted above
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_dog10k_cyp1a2.sh
bash build_dog10k_cyp1a2.sh # writes ./dog10k_cyp1a2_build/build_dog10k_cyp1a2_cn.sh
builds the copy-number tracks. It:
- reads depth over the gene and the sequence around it from each of the 15 published CRAMs, counting only positions RepeatMasker leaves unmasked, the same restriction QuicK-mer2's unique k-mers make
- scales each window against that dog's own flanks, taken as two copies, and prints each dog's copy number over the element beside the spread of its flanks around two
- takes the same ratio from the callset's per-sample depth for every canid, and prints how closely it agrees with the CRAM estimate in the 15 dogs
- drops any flanking window whose median across the animals is off two, rounds the rest to whole copies, and writes one painted row per animal
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_dog10k_cyp1a2_cn.sh
bash build_dog10k_cyp1a2_cn.sh # writes ./dog10k_cyp1a2_cn_build/See also
- Tutorial: Structural variants (Dog10K)
- Tutorial: A selected haplotype (Dog10K)
- Tutorial: Local ancestry (Dog10K)
- User guide: Multi-sample variant display
- Config guide: Variant track
Citations
- Meadows et al. (2023). Genome sequencing of 2000 canids by the Dog10K consortium advances the understanding of demography, genome function and architecture
- Court (2013). Canine cytochrome P450 pharmacogenetics
Notes
-
The widest stripe sits on a CpG island, whose high GC content lowers read depth in every canid, and each 5 kb window spreads that over the blocks around it. ↩
Feedback on this tutorial is welcome: contact us.