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.15. Desktop beta builds are coming soon.

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 its frequency differs sharply between ancestries. We put all of that in one JBrowse session beside the residue it changes, in a crystal structure of sickle haemoglobin. We'll:

  • open HBB on codon 7 in the protein browser, with structure 2HBS, and hover the codon's alignment column to find its residue
  • add the GWAS Catalog and read what has been associated with rs334
  • add 1000 Genomes genotypes grouped by ancestry, with a common synonymous variant four codons upstream as the control

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 config 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/
The population and superpopulation of each of those genomes20130606_g1k_3202_samples_ped_population.txt
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 both haemoglobin tetramers in the asymmetric unit
  • 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

Open the track selector from the genome view's menu, Open track selector, search GWAS Catalog, and tick it. The catalog has one record per study and trait, so a variant studied many times draws as a stack. Click a box to open its record.

The GWAS Catalog under HBB. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The records at rs334 include malaria, haemolytic anaemias, glycated haemoglobin (HbA1c) and kidney function.

Sickle allele by ancestry in 1000 Genomes

The 1000 Genomes callset is a VCF with a column per genome, and JBrowse can draw one row per genome, colored by a column of a sample sheet. The release's ped file holds each genome's population and superpopulation. This turns it into the sheet JBrowse reads, sorted so each superpopulation is one block of rows:

awk 'BEGIN{OFS="\t"} NR>1 {print $2,$6,$7}' samples_ped_population.txt \
  | sort -t$'\t' -k3,3 -k2,2 \
  | (printf 'name\tpopulation\tsuperpopulation\n'; cat) > samples.tsv

Add the callset as a track, with a hosted copy of that sheet:

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"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/1000g/1kgp_3202_superpopulations.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearMultiSampleVariantDisplay",
      "rowColor": "superpopulation"
    }
  ]
}
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"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/1000g/1kgp_3202_superpopulations.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearMultiSampleVariantDisplay",
      "rowColor": "superpopulation"
    }
  ]
}'

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"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/1000g/1kgp_3202_superpopulations.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "kgp3202_chr11_snv-LinearMultiSampleVariantDisplay",
      "rowColor": "superpopulation"
    }
  ]
}

Open this track in JBrowse Web ↗

Read two columns. rs334 sits under the highlighted codon, and rs713040, a synonymous change in the His3 codon four codons upstream, sits to its left on screen. HbC, rs33930165, a second allele of the same codon, shares the highlighted column.

1000 Genomes genotypes under HBB, one row per genome, colored down the left edge by superpopulation. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

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. Part 2 takes dog SOD1 onto the assembly its population callset uses, and part 3 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 VCF with a sample sheet adds the population layer, colored by whichever column of the sheet you name in rowColor

See also

Citations