JBrowseMSA

A protein family from a list of accessions

Human NLRP1, an immune sensor, starts with about eighty residues, a pyrin domain, that the mouse protein lacks. Aligning the two with their relatives puts the shared domains in the same columns, and the pyrin domain’s absence shows as blank space in the rows without it. This page starts from twelve UniProt accessions for NLRP1 and ends on a link that opens the twelve sequences aligned, with the tree inferred from them and their Pfam domains drawn in the alignment’s coordinates. Four commands produce the three files.

Prerequisites

Every figure below links to the live view it captured.

Where the data comes from

Twelve NLRP1 orthologs as UniProtKB holds them, with Pfam matches from InterPro release 110.0.

Both endpoints serve one protein per request, so the whole family takes twelve requests to each.

The family

NLRP1 is an inflammasome sensor. Every vertebrate ortholog shares the same core in the same order: a NACHT nucleotide-binding domain, a winged helix, a helical domain, then FIIND and CARD at the C terminus. The orthologs differ at the N terminus, where primates carry a pyrin (PYD) death-fold domain and rodents do not.

The proteins run from 1143 to 1537 residues, so residue 328 is a different place in every row. Drawing the domains in each protein’s own coordinates puts the shared core in twelve different places; drawing them in the alignment’s coordinates puts it in one.

1. Name the rows

Write one accession per line, followed by the label you want down the side of the viewer. The label becomes the FASTA defline, the tree tip, and the first column of the domain GFF.

Q9C000	Human
H2QC06	Chimp
A0A1D5QWR0	Rhesus
A0A8I3MJ75	Dog
A0ABM3YGI7	Hedgehog
E1BNN6	Cow
K9IW94	Pig
A0A9L0RFW2	Horse
Q2LKU9	Mouse
D9I2G4	Rat
A0ABM2XMM7	Hamster
A0A386CAB9	Zebrafish

Save that as accessions.tsv. The separator is a tab, and a line starting with # is a comment, so the file can carry its own notes about why each row is in it.

2. Fetch the sequences

# read each accession and write one FASTA record named by its label, which is
# what the viewer draws down the side
while IFS=$'\t' read -r accession label; do
  printf '>%s\n' "$label" >> family.fasta
  curl -sf "https://rest.uniprot.org/uniprotkb/$accession.fasta" |
    tail -n +2 | tr -d '\n' >> family.fasta
  printf '\n' >> family.fasta
done < accessions.tsv

tail -n +2 drops UniProt’s own defline, which carries the accession, the species and the protein name; tr -d '\n' unwraps the sequence onto one line. Check what you got before aligning it:

awk '/^>/{name=$0; next}{print name, length($0)}' family.fasta
>Human 1473
>Chimp 1471
>Rhesus 1537
>Dog 1432
>Hedgehog 1360
>Cow 1410
>Pig 1294
>Horse 1465
>Mouse 1182
>Rat 1218
>Hamster 1143
>Zebrafish 1355

All twelve records have plausible lengths for a full-length NLRP1. A truncated fragment or a stray isoform would show up here as an outlier length.

The viewer opens unaligned sequences as readily as an alignment, once every row is the same length. Right-pad the shorter ones and family.fasta loads as a block of twelve rows.

The twelve sequences before alignment, one row each, colored by residue. Every row starts at residue 1 and stops at its own length, which is the ragged right edge. The conservation track above them is near flat.

3. Align them

clustalw -INFILE=family.fasta -ALIGN -TYPE=PROTEIN \
  -OUTPUT=FASTA -OUTFILE=family.afa

The alignment takes under two seconds. ClustalW reports Alignment Score 258590 and 1666 columns, 129 more than the longest input, to make room for insertions.

ClustalW is a progressive aligner. It is fast and deterministic, needs no configuration, and places the domain architecture below in the right columns. For a figure about the phylogeny, use MAFFT or MUSCLE for the alignment and IQ-TREE or RAxML for the tree.

The same twelve after ClustalW, at 1666 columns. Vertical bands of color run through every row where the aligner found the same residues, and the pale stretches are gaps it inserted to keep them there. The conservation track has structure now.

4. Infer a tree

# -TREE reads an existing alignment and writes neighbor-joining Newick
clustalw -INFILE=family.afa -TREE -TYPE=PROTEIN -OUTPUTTREE=phylip

# ClustalW wraps the Newick across many lines; the viewer wants one string
tr -d '[:space:]' < family.ph > family.nwk

Read the result before trusting it:

