Where p53’s damaging variants fall

TP53 is the gene most often mutated in human cancer. Clinical laboratories report its damaging mutations in patients, a model predicts the effect of every possible substitution, and a laboratory screen measured thousands of variants in cells. For fifteen vertebrate p53 proteins, human to zebrafish, we:

• fetch the proteins, align them and infer a tree
• add ClinVar, AlphaMissense and a saturation screen as bar tracks
• open the three tracks from one link
• zoom to the hotspots, the N-terminal control and the oligomerization domain
• check one residue against the source files

Prerequisites

• curl, jq and awk
• MAFFT and FastTree. apt install mafft fasttree on Debian or Ubuntu, brew install mafft fasttree on macOS, or run the biocontainers with docker: quay.io/biocontainers/mafft:7.525--h031d066_1 and quay.io/biocontainers/fasttree:2.1.11--h031d066_4. The build script takes either.

The data

Sequences from NCBI RefSeq, the clinical classifications from ClinVar, the predictions from AlphaMissense as the AlphaFold entry publishes them, and one saturation screen from MaveDB. Nothing to download: the commands below fetch these.

The build script writes the domain bands in as residue numbers, taken from UniProt’s annotation of P04637.

1. Choose the rows

NCBI lists the vertebrate orthologs of TP53:

curl -s 'https://api.ncbi.nlm.nih.gov/datasets/v2alpha/gene/id/7157/orthologs?taxon_filter=vertebrates&page_size=1000' |
  jq -r '.reports[].gene | [.taxname, .common_name, .gene_id] | @tsv'

The query returns 660 species, and the list below takes fifteen. Human is NP_000537.3, the product of NM_000546, the transcript ClinVar numbers its variants on.

NP_000537.3	Human
XP_045231359.1	Macaque
XP_027375188.1	Cow
XP_047612858.1	Pig
NP_001189334.1	Horse
XP_025284307.1	Dog
XP_010594888.1	Elephant
XP_030101782.1	Mouse
NP_112251.2	Rat
XP_056673571.1	Opossum
NP_990595.1	Chicken
XP_065426476.1	Turtle
XP_062840611.1	Anole
NP_001001903.1	Frog
XP_073805887.1	Zebrafish

Save that as accessions.tsv, tab separated.

2. Fetch the sequences

One request takes the whole list, and an awk pass renames each record from its accession to its label:

IDS=$(cut -f1 accessions.tsv | paste -sd,)
curl -sf "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&id=$IDS&rettype=fasta&retmode=text" \
  > p53-refseq.fa

# every layer below names a row by the label the viewer draws
awk 'NR == FNR {split($0, row, "\t"); label[row[1]] = row[2]; next}
  /^>/ {split($0, field, " "); print ">" label[substr(field[1], 2)]; next}
  NF {print}' accessions.tsv p53-refseq.fa > p53.fasta

Check the lengths before aligning them:

awk '/^>/ {if (name) print name, n; name = $0; n = 0; next} {n += length($0)}
  END {print name, n}' p53.fasta

>Human 393
>Macaque 393
>Cow 386
>Pig 386
>Horse 381
>Dog 381
>Elephant 390
>Mouse 390
>Rat 391
>Opossum 361
>Chicken 367
>Turtle 409
>Anole 395
>Frog 362
>Zebrafish 374

The lengths run from 361 to 409, with no outlier that would mark a truncated model.

3. Align them and infer a tree

# --anysymbol accepts a record with an X or a U in it
mafft --auto --anysymbol p53.fasta > p53.afa

# -lg is the amino-acid substitution model; FastTree reads the alignment and
# writes Newick with a support value on each internal node
FastTree -lg -quiet p53.afa > p53.nh

Over the 478 alignment columns, FastTree joins the two rodents at support 1.000 and human with macaque at 0.997, but its deep splits reach only 0.295 and 0.310, so the tree orders the rows for reading and its deep topology is unsupported.

4. Count ClinVar’s pathogenic variants per residue

The script reads the classification from each record and keeps the pathogenic and likely pathogenic ones, because a search term matching “pathogenic” also matches “Conflicting classifications of pathogenicity”:

EUTILS=https://eutils.ncbi.nlm.nih.gov/entrez/eutils
TERM='TP53[gene] AND "missense variant"[molecular consequence]'
curl -sf "$EUTILS/esearch.fcgi?db=clinvar&retmode=json&retmax=5000&term=$(
  jq -rn --arg t "$TERM" '$t | @uri'
)" | jq -r '.esearchresult.idlist[]' > clinvar.ids

