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.16. 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, then a copy filtered to rs334's strongest record for each trait tied to sickle cell trait, and read what each record reports
  • add 1000 Genomes genotypes grouped by ancestry, with a common synonymous variant four codons upstream as the control
  • 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/
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 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

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 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 own 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"
    },
    "showLabels": "nameAndDescription"
  }
}
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"
    },
    "showLabels": "nameAndDescription"
  }
}'

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"
    },
    "showLabels": "nameAndDescription"
  }
}

Open this track in JBrowse Web ↗

The subset is a filter, which Filter by... → Edit filters... in the track menu edits to keep other traits.

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 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)

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 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 numbered callouts follow rs334 from a statistic to a mechanism, one view at a time:

  1. 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
  2. the 1000 Genomes carriers sit in the African-ancestry block of rows, nearly all with one copy, since malaria favoured carriers there
  3. the alignment and the domain track mark no heme-binding or tetramer-interface site at residue 6, so the change breaks none of the sites NCBI records1
  4. on the structure, valine 6 trades glutamate's charged side chain for a hydrophobic one on the β chain's surface
  5. in the 2HBS crystal, valine 6 of the second tetramer fills a pocket of Phe85 and Leu88 on this tetramer's chain B. Deoxygenated haemoglobin S stacks into fibres through that contact, and the fibres bend red cells into sickles

The protein view draws the crystal's first tetramer, and both the valine and the pocket it docks into sit on it.

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

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