(((((Mouse:0.11672,Rat:0.11797):0.05644,Hamster:0.14779):0.05664,
Zebrafish:0.56540):0.02478,Hedgehog:0.22128):0.01141,...

Zebrafish sits inside the rodents. It should be the outgroup, and its branch length is 0.56540 against 0.11 to 0.22 for everything else, the longest branch in the tree by a factor of three. Neighbor joining pulls the longest branch towards whichever other branch is longest (long-branch attraction). This page uses the tree to order the rows; its topology is not a phylogenetic result.

With the tree open beside the alignment, the viewer draws the rows in tree order.

The alignment with the ClustalW tree drawn beside it. The three primates are adjacent rows and so are the three rodents, and Zebrafish is drawn next to the rodents on the long branch that put it there.

Clicking a node in the tree collapses the clade under it into one triangle, labelled with how many tips it holds, and hides those rows from the alignment.

Mouse, Rat and Hamster collapsed into the triangle marked 3. Eleven rows are drawn where there were twelve, and the columns of the rows that remain do not move.

5. Ask InterPro what the domains are

InterPro precomputes matches for every UniProtKB sequence, so for accessions the CLI looks them up:

react-msaview-cli interpro accessions.tsv -o family-domains.gff
InterPro release 110.0; reading precomputed pfam matches...
  [1/12] Q9C000: 9 pfam entries
  ...
12 fetched, 0 from /home/you/.cache/react-msaview-cli/interpro

The lookup takes four seconds for twelve proteins. The CLI writes the release number into the GFF header, which pins the coordinates in the figures to one InterPro version. A rerun reads from the disk cache and makes a single request.

The other command, react-msaview-cli interproscan, submits sequences to the EBI job queue and waits about fifteen minutes. Use it for de novo assemblies or predicted proteins that InterPro has never seen.

Count what came back before drawing it:

grep IPR004020 family-domains.gff | cut -f1
Human
Chimp
Rhesus
Dog
Hedgehog

Five rows out of twelve carry the PYD, at residues 9-83 in each.

6. Open the three files

family.afa, family.nwk and family-domains.gff go in the alignment, tree and annotation slots of the import form, either as local files or as URLs. Each slot has a FILE and a URL toggle, and the GFF slot is the one marked optional.

All three files in one view. The domain boxes replace the residue colors: NACHT first, the winged helix and the helical domain next to it, the leucine-rich repeats scattered through the middle, then the two FIIND blocks and the CARD at the right. Each of those spans one band of columns across all twelve rows. The PYD at the far left is drawn on five rows, and the other seven are blank there.

The legend in the top right names every accession in the file, and File → Annotations → Filter annotations opens the same list with a checkbox and a count per accession.

The filter dialog over the view it filters, one row per InterPro accession, in the color the overlay draws it. IPR004020, the pyrin domain, has a count of 5 where NACHT, the winged helix, the helical domain and the CARD each have 12. IPR025307 counts 24 because it matches twice in every row, once as FIIND and once as the UPA-FIIND block beside it.

Rows with the pyrin domain

Unchecking everything except IPR004020 leaves only the pyrin domain on screen.

The overlay reduced to the pyrin domain. The three primates carry it as one block of rows, Dog carries it in the row right under them, and Hedgehog carries it in the last row of the frame. Cow, Pig and Horse sit between those two without it.

Domain coordinates and columns

The NACHT domain is annotated in all twelve rows, and InterPro reports it at a different residue in each:

grep IPR007111 family-domains.gff | cut -f1,4,5
Human	328	497
Chimp	328	497
Rhesus	332	501
Dog	299	467
Hedgehog	296	465
Cow	294	463
Pig	162	330
Horse	328	493
Mouse	134	301
Rat	176	344
Hamster	93	258
Zebrafish	258	428

Human 328 and Hamster 93 are 235 residues apart. In the viewer, Human, Chimp, Rhesus, Dog, Hedgehog, Cow, Pig and Horse all begin that domain in column 371, and Mouse, Rat, Hamster and Zebrafish begin it in column 372, because the overlay projects each row’s residue coordinates through that row’s own gaps before drawing. The same projection draws the PYD in one place in the five rows that carry it, and the seven rows without it are empty there.

Missing annotations

Cow, Hedgehog and Zebrafish have no leucine-rich repeats called at all, where Human has three and Rhesus four:

grep IPR001611 family-domains.gff | cut -f1 | sort | uniq -c
      2 Chimp
      2 Dog
      1 Hamster
      2 Horse
      3 Human
      1 Mouse
      3 Pig
      2 Rat
      4 Rhesus

A region that draws blank can mean the protein has none, or that nobody annotated the ones it has. Several of these rows are unreviewed entries whose gene models under-call. At whole-protein scale the two look identical; at base resolution a row with residues under an unannotated stretch looks different from a row that is gap there.

The left edge of the PYD block at one letter per column. The five rows with the domain are green. Cow fills all 75 of those columns with residues and carries no pyrin domain over them; Horse fills 73 of them. Mouse, Rat and Hamster are mostly dashes there, and Hamster has 24 residues in the 75 columns.

Share the view

A file you opened from your own computer travels inside the page URL, and an alignment this size does not fit, so the header says Not in the link. Put the three files on a server that sends them over HTTPS with CORS (a GitHub repo through raw.githubusercontent.com works), and name the URLs instead:

{
  "msaview": {
    "type": "MsaView",
    "colWidth": 0.7,
    "colorSchemeName": "clustalx_protein_dynamic",
    "msaFilehandle": { "uri": "https://example.org/family.afa" },
    "treeFilehandle": { "uri": "https://example.org/family.nwk" },
    "gffFilehandle": { "uri": "https://example.org/family-domains.gff" },
    "highlights": [{ "start": 39, "end": 113, "label": "PYD, Pfam PF02758" }]
  }
}

URL-encode that and append it to the app as ?data=. The link carries the addresses of the files, so it works at any size. The optional highlights takes 1-based inclusive columns and puts the PYD band in the link.

What that link opens: the three hosted files, and columns 39 to 113 banded and labelled PYD across all twelve rows.

Reproduce it end to end

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

With no arguments the script writes the accession list above and builds all three files beside it. Point it at your own list to do the same for another family:

bash build_protein_family.sh my-accessions.tsv out/

See also

References