# esummary takes 200 ids at a time. Keep only what the submitters classify as
# pathogenic, and take the variant name, which carries the protein change
split -l 200 clinvar.ids clinvar.batch.
for batch in clinvar.batch.*; do
  curl -sf "$EUTILS/esummary.fcgi?db=clinvar&retmode=json&id=$(paste -sd, "$batch")" |
    jq -r '.result | del(.uids) | .[]
      | select(.germline_classification.description
        | IN("Pathogenic", "Likely pathogenic", "Pathogenic/Likely pathogenic"))
      | .variation_set[].variation_name' >> clinvar.changes
  # NCBI asks unauthenticated callers to stay under three requests a second
  sleep 0.4
done

Of 1,519 missense records, 253 pass the filter. The count per residue reads p.Arg248Gln out of NM_000546.6(TP53):c.743G>A (p.Arg248Gln) and drops the nonsense changes with grep -v 'Ter)$', since a stop truncates the protein instead of exchanging one residue:

grep '^NM_000546' clinvar.changes |
  grep -oE '\(p\.[A-Z][a-z]{2}[0-9]+[A-Z][a-z]{2}\)' |
  grep -v 'Ter)$' |
  grep -oE '[0-9]+' |
  awk -v n=393 '{count[$1]++}
    END {for (i = 1; i <= n; i++) printf "%s%d", (i > 1 ? "," : ""), count[i]}'

The 253 variants fall on 104 of the 393 residues, 94% of them inside the DNA-binding domain.

5. Average AlphaMissense scores per residue

AlphaMissense scores every possible substitution, and the AlphaFold entry for an accession names its file:

UA='msaview-tutorial (https://gmod.org/JBrowseMSA)'
AM_URL=$(curl -sf -A "$UA" https://alphafold.ebi.ac.uk/api/prediction/P04637 |
  jq -r '.[0].amAnnotationsUrl')
curl -sf -A "$UA" "$AM_URL" -o alphamissense.csv

The command below averages the 19 scores, each 0 to 1, at each of the 393 residues, reading the position out of the M1A in the first column:

tail -n +2 alphamissense.csv |
  awk -F, -v n=393 '{
    pos = substr($1, 2, length($1) - 2)
    sum[pos] += $2
    seen[pos]++
  }
  END {for (i = 1; i <= n; i++)
    printf "%s%.2f", (i > 1 ? "," : ""), (seen[i] ? sum[i] / seen[i] : 0)}'

6. Average the saturation screen per residue

Giacomelli and colleagues grew A549 cells with p53 variants beside a wild-type copy under nutlin-3, so a variant that enriches has disabled p53 with the wild-type protein present:

curl -sfL https://api.mavedb.org/api/v1/score-sets/urn:mavedb:00000068-a-1/scores \
  -o mavedb.csv

# the score column against the protein change; p.Met384= is synonymous and
# p.Arg248Ter is a stop, and neither is one residue exchanged for another
tail -n +2 mavedb.csv |
  awk -F, -v n=393 '$4 ~ /^p\.[A-Z][a-z][a-z][0-9]+[A-Z][a-z][a-z]$/ &&
    $4 !~ /Ter$/ && $5 != "NA" {
    match($4, /[0-9]+/)
    pos = substr($4, RSTART, RLENGTH)
    sum[pos] += $5
    seen[pos]++
  }
  END {for (i = 1; i <= n; i++) {
    mean = seen[i] ? sum[i] / seen[i] : 0
    printf "%s%.2f", (i > 1 ? "," : ""), (mean > 0 ? mean : 0)
  }}'

The screen used a P72R background, so its residue 72 is an arginine where the aligned RefSeq protein has a proline. A bar track clamps at zero, so about half the residues draw nothing.

7. Open the three tracks from one link

Each columnTracks entry names the row its values index, so the viewer maps human residue numbers onto the alignment columns:

jq -rj --arg msa data/p53/p53-vertebrates.afa \
  --arg tree data/p53/p53-vertebrates.nh '
  {
    msaview: {
      type: "MsaView",
      relativeTo: "Human",
      colorSchemeName: "clustalx_protein_dynamic",
      msaFilehandle: {uri: $msa},
      treeFilehandle: {uri: $tree},
      highlights: .highlights,
      columnTracks: .columnTracks
    }
  } | @uri' p53-layers.json

The whole snapshot travels in the link, and the server stops serving a request line past 8,192 characters. This link stays at 7,820 because AlphaMissense goes in as a percent against domainMax: 100.

