JBrowseMSA

TRIM5 and the primate antiviral arms race

TRIM5 is a restriction factor: a protein that recognizes an incoming retrovirus capsid and destroys it before the virus can establish an infection. Every primate lineage carries a copy, and the copies restrict different viruses. Human TRIM5 barely restricts HIV-1, and rhesus macaque TRIM5 restricts it strongly. A single substitution (R332P) converts the human protein into an HIV-1 restrictor (Stremlau et al. 2005; Yap et al. 2005). R332 sits in a patch of the protein that has been under positive selection across the primates for tens of millions of years (Sawyer et al. 2005): as viral capsids change, the restriction factor’s binding surface changes with them, and a codon alignment records that as an excess of amino-acid-changing substitutions over silent ones.

This page builds that alignment for 32 TRIM5 orthologs, infers a tree, adds the 7-exon gene structure and a per-codon dN/dS track from HyPhy, and ends on a link that opens all of it together: the alignment, the tree, the exons, the selection track and labels on the variable patch, on the single residue R332, and on the RING/B-box zinc fingers, a built-in control that the same analysis should call conserved.

Prerequisites

Where the data comes from

1. Every ortholog NCBI knows about

datasets download gene gene-id 85363 --ortholog all \
  --include gene,cds,protein --filename trim5_orthologs.zip
unzip -o trim5_orthologs.zip -d ortho

The download holds 33 records: 28 primates and, as the outgroup, 5 rodents/squirrels (mouse relative Mus pahari, rat, and three ground squirrels). cds.fna in the package carries every predicted transcript for every gene; human alone has 21.

2. One coding sequence per species

Pick the curated (NM_) RefSeq transcript where the species has one; where every transcript is a predicted (XM_) gene model, pick whichever is closest in length to the well-annotated human isoform (1482 nt, 493 amino acids plus the stop codon). A gene with dozens of predicted splice variants offers no better criterion for choosing one transcript to align:

nm = [t for t in transcripts if t.accession.startswith('NM_')]
pool = nm if nm else transcripts
chosen = min(pool, key=lambda t: abs(len(t.seq) - 1482))

All 33 chosen sequences have a length that is a multiple of 3. Translating each and scanning for a stop codon before the last one finds one problem: the black snub-nosed monkey (Rhinopithecus bieti) has only one transcript on file, XM_017887309.1, and it translates with two internal stops, at codons 260 and 265. The species has no other transcript, so the script drops it, and 32 orthologs go forward.

3. A codon-aware alignment

docker run --rm -v "$PWD":/data -w /data \
  quay.io/biocontainers/macse:2.07--hdfd78af_0 \
  macse -prog alignSequences -seq trim5_cds.fasta \
  -out_NT trim5_macse_NT.fasta -out_AA trim5_macse_AA.fasta

MACSE aligns in amino-acid space and writes back nucleotides, so every codon occupies three whole columns and the aligner introduces no frameshift. The result has 1800 nucleotide columns (600 codons) and takes a few minutes for 32 sequences.

Every sequence’s CDS includes its own stop codon, and the gaps put each one in a different alignment column. The script therefore checks and blanks each row’s last codon individually. All 32 rows carry exactly one stop codon, in the last position.

4. A tree

docker run --rm -v "$PWD":/data -w /data \
  quay.io/biocontainers/fasttree:2.2.0--h7b50bb2_1 \
  sh -c "FastTree trim5_macse_AA.fasta > trim5.nwk"

FastTree runs on the protein translation, which carries more signal than the codons at this depth of divergence. Read the support values before relying on the topology (see also a protein family from a list of accessions). The weakest node in this tree has a local support of 0.135 and the second-weakest 0.177, both around the boundary between Old World monkey lineages, where a single gene’s tree is expected to be uncertain. The selection tests below depend on branch lengths, not on that part of the topology being right.

The 32-ortholog codon alignment with its tree, colored by nucleotide. Apes, Old World monkeys, New World monkeys, a tarsier, two strepsirrhines and the rodent/squirrel outgroup, all in one frame-correct alignment.

5. The exon structure

genestructure projects TRIM5’s MANE Select transcript, NM_033034.3, from the human reference row onto every row:

react-msaview-cli genestructure trim5_cds_aln.fasta \
  --gene-id 85363 --transcript NM_033034.3 --ref human \
  -o trim5-exons.gff

The transcript has 7 coding exons. Both zinc fingers this page uses as a control, RING and B-box, fall in the gene’s large first coding exon, and the variable patch in the SPRY domain falls in the last one.

The same alignment with the 7-exon structure overlaid, one color per exon, projected through every row’s own gaps. Exon 7, the longest, carries the whole SPRY domain.

6. Per-codon selection with HyPhy

