AlphaMissense predictions on pathogenicity of amino acid changes can be combined with AlphaFold (or other) predictions of protein structure.
Fast path
Figure 3F of the AlphaMissense publication visualizes mean pathogenicity for UniProt id P35557. Filter amino acid data for that identifier
P35557_aa <-
am_data("aa_substitutions") |>
dplyr::filter(uniprot_id == "P35557")
and visualization median pathogenicity with
af_prediction_view(P35557_aa)
The image is interactive, including rotation and zoom. The following sections explore this visualization in more detail.
UniProt identifiers
Both AlphaMissense and AlphaFold use UniProt identifiers. Find all AlphaMissense amino acid substitutions with UniProt identifiers starting with P3555; the choice of this identifier is so that results can be compared with Figure 3F of the AlphaMissense publication.
uniprot_ids <-
am_data("aa_substitutions") |>
dplyr::filter(uniprot_id %like% "P3555%") |>
dplyr::distinct(uniprot_id) |>
pull(uniprot_id)
uniprot_ids
#> [1] "P35558" "P35556" "P35557" "P35555"
The AlphaMissenseR package includes several functions that facilitate interaction with AlphaFold; these functions start with af_*(). Use af_predictions() to discover AlphaFold predictions (via the AlphaFold API) associated with UniProt identifiers.
prediction <- af_predictions(uniprot_ids)
#> * [17:35:02][info] 2 of 4 uniprot accessions not found
#> 'P35556', 'P35555'
glimpse(prediction)
#> Rows: 2
#> Columns: 25
#> $ entryId <chr> "AF-P35558-F1", "AF-P35557-F1"
#> $ gene <chr> "PCK1", "GCK"
#> $ sequenceChecksum <chr> "78D309E0845CC181", "094D4A2F78096724"
#> $ sequenceVersionDate <chr> "2006-03-07", "1994-06-01"
#> $ uniprotAccession <chr> "P35558", "P35557"
#> $ uniprotId <chr> "PCKGC_HUMAN", "HXK4_HUMAN"
#> $ uniprotDescription <chr> "Phosphoenolpyruvate carboxykinase, cytosolic […
#> $ taxId <int> 9606, 9606
#> $ organismScientificName <chr> "Homo sapiens", "Homo sapiens"
#> $ uniprotStart <int> 1, 1
#> $ uniprotEnd <int> 622, 465
#> $ uniprotSequence <chr> "MPPQLQNGLNLSAKVVQGSLDSLPQAVREFLENNAELCQPDHIHIC…
#> $ modelCreatedDate <chr> "2022-06-01", "2022-06-01"
#> $ latestVersion <int> 4, 4
#> $ allVersions <list> [1, 2, 3, 4], [1, 2, 3, 4]
#> $ isReviewed <lgl> TRUE, TRUE
#> $ isReferenceProteome <lgl> TRUE, TRUE
#> $ cifUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ bcifUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ pdbUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ paeImageUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ paeDocUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ amAnnotationsUrl <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ amAnnotationsHg19Url <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
#> $ amAnnotationsHg38Url <chr> "https://alphafold.ebi.ac.uk/files/AF-P35558-F1…
Note the message indicating that some UniProt identifiers (accessions) are not found in the AlphaFold database. The query returns a tibble containing columns with information on organism and UniProt characteristics (including protein sequence) , as well as URLs for files representing three-dimensional protein structure. We will use pdbUrl.
Protein structure
Focus on a particular UniProt identifier and the PDB url.
pdb_url <-
prediction |>
dplyr::filter(uniprotAccession == "P35557") |>
dplyr::pull(pdbUrl)
Cache the PDB file using BiocFileCache, and read the PDB file using bio3d.
pdb_file <- BiocFileCache::bfcrpath(rnames = basename(pdb_url), fpath = pdb_url)
pdb <- bio3d::read.pdb(pdb_file)
pdb
#>
#> Call: bio3d::read.pdb(file = pdb_file)
#>
#> Total Models#: 1
#> Total Atoms#: 3642, XYZs#: 10926 Chains#: 1 (values: A)
#>
#> Protein Atoms#: 3642 (residues/Calpha atoms#: 465)
#> Nucleic acid Atoms#: 0 (residues/phosphate atoms#: 0)
#>
#> Non-protein/nucleic Atoms#: 0 (residues: 0)
#> Non-protein/nucleic resid values: [ none ]
#>
#> Protein sequence:
#> MLDDRARMEAAKKEKVEQILAEFQLQEEDLKKVMRRMQKEMDRGLRLETHEEASVKMLPT
#> YVRSTPEGSEVGDFLSLDLGGTNFRVMLVKVGEGEEGQWSVKTKHQMYSIPEDAMTGTAE
#> MLFDYISECISDFLDKHQMKHKKLPLGFTFSFPVRHEDIDKGILLNWTKGFKASGAEGNN
#> VVGLLRDAIKRRGDFEMDVVAMVNDTVATMISCYYEDHQCEVGMI...<cut>...MLGQ
#>
#> + attr: atom, xyz, seqres, calpha, call
Visualize the protein using r3dmol, using the ‘cartoon’ style.
r3dmol::r3dmol() |>
## use the PDB representation
r3dmol::m_add_model(r3dmol::m_bio3d(pdb)) |>
## visualize as a 'cartoon' with alpha helices and beta sheets
r3dmol::m_set_style(style = r3dmol::m_style_cartoon(arrow = TRUE)) |>
## fit molecule into display area
r3dmol::m_zoom_to()
Average pathogenicity
Our goal is to visualize some measure of ‘average’ pathogenicity on the three-dimensional protein structure provided by AlphaFold. Start with a specific genome sequence (e.g., hg38). Filter to the amino acids in our UniProt region of interest.
P35557 <-
am_data("hg38") |>
dplyr::filter(uniprot_id == "P35557")
At each chromosome position, the AlphaMissense predictions contain several alternative alleles and hence protein variants. The (arithmetic) average pathogenicity (this is an extremely naive computation) at each amino acid position is
pathogenicity <- am_aa_pathogenicity(P35557)
pathogenicity
#> # A tibble: 464 × 9
#> uniprot_id aa_pos aa_ref aa_pathogenicity_n aa_pathogenicity_mean
#> <chr> <int> <chr> <int> <dbl>
#> 1 P35557 2 L 5 0.0818
#> 2 P35557 3 D 8 0.184
#> 3 P35557 4 D 8 0.147
#> 4 P35557 5 R 6 0.250
#> 5 P35557 6 A 6 0.138
#> 6 P35557 7 R 7 0.257
#> 7 P35557 8 M 9 0.142
#> 8 P35557 9 E 7 0.212
#> 9 P35557 10 A 6 0.142
#> 10 P35557 11 A 6 0.142
#> # ℹ 454 more rows
#> # ℹ 4 more variables: aa_pathogenicity_median <dbl>,
#> # aa_pathogenicity_min <dbl>, aa_pathogenicity_max <dbl>,
#> # aa_pathogenicity_mode <fct>
Coloring amino acids by position
Individual amino acids can be colored using the colorfunc= argument to r3dmol::m_style_cartoon(). This is a Javascript function that takes each atom position and returns the corresponding color. The approach taken in AlphaMissenseR is to use a template, ultimately replacing ... with a vector of residue colors.
cat(
AlphaMissenseR:::js_template("colorfunc", colors = "..."),
"\n"
)
#> function(atom) {
#> const residue_colors = [ ... ];
#> return residue_colors[atom.resi];
#> }
The function af_colorfunc_by_position() provides a mechanism for translating a vector of scores between zero and one into a vector of colors. This is illustrated for a 12-amino acid sequence where the first and last residues are uncolored.
df <- tibble(
pos = 1 + 1:10, # no color information for position 1
value = 10:1 / 10
)
colorfunc <- af_colorfunc_by_position(
df,
"pos", "value",
pos_max = 12 # no color information for position 12
)
cat(colorfunc, "\n")
#> function(atom) {
#> const residue_colors = [ 'gray', '#8E063B', '#AB5468', '#C18692', '#D2B0B6', '#DDD0D2', '#D2D3DC', '#B3B7CF', '#8C94BF', '#5D6CAE', '#023FA5', 'gray' ];
#> return residue_colors[atom.resi];
#> }
The following color function is similar to that used in af_prediction_view(), but uses the mean rather than median pathogenicity and scales the palette between the minimum and maximum values of the mean pathogenicity vector, rather than between 0 and 1.
colorfunc <-
pathogenicity |>
af_colorfunc_by_position(
"aa_pos", "aa_pathogenicity_mean",
length(pdb$seqres)
)
Add this as the colorfunc= argument to m_style_cartoon() for visualization.
r3dmol::r3dmol() |>
## use the PDB representation
r3dmol::m_add_model(r3dmol::m_bio3d(pdb)) |>
## visualize as a 'cartoon' with alpha helices and beta sheets
r3dmol::m_set_style(
style = r3dmol::m_style_cartoon(
arrow = TRUE,
## color residue according to colorfunc
colorfunc = colorfunc
)
) |>
## fit molecule into display area
r3dmol::m_zoom_to()