WebR

ggmsaview passes the highlights and columnTracks of p53-layers.json to msaview() as they stand, and msaview_plot() draws a ggplot2 figure:

library(ggmsaview)
data <- "https://gmod.org/JBrowseMSA/demo/data/p53/"
layers <- jsonlite::fromJSON(paste0(data, "p53-layers.json"),
  simplifyVector = FALSE)
w <- msaview(
  msa = paste0(data, "p53-vertebrates.afa"),
  tree = paste0(data, "p53-vertebrates.nh"),
  relative_to = "Human",
  color_scheme = "clustalx_protein_dynamic",
  highlights = layers$highlights,
  column_tracks = layers$columnTracks
)
w
ggplot2::ggsave("p53-variant-three-tracks-r.png",
  msaview_plot(w, tracks = c("clinvar", "alphamissense", "mavedb",
    "conservation")),
  width = 12, height = 7)

Averaged over the regions UniProt annotates:

8. Zoom to the hotspots

The six residues most often mutated in cancer all lie in the DNA-binding domain:

The rows with the human residue at each hotspot, beside three N-terminal residues:

9. Compare the N-terminal control

UniProt annotates the first hundred residues as disordered. The figure below shows them at the hotspot figure’s zoom:

10. Zoom to the oligomerization domain

Across the 32 residues of the helix that forms the p53 tetramer, ClinVar has 13 pathogenic missense variants and AlphaMissense averages 0.59, while the screen comes out at zero, because every raw mean there is zero or negative:

11. Check a residue against the files

Every bar is a count or a mean over a source file, so a grep recovers it. R248 against P47:

grep 'p.Arg248' clinvar.changes                          # 6 variants
grep -E '^R248' alphamissense.csv |
  awk -F, '{s += $2; n++} END {printf "%.2f over %d\n", s / n, n}'   # 1.00 over 19
grep -E 'p\.Arg248[A-Z]' mavedb.csv | grep -v Ter |
  awk -F, '{s += $5; n++} END {printf "%.2f over %d\n", s / n, n}'   # 1.50 over 19

At P47 the same commands return 0 variants, 0.15 over 19, and a negative MaveDB mean, which the build script writes as 0.

Two of the six ClinVar records at 248, c.742C>T and c.741_742delinsTT, reach the same p.Arg248Trp, so the track counts records rather than patients.

Reproduce it end to end

build_p53_variant_effects.sh runs the steps above in order:

1. It fetches the fifteen proteins from NCBI, aligns them with MAFFT and infers a tree with FastTree.
2. It counts the pathogenic and likely pathogenic ClinVar missense alleles at each residue of the human row.
3. It averages the AlphaMissense and screen scores at each residue.
4. It prints the hotspot and region tables and writes the layers file and the link.

The script needs the tools under Prerequisites.

curl -O https://raw.githubusercontent.com/GMOD/JBrowseMSA/main/docs/tutorials/scripts/build_p53_variant_effects.sh
bash build_p53_variant_effects.sh out/

The script leaves the link in out/p53-link.url and runs the biocontainers when MAFFT and FastTree are not on PATH. ClinVar updates weekly, so its counts change between runs where the other sources do not.

Here is what that link opens, against the hosted copy of this alignment.

See also

• Data layers
• A protein family from a list of accessions
• AlphaFold confidence across a protein family
• User guide
• CLI

References

• Landrum MJ, et al. ClinVar: improving access to variant interpretations and supporting evidence. Nucleic Acids Research 46:D1062-D1067 (2018).
• Cheng J, et al. Accurate proteome-wide missense variant effect prediction with AlphaMissense. Science 381:eadg7492 (2023).
• Varadi M, et al. AlphaFold Protein Structure Database in 2024. Nucleic Acids Research 52:D368-D375 (2024), which serves the AlphaMissense file above.
• Giacomelli AO, et al. Mutational processes shape the landscape of TP53 mutations in human cancer. Nature Genetics 50:1381-1387 (2018), the screen behind urn:mavedb:00000068-a-1.
• Esposito D, et al. MaveDB: an open-source platform to distribute and interpret data from multiplexed assays of variant effect. Genome Biology 20:223 (2019).
• Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7. Molecular Biology and Evolution 30:772-780 (2013).
• Price MN, Dehal PS, Arkin AP. FastTree 2, approximately maximum-likelihood trees for large alignments. PLoS ONE 5:e9490 (2010).