Variants on protein structures part 2: degenerative myelopathy in dogs

Degenerative myelopathy is a late-onset paralysis of dogs that resembles amyotrophic lateral sclerosis in people. A genome-wide association study found the risk allele in SOD1: c.118G>A changes glutamate 40 of superoxide dismutase to lysine (E40K), and most affected dogs are homozygous for it (Awano et al. 2009). We put the codon beside a crystal structure of the E40K protein itself, then read the allele across the Dog10K collection, with wolves and coyotes as the control, since an allele they share would predate domestication. We'll:

• open SOD1 on codon 40 of the canFam4 dog assembly and launch structure 7WX0 and an ortholog alignment from the gene's right-click menu
• slice SOD1 out of the Dog10K SNV callset and make a sample sheet of each animal's category
• add the slice as a track, one row per animal, colored by category

Part 1 opened its gene from the protein browser. The protein browser opens dog genes on a Labrador assembly, and Dog10K is called against canFam4, so here we launch the same three views ourselves on canFam4.

Prerequisites#

• for your own JBrowse, the protein3d and msaview plugins from the plugin store; genomes.jbrowse.org already loads them
• bcftools built with libcurl, and htslib (tabix)
• curl

# Debian/Ubuntu
sudo apt install bcftools tabix
# or conda, which has every tool above
conda install -c bioconda -c conda-forge bcftools htslib

The data#

The Dog10K consortium's public share (Meadows et al. 2023). Nothing to download: the commands below read these files by URL. The slicing step reads any tabix-indexed VCF over HTTP, so a callset of your own slices the same way.

Opening SOD1 on codon 40 with the E40K structure#

Codon 40 of NM_001003035.1 sits in SOD1's second coding exon, at chr31:27,123,057-27,123,059 on canFam4. Open the canFam4 hub on genomes.jbrowse.org and:

• type chr31:27123013-27123103 into the location box
• in the track selector, tick NCBI RefSeq - RefSeq Curated
• right-click the SOD1 transcript and choose Launch protein view. On the PDB search tab pick 7WX0, then choose Open side by side from the arrow beside Launch
• right-click the transcript again and choose Launch MSA view, which opens on its Orthologs tab; press Submit
• click Reset Zoom, the circular arrow at the top right of the structure, to fit the whole structure
• hover column 40 of the dog_query row in the alignment

The hover lights the codon on the gene, the column in the alignment and the residue on the structure, and the protein view's header compares the two letters: glutamate in the dog genome, lysine in the crystal.

Slicing SOD1 out of Dog10K#

The Dog10K callset has a tabix index, so bcftools reads one gene over HTTP. This keeps the genotypes and drops the rest, which leaves a file of a few tens of kilobytes:

SNVS=https://kiddlabshare.med.umich.edu/dog10K/SNP_and_indel_calls_2021-10-17/AutoAndXPAR.SNPs.vqsr99.vcf.gz
# -x INFO,^FORMAT/GT: drop every INFO field and every FORMAT field but GT
bcftools annotate -r chr31:27118000-27128000 -x INFO,^FORMAT/GT \
  -Oz -o dog10k_sod1.vcf.gz "$SNVS"
tabix -p vcf dog10k_sod1.vcf.gz

The sample table files each animal under a category:

curl -o sample-table.txt https://kiddlabshare.med.umich.edu/dog10K/sample-information/dog10K-alignment-sample-table.2022-02-23-v7.txt

An awk command turns the table into the sheet JBrowse reads, with wolves first, then coyotes, village dogs, mixed and other dogs, and breed dogs, sorted by breed within each:

awk -F'\t' 'BEGIN { OFS = "\t"
    n = split("Wolf Coyote Village_Dogs Mixed/Other Breed_Dogs", o, " ")
    for (i = 1; i <= n; i++) rank[o[i]] = i }
  NR > 1 { print rank[$3], $1, $3, $2 }' sample-table.txt \
  | sort -t$'\t' -k1,1n -k4,4 -k2,2 | cut -f2- \
  | (printf 'name\tcategory\tbreed\n'; cat) > samples.tsv

E40K across breeds, village dogs and wild canids#

Add the slice as a track, with hosted copies of the two files:

{
  "type": "VariantTrack",
  "trackId": "dog10k_sod1",
  "name": "Dog10K SNVs at SOD1",
  "assemblyNames": ["canFam4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz.tbi"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "dog10k_sod1-LinearMultiSampleVariantDisplay",
      "rowColor": "category"
    }
  ]
}

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "dog10k_sod1",
  "name": "Dog10K SNVs at SOD1",
  "assemblyNames": ["canFam4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz.tbi"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "dog10k_sod1-LinearMultiSampleVariantDisplay",
      "rowColor": "category"
    }
  ]
}'

{
  "type": "VariantTrack",
  "trackId": "dog10k_sod1",
  "name": "Dog10K SNVs at SOD1",
  "assemblyNames": ["canFam4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "vcfGzLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz"
    },
    "index": {
      "location": {
        "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_sod1.vcf.gz.tbi"
      }
    },
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/dog10k_sod1/dog10k_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "displayId": "dog10k_sod1-LinearMultiSampleVariantDisplay",
      "rowColor": "category"
    }
  ]
}

To put a site where the wolves have the alternate allele beside codon 40, type chr31:27122805-27122834 chr31:27123013-27123103 into the location box. The first region is 30 bp of the upstream intron, around an SNV at chr31:27,122,819.

Counting E40K carriers by category#

The command below joins each animal's genotype at chr31:27,123,057 to its category in samples.tsv and prints carriers, then animals, per category:

bcftools query -r chr31:27123057 -f '[%SAMPLE\t%GT\n]' dog10k_sod1.vcf.gz \
  | awk -F'\t' 'NR == FNR { cat[$1] = $2; next }
      { n[cat[$1]]++; if ($2 ~ /1/) carriers[cat[$1]]++ }
      END { for (c in n) print c, carriers[c] + 0, n[c] }' samples.tsv - \
  | sort

Breed_Dogs 165 1593
Coyote 0 4
Mixed/Other 3 18
Village_Dogs 30 309
Wolf 0 63

See also#

• Tutorial: Variants on protein structures part 1: sickle cell from GWAS hit to haemoglobin S
• Tutorial: Variants on protein structures part 3: seedless grapes on an AlphaFold model
• Tutorial: A loss-of-function allele across breeds (Dog10K)
• Tutorial: Proteins on genomes.jbrowse.org

External links#

• Dog10K data share
• jbrowse-plugin-protein3d
• jbrowse-plugin-msaview

Citations#

• Awano T, Johnson GS, Wade CM, et al. Genome-wide association analysis reveals a SOD1 mutation in canine degenerative myelopathy that resembles amyotrophic lateral sclerosis. Proc Natl Acad Sci U S A 2009. https://doi.org/10.1073/pnas.0812297106
• Hashimoto K, Watanabe S, Akutsu M, et al. Intrinsic structural vulnerability in the hydrophobic core induces species-specific aggregation of canine SOD1 with degenerative myelopathy-linked E40K mutation. J Biol Chem 2023. https://doi.org/10.1016/j.jbc.2023.104798
• Meadows JRS, Kidd JM, Wang GD, et al. Genome sequencing of 2000 canids by the Dog10K consortium advances the understanding of demography, genome function and architecture. Genome Biol 2023. https://doi.org/10.1186/s13059-023-03023-7

Feedback on this tutorial is welcome: contact us.