Skip to main content

Variants on protein structures part 1: sickle cell from GWAS hit to haemoglobin S

JBrowse v5 beta

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.20. Desktop beta builds for Linux and macOS are on the releases page; the Windows build waits on code signing.

Sickle cell disease comes from one base change in HBB, the gene for the β chain of haemoglobin. Codon 7 changes from GAG to GTG, so glutamate 6 of the mature protein becomes valine. The variant, rs334, turns up again and again in genome-wide association studies, and population callsets show who has one copy and who has two. We put all of that in one JBrowse session beside the residue it changes, in a crystal structure of sickle haemoglobin. We:

  • open HBB on codon 7 from the protein browser with structure 2HBS, pin the human row of the alignment, and highlight and label the codon
  • add the GWAS Catalog track, then a copy whose filter and labels keep rs334's strongest record for each trait
  • add the 1000 Genomes VCF as a variant track, one row per genome, and sort the rows by genotype at rs334
  • color the genome, the alignment and the structure by the alignment's protein domains

Prerequisites

  • for your own JBrowse, the protein3d and msaview plugins from the plugin store; genomes.jbrowse.org already loads them

The data

Nothing to download: the protein browser, the plugins and the track configs below read these files by URL. Your own gene covers a gene of your own.

The files
What it isFile
Sickle haemoglobin, deoxy form, PDB 2HBS2HBS.cif
The NHGRI-EBI GWAS Catalog, which UCSC distributes as its GWAS Catalog track and the hosted hg38 config loadsgwas/
The 1000 Genomes 30x callset for chr11, 3,202 genomes20220422_3202_phased_SNV_INDEL_SV/
NCBI's ortholog report, the rows of the alignmentrest-api/

Opening HBB on the sickle cell codon

The protein browser opens a gene, a structure and an alignment as one linked session (Proteins on genomes.jbrowse.org walks through its page). Open it on HBB residue 7 with structure 2HBS:

  • open genomes.jbrowse.org/protein-browser/?gene=HBB&ref=9606&residue=7&structure=2hbs and click Open in JBrowse ↗
  • close the ClinVar and AlphaMissense tracks with the × on each track's label
  • type chr11:5226956-5227046[rev] into the location box. [rev] keeps the minus-strand gene reading left to right, the way it opened
  • click Reset Zoom, the circular arrow at the top right of the structure, to fit the haemoglobin tetramer
  • pick Use as reference row on the human_query tip of the alignment's tree: the human sequence pins above the alignment, and a residue matching it draws as a dot
  • drag across codon 7, the E of the gene's translation, in the genome view's ruler, and pick Highlight region
  • click the link chip at the top of the new band, pick Edit label... and type E6V. The genome view writes the label on the band, the alignment on column 7 and the structure beside Val6
  • hover column 7 of the human_query row in the alignment

Residue 7 counts the initiator methionine, which the mature protein loses, so the literature's Glu6 is the transcript's residue 7. The protein view's header reads both numbers off the hover and compares the two letters, glutamate in the transcript and valine in the crystal.

HBB at codon 7 with its ortholog alignment and 2HBS. The codon band, the alignment column and Val6 mark the same position. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The GWAS Catalog at rs334

The hosted hg38 config already has the GWAS Catalog, one record per study and trait:

  • choose Open track selector in the genome view's menu, search GWAS Catalog, and tick it. A variant studied many times draws as a stack
  • click a box to open its record. The records at rs334 include malaria, haemolytic anaemias, glycated haemoglobin (HbA1c) and kidney disease

A second copy of the catalog keeps rs334's strongest record for each of those traits, the one with the smallest p-value, and labels it with its trait, p-value, first author and year, read from the file's columns:

Goes in the tracks array of config.json. See Tracks.

{
  "type": "FeatureTrack",
  "trackId": "gwasCatalog_sickle_cell_traits",
  "name": "GWAS Catalog at rs334, one record per trait",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz"
    },
    "index": {
      "indexType": "CSI",
      "location": {
        "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz.csi"
      }
    }
  },
  "displayDefaults": {
    "filter": [
      "jexl:feature.name == 'rs334'",
      "jexl:feature.trait in ['Malaria', 'Hereditary hemolytic anemias (PheCode 282)', 'Glycated hemoglobin levels', 'Chronic Kidney Disease, Stage III (PheCode 585.33)']",
      "jexl:feature.pValue in ['4E-69', '1E-323', '3E-39', '9E-31']"
    ],
    "labels": {
      "name": "trait",
      "description": "jexl:'p=' + feature.pValue + ', ' + substring(feature.pubDate, 0, 4) + ', ' + feature.author"
    }
  }
}
jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "gwasCatalog_sickle_cell_traits",
  "name": "GWAS Catalog at rs334, one record per trait",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz"
    },
    "index": {
      "indexType": "CSI",
      "location": {
        "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz.csi"
      }
    }
  },
  "displayDefaults": {
    "filter": [
      "jexl:feature.name == '\''rs334'\''",
      "jexl:feature.trait in ['\''Malaria'\'', '\''Hereditary hemolytic anemias (PheCode 282)'\'', '\''Glycated hemoglobin levels'\'', '\''Chronic Kidney Disease, Stage III (PheCode 585.33)'\'']",
      "jexl:feature.pValue in ['\''4E-69'\'', '\''1E-323'\'', '\''3E-39'\'', '\''9E-31'\'']"
    ],
    "labels": {
      "name": "trait",
      "description": "jexl:'\''p='\'' + feature.pValue + '\'', '\'' + substring(feature.pubDate, 0, 4) + '\'', '\'' + feature.author"
    }
  }
}'

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": "gwasCatalog_sickle_cell_traits",
  "name": "GWAS Catalog at rs334, one record per trait",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "BedTabixAdapter",
    "bedGzLocation": {
      "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz"
    },
    "index": {
      "indexType": "CSI",
      "location": {
        "uri": "https://jbrowse.org/ucsc/hg38/gwasCatalog.bed.gz.csi"
      }
    }
  },
  "displayDefaults": {
    "filter": [
      "jexl:feature.name == 'rs334'",
      "jexl:feature.trait in ['Malaria', 'Hereditary hemolytic anemias (PheCode 282)', 'Glycated hemoglobin levels', 'Chronic Kidney Disease, Stage III (PheCode 585.33)']",
      "jexl:feature.pValue in ['4E-69', '1E-323', '3E-39', '9E-31']"
    ],
    "labels": {
      "name": "trait",
      "description": "jexl:'p=' + feature.pValue + ', ' + substring(feature.pubDate, 0, 4) + ', ' + feature.author"
    }
  }
}