FEL (Fixed Effects Likelihood) fits a synonymous rate (alpha) and a non-synonymous rate (beta) at every codon and tests beta against alpha with a likelihood ratio:

docker run --rm -v "$PWD":/data -w /data \
  quay.io/biocontainers/hyphy:2.5.101--h526e2cb_0 \
  hyphy fel --alignment trim5_fel_input.fasta --tree trim5.nwk \
  --output trim5_fel.json --branches All

HyPhy’s single ratio across the whole tree and every codon is 0.82, so most of TRIM5 is under mild purifying selection. FEL’s per-site test (default p<0.1, 32 taxa, all 600 codon columns, about 3 minutes with no subsampling) calls 33 of 600 sites (5.5%) positively selected and 80 (13.3%) negatively selected. The strongest positively selected site in the gene is human residue 337 (alpha near zero, beta = 23.5, p = 0.0063), inside the variable patch that step 7 zooms into.

MEME tests the same alignment for episodic rather than pervasive selection, one branch at a time per site:

docker run --rm -v "$PWD":/data -w /data \
  quay.io/biocontainers/hyphy:2.5.101--h526e2cb_0 \
  hyphy meme --alignment trim5_fel_input.fasta --tree trim5.nwk \
  --output trim5_meme.json --branches All

27 of 600 sites (p<0.05) show episodic diversifying selection, and 18 of those 27 (two-thirds) fall inside the B30.2/SPRY domain. The second-strongest MEME hit in the whole gene (p = 0.0022) is again human residue 337, agreeing with FEL. The strongest hit (p = 0.0011) is human residue 45, one of the 45 codons of the RING domain, which this page uses as the purifying control.

The track draws the per-codon beta/alpha ratio. A site with almost no synonymous substitutions puts alpha near zero, and the ratio grows without bound: at human residue 371 the optimizer returns alpha = 5098. The script therefore clamps the values to 5 and writes them as one nucleotide-column-wide JSON array (packages/app/public/data/trim5/trim5-dnds-values.json), repeated three times per codon and passed as a columnTracks bar (see data layers):

Red bars are FEL’s per-codon dN/dS, clamped at 5, above the exon-colored alignment. The spikes cluster toward the SPRY-domain half on the right and thin out over the first three exons on the left. The next two figures zoom into each end.

7. The SPRY variable patch

Sawyer et al. (2005) found a small patch of the SPRY domain carrying most of the gene’s positively selected sites, overlapping the “variable region 1” (V1) loop later work narrowed the HIV-1 restriction specificity down to (Stremlau et al. 2005; Yap et al. 2005; Ohkura et al. 2006), human residues 326-341. R332 sits in the middle of it, and R332P alone is the substitution that converts human TRIM5 into an HIV-1 restrictor:

"highlights": [
  { "row": "human", "start": 976, "end": 1023, "label": "SPRY V1 patch" },
  { "row": "human", "start": 994, "end": 996, "label": "R332" }
]

(start/end are the human row’s own nucleotide positions: residue 326 is nucleotide 976, (326-1)*3+1, and R332 is nucleotides 994-996.)

The V1 patch (orange) and R332 (red) at base resolution, under the dense red spikes of the dN/dS track. The patch is 16 codons long. FEL calls 3 of them (18.8%) positively selected and none negatively selected, the highest positive fraction and the only region of comparable size in the gene with no negatively selected site. R332 is not one of the 3: alpha = 1.77 and beta = 3.80 give an elevated ratio of about 2.1, but p = 0.52 misses FEL’s significance threshold in this 32-species set, although R332 is the residue the published point mutation targets.

8. The purifying control

RING and B-box are the zinc-finger domains that fold the protein and drive its E3 ubiquitin ligase activity, so the same analysis should call them conserved. The figure above shows both as a combined highlight, and this figure zooms into them:

RING (human aa 15-59, 45 codons) and B-box (aa 90-132, 43 codons), same color scale as every other figure on this page. FEL calls 1 RING codon and 2 B-box codons positively selected (2.2% and 4.7%), against 16 and 12 negatively selected (35.6% and 27.9%), the reverse of the V1 patch’s ratio. MEME finds no episodic hits in B-box, and its strongest hit in the whole gene is one RING codon, human residue 45.

9. Open it

This snapshot combines the last three figures: alignment, tree, exon structure, the dN/dS track and both highlights.

The full gene, with the tree and the 7-color exon structure. The RING+B-box control on the left has few dN/dS bars above it, and R332 and the SPRY V1 patch on the right sit under the densest bars.

Reproduce it end to end

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

With no arguments the script writes every file into the current directory. Point it at another gene by editing the GENE_ID and the --transcript argument to genestructure; the CDS-selection and stop-codon QC steps are generic.

See also

References