JBrowseMSA

A protease family in R

Chymotrypsin cuts a peptide bond with three residues, His57, Asp102 and Ser195, that sit far apart in the chain and side by side in the folded protein. Several human proteins have the same peptidase S1 domain with one of those three replaced, and cut nothing. Working in R with ggmsaview from fourteen UniProtKB accessions, we:

  • cut each chain to its peptidase S1 domain and align the domains
  • build a tree from the alignment
  • add a BLOSUM62 track, disulfide arcs and bands on the three catalytic columns
  • read the letter each row has in those columns
  • save the widget and write a link for the web viewer

Prerequisites

  • R 4.3 or newer
  • ggmsaview: remotes::install_github("GMOD/ggmsaview")
  • Bioconductor: BiocManager::install(c("Biostrings", "DECIPHER"))
  • CRAN: install.packages(c("ape", "jsonlite", "htmlwidgets"))

The data

UniProtKB supplies the feature table, with each domain’s bounds and charge-relay residues, and the sequences. Nothing to download: the commands below fetch these.

What it isFile
The feature table for all fourteen accessions in one request, shown here with two of themuniprotkb-accessions
The sequences, from the same endpointuniprotkb-accessions
The alignment the commands below write, which the link in step 10 opensproteases.aln
Its treeproteases.nh
The files this page writes, hosted for the figure links
What it isFile
The layers R computedproteases-layers.json

1. List the proteins

UniProt annotates the first nine as active peptidases and the last five with the charge-relay system broken, which the page checks at the end.

proteins <- data.frame(
  accession = c(
    "P00766", "P07477", "P08246", "P10144", "P07288",
    "P00734", "P00742", "P00747", "P00750",
    "P00738", "P20160", "P22891", "P14210", "P26927"
  ),
  label = c(
    "Chymotrypsinogen", "Trypsin_1", "Elastase", "Granzyme_B", "PSA",
    "Prothrombin", "Factor_X", "Plasminogen", "tPA",
    "Haptoglobin", "Azurocidin", "Protein_Z", "HGF", "MST1"
  )
)

The label becomes the FASTA header, the tree tip and the row the viewer draws down the side. Underscores keep the names through Newick, which reads a space as the end of a name.

2. Fetch the feature table

One TSV request covers all fourteen entries and three feature types.

REST <- "https://rest.uniprot.org/uniprotkb/accessions?accessions="
ids <- paste(proteins$accession, collapse = ",")

fields <- "accession,id,protein_name,ft_domain,ft_act_site,ft_disulfid"
features <- read.delim(
  url(paste0(REST, ids, "&fields=", fields, "&format=tsv")),
  quote = "", check.names = FALSE
)
features <- features[match(proteins$accession, features$Entry), ]

Two small parsers turn each feature string into numbers:

# "DOMAIN 364..618; /note="Peptidase S1"" -> c(364, 618)
s1_domain <- function(text) {
  m <- regmatches(text, regexpr('DOMAIN (\\d+)\\.\\.(\\d+); /note="Peptidase S1"', text))
  as.integer(regmatches(m, gregexpr("\\d+", m))[[1]])
}
# "ACT_SITE 406; ... ACT_SITE 462; ..." -> c(406, 462, 568)
positions <- function(text, key) {
  m <- regmatches(text, gregexpr(paste0(key, " (\\d+)"), text))[[1]]
  as.integer(sub(paste0(key, " "), "", m))
}

domains <- lapply(features[["Domain [FT]"]], s1_domain)
proteins$start <- vapply(domains, `[`, integer(1), 1)
proteins$end <- vapply(domains, `[`, integer(1), 2)
proteins$act <- lapply(features[["Active site"]], positions, "ACT_SITE")
          protein   domain residues   active_site
 Chymotrypsinogen  16..243      228  57, 102, 195
        Trypsin_1  24..244      221  63, 107, 200
         Elastase  30..247      218  70, 117, 202
       Granzyme_B  21..245      225  64, 108, 203
              PSA  25..258      234  65, 120, 213
      Prothrombin 364..618      255 406, 462, 568
         Factor_X 235..467      233 276, 322, 419
      Plasminogen 581..808      228 622, 665, 760
              tPA 311..561      251 357, 406, 513
      Haptoglobin 162..404      243          none
       Azurocidin  27..244      218          none
        Protein_Z 175..400      226          none
              HGF 495..721      227          none
             MST1 484..709      226          none

The domain starts at residue 16 of chymotrypsinogen and at residue 581 of plasminogen, so residue 195 is a different place in each row.

3. Cut each chain to its domain

fasta <- "proteases-full.fasta"
download.file(paste0(REST, ids, "&format=fasta"), fasta, quiet = TRUE)
full <- readAAStringSet(fasta)
names(full) <- sub("^[a-z]+\\|([A-Z0-9]+)\\|.*$", "\\1", names(full))
full <- full[proteins$accession]