Open this track in JBrowse Web ↗

To keep other traits, edit the filter with Filter by → Edit filters... in the track menu.

The GWAS Catalog under HBB, and the copy filtered to rs334's strongest record for each trait, labeled with its trait. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

Sickle allele in 1000 Genomes

The 1000 Genomes callset is a VCF with a column per genome, which JBrowse draws as one row per genome:

Goes in the tracks array of config.json. See Tracks.

{
  "type": "VariantTrack",
  "trackId": "kgp3202_chr11_snv",
  "name": "1000 Genomes 30x, chr11 SNVs and indels",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearVariantDisplay"
    }
  ]
}
jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "kgp3202_chr11_snv",
  "name": "1000 Genomes 30x, chr11 SNVs and indels",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearVariantDisplay"
    }
  ]
}'

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": "kgp3202_chr11_snv",
  "name": "1000 Genomes 30x, chr11 SNVs and indels",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr11.filtered.SNV_INDEL_SV_phased_panel.vcf.gz.tbi"
      }
    }
  },
  "displays": [
    {
      "type": "LinearVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearVariantDisplay"
    }
  ]
}

Open this track in JBrowse Web ↗

rs334 sits under the highlighted codon. Right-click its column and pick Sort rows by genotype here, which gathers the genomes with the allele into one block, nearly all of them with one copy. Two neighbours help read it:

  • rs713040, a synonymous change in the His3 codon four codons upstream, sits to its left on screen as the control, a common variant with one copy in some genomes and two in others
  • HbC, rs33930165, a second allele of the same codon, shares the highlighted column
1000 Genomes genotypes under HBB, one row per genome, sorted by genotype at rs334. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The domains across all three views

The alignment view fills its rows with NCBI's conserved domains and sites. The genome view and the structure read them from the alignment, so each domain keeps one color in every view:

  • in the alignment view's menu, click Show domains on the genome, which adds a domain track to the genome view
  • in the genome view's menu, untick Show → CDS reading frame colors, whose frame colors would mix with the domain colors
  • type chr11:5226956-5227046[rev] chr11:5226750-5226799[rev] into the location box, which adds the start of exon 2 beside codon 7
  • in the protein view's palette menu, pick Domains from the alignment
  • in the protein view's menu, pick Assembly → Deposited model (every chain), which draws both tetramers the crystal holds, then click Reset Zoom

In the protein view's header, the From the alignment rows list the same domains and sites.

HBB in three linked views: the genome's GWAS, 1000 Genomes and domain tracks, the alignment's rows and the 2HBS cartoon, colored by the alignment's domains. On the structure a cartoon color marks a domain and side-chain sticks mark a site. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The callouts follow rs334 from a statistic to a mechanism, one view at a time:

  • the GWAS records tie the allele to malaria and to haemolytic anaemia and kidney disease: one copy protects against severe malaria, and two copies cause sickle cell disease
  • in the 1000 Genomes track, nearly every carrier has one copy, the protected genotype
  • the alignment and the domain track mark no heme-binding or tetramer-interface site at Glu6, which sits on the protein's surface1
  • on the structure, the hydrophobic valine 6 of the second tetramer's chain H fills a pocket of Phe85 and Leu88 on the first tetramer's chain B. Deoxygenated haemoglobin S stacks into fibres through this contact, and the fibres bend red cells into sickles

Your own gene

Each step on this page works on any gene the protein browser opens:

  • the chips under the search box open other organisms on a residue each, from M. tuberculosis rpoB to dog SOD1, on genomes.jbrowse.org's GenArk assemblies. Proteins on genomes.jbrowse.org reads a grapevine gene on an AlphaFold model
  • Launch protein view in a gene's right-click menu picks the structure yourself, and its File or URL tab takes a PDB or mmCIF of your own, such as a ColabFold or AlphaFold model
  • any multi-sample VCF, bgzipped and tabix-indexed, adds the genotype layer as a variant track configured like the 1000 Genomes track above, and Sort rows by genotype here gathers the genomes with whichever variant you right-click. Part 2 slices one gene out of a large callset and colors the rows from a sample sheet

See also

Citations

Notes

  1. The Hb-beta-like domain, NCBI's cd08925, starts at the transcript's residue 8, one past the sickle residue. In 2HBS, helix A runs from Thr4, so Glu6 sits inside the globin fold; the bar's start marks where cd08925's model begins. ↩