s1 <- subseq(full, start = proteins$start, end = proteins$end)
names(s1) <- proteins$label
full chains 245-810 aa, S1 domains 218-255 aa

The cut removes the Gla and kringle domains that prothrombin and plasminogen have ahead of the peptidase.

4. Align the domains

aln <- AlignSeqs(s1, verbose = FALSE)
writeXStringSet(aln, "proteases.aln")
alignment: 14 rows x 297 columns
The fourteen domains after AlignSeqs in accession-table order, the nine active peptidases above the five without a charge-relay system. The gaps fall in the loops, which differ in length between families.

5. Build a tree

nj builds the tree from the alignment’s distances, and ladderize sorts each node so the smaller clade comes first.

tree <- ladderize(nj(as.dist(DistanceMatrix(aln, verbose = FALSE))))
write.tree(tree, "proteases.nh")

An ape::phylo goes into msaview() as it is, along with the AAStringSet:

msaview(msa = aln, tree = tree, color_scheme = "clustalx_protein_dynamic")
The tree beside the alignment, with the five rows that have no charge-relay system outlined. Haptoglobin and protein Z sit with prothrombin and factor X, azurocidin pairs with elastase, and HGF pairs with MST1 beside plasminogen, so the five fall in four places.

6. Score each column with BLOSUM62

Mean pairwise BLOSUM62 weights agreement in a column by how rare the residue is, so a column of conserved cysteines scores 9 and one of conserved alanines 4.

data(BLOSUM62, package = "Biostrings")
m <- as.matrix(aln)
pairs <- combn(nrow(m), 2)
blosum <- apply(m, 2, function(column) {
  a <- column[pairs[1, ]]
  b <- column[pairs[2, ]]
  keep <- a %in% rownames(BLOSUM62) & b %in% rownames(BLOSUM62)
  if (!any(keep)) 0 else mean(BLOSUM62[cbind(a[keep], b[keep])])
})

column_tracks takes the vector as one bar per column, scaled by scales$y$domainMax:

blosum_track <- list(
  id = "blosum", name = "Mean pairwise BLOSUM62", kind = "bar",
  values = pmax(blosum, 0), scales = list(y = list(domainMax = 9)),
  color = "#6a51a3", height = 60
)

7. Draw the disulfide bonds as arcs

The feature table lists five disulfide bonds in chymotrypsinogen, and four lie inside the aligned domain.

disulfides <- function(text) {
  m <- regmatches(text, gregexpr("DISULFID (\\d+)\\.\\.(\\d+)", text))[[1]]
  matrix(as.integer(unlist(regmatches(m, gregexpr("\\d+", m)))), ncol = 2, byrow = TRUE)
}

ref <- which(proteins$label == "Chymotrypsinogen")
ss <- disulfides(features[["Disulfide bond"]][ref])
ss <- ss[ss[, 1] >= proteins$start[ref] & ss[, 2] <= proteins$end[ref], , drop = FALSE]
bonds <- data.frame(
  start = ss[, 1] - proteins$start[ref] + 1,
  end = ss[, 2] - proteins$start[ref] + 1
)

ss_track <- list(
  id = "ss", name = "Disulfide bonds (chymotrypsinogen)", kind = "arc",
  arcs = bonds, row = "Chymotrypsinogen", color = "#b35806", height = 50
)

The subtraction converts a chain position to a position in the cut domain, and row = "Chymotrypsinogen" has the viewer walk that row’s gaps to find each column.

The alignment with purple bars of mean pairwise BLOSUM62 above it and four arcs joining the columns of chymotrypsinogen’s bonded cysteines. The tallest bar, 11, stands on a column every row reads as tryptophan, and two of the arc ends stand at 9, which is a column of cysteines.

8. Band the catalytic columns

Subtracting the domain start moves the charge-relay residues to domain numbering too, and the viewer projects them onto columns:

triad <- proteins$act[[ref]] - proteins$start[ref] + 1
names(triad) <- c("His57", "Asp102", "Ser195")

bands <- Map(
  function(residue, label) {
    list(row = "Chymotrypsinogen", start = residue, end = residue, label = label)
  },
  triad, names(triad)
)

widget <- msaview(
  msa = aln,
  tree = tree,
  color_scheme = "clustalx_protein_dynamic",
  column_tracks = list(blosum_track, ss_track),
  highlights = bands,
  height = 520
)
widget

9. Read each row’s catalytic residues

A cumulative sum over the reference row finds the three banded columns in the matrix, the same walk the viewer does:

residue_column <- cumsum(m[ref, ] != "-")
triad_column <- vapply(triad, function(r) which(residue_column == r)[1], integer(1))
data.frame(
  protein = proteins$label,
  charge_relay = ifelse(lengths(proteins$act) > 0, "annotated", "none"),
  His57 = m[, triad_column[1]],
  Asp102 = m[, triad_column[2]],
  Ser195 = m[, triad_column[3]]
)
          protein charge_relay His57 Asp102 Ser195
 Chymotrypsinogen    annotated     H      D      S
        Trypsin_1    annotated     H      D      S
         Elastase    annotated     H      D      S
       Granzyme_B    annotated     H      D      S
              PSA    annotated     H      D      S
      Prothrombin    annotated     H      D      S
         Factor_X    annotated     H      D      S
      Plasminogen    annotated     H      D      S
              tPA    annotated     H      D      S
      Haptoglobin         none     K      D      A
       Azurocidin         none     S      D      G
        Protein_Z         none     K      D      M
              HGF         none     Q      D      Y
             MST1         none     Q      Q      Y

All nine active peptidases read H, D and S. Each of the other five reads something else at His57 and at Ser195, and four of them keep the aspartate.

The alignment at the catalytic serine, 13 pixels per column, with the five rows boxed. Nine rows read the GDSGGP motif; haptoglobin reads A, protein Z reads M, MST1 and HGF read Y, and azurocidin reads G.

The BLOSUM62 track, the arcs and the bands go into a URL that opens the hosted alignment in the web viewer:

snapshot <- list(msaview = list(
  type = "MsaView",
  height = 560,
  treeAreaWidth = 150,
  colorSchemeName = "clustalx_protein_dynamic",
  msaFilehandle = list(uri = "data/proteases/proteases.aln"),
  treeFilehandle = list(uri = "data/proteases/proteases.nh"),
  columnTracks = list(blosum_track, ss_track),
  highlights = unname(bands)
))
json <- as.character(jsonlite::toJSON(snapshot, auto_unbox = TRUE, digits = 3))
paste0("https://gmod.org/JBrowseMSA/demo/#data=", URLencode(json, reserved = TRUE))
the link is 3075 characters

unname() is required, because Map over a named vector names its result and a named list serializes as a JSON object where the viewer reads an array.

saveWidget writes a page that opens without R, and msaview_plot() draws the widget as a ggplot2 figure. In a fresh session, the hosted alignment, tree and proteases-layers.json rebuild widget:

library(ggmsaview)
data <- "https://gmod.org/JBrowseMSA/demo/data/proteases/"
layers <- jsonlite::fromJSON(paste0(data, "proteases-layers.json"),
  simplifyVector = FALSE)
widget <- msaview(
  msa = paste0(data, "proteases.aln"),
  tree = paste0(data, "proteases.nh"),
  color_scheme = "clustalx_protein_dynamic",
  column_tracks = layers$columnTracks,
  highlights = layers$highlights
)
htmlwidgets::saveWidget(widget, "proteases.html", selfcontained = FALSE)
ggplot2::ggsave("r-protease-triad-r.png",
  msaview_plot(widget, tracks = c("blosum", "ss")),
  width = 12, height = 6)
The fourteen S1 domains over 297 columns in ClustalX colors, under the BLOSUM62 bars and the chymotrypsinogen disulfide arcs, with the His57, Asp102 and Ser195 columns boxed and the tree on the left.

Reproduce it end to end

build_r_protease_triad.R runs the R above in order:

  1. It reads the domain bounds, charge-relay residues and disulfide bonds from UniProt in one request, then cuts each chain to its S1 domain.
  2. DECIPHER aligns the domains, and ape builds a neighbor-joining tree.
  3. It scores every column by mean pairwise BLOSUM62 and prints each row’s letters at His57, Asp102 and Ser195 beside UniProt’s annotation.
  4. It draws the result with msaview(), saves the widget, and writes the layers as JSON and as a URL for the web viewer.

It needs the tools under Prerequisites.

curl -O https://raw.githubusercontent.com/GMOD/JBrowseMSA/main/docs/tutorials/scripts/build_r_protease_triad.R
Rscript build_r_protease_triad.R out/

The script writes out/proteases.aln, out/proteases.nh, out/proteases-layers.json and out/proteases.html.

See also

References

  • Blow DM, Birktoft JJ, Hartley BS. Role of a buried acid group in the mechanism of action of chymotrypsin. Nature 221:337-340 (1969).
  • Rawlings ND, et al. The MEROPS database of proteolytic enzymes, their substrates and inhibitors in 2017. Nucleic Acids Research 46:D624-D632 (2018), which classifies the five rows above as non-peptidase homologues.
  • Wright ES. DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics 16:322 (2015).
  • Paradis E, Schliep K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35:526-528 (2019).
  • Henikoff S, Henikoff JG. Amino acid substitution matrices from protein blocks. PNAS 89:10915-10919 (1992).
  • The UniProt Consortium. UniProt: the Universal Protein Knowledgebase in 2025. Nucleic Acids Research 53:D609-D617 (2025).