JBrowseMSA

TRIM5 and the primate antiviral arms race

TRIM5 is a primate restriction factor that binds incoming retrovirus capsids, and a single residue, R332, decides whether the human copy restricts HIV-1 (Stremlau et al. 2005; Yap et al. 2005). Positive selection has reshaped its capsid-binding patch across the primates (Sawyer et al. 2005), and a codon alignment records that selection as an excess of amino-acid-changing substitutions over silent ones. For 32 TRIM5 orthologs, primates with a rodent outgroup, we:

  • download the orthologs and keep one coding sequence per species
  • align the coding sequences codon by codon and infer a tree
  • add the 7-exon gene structure
  • add a per-codon dN/dS track from HyPhy
  • label the variable patch, R332 and the conserved RING/B-box control, and open the result from one link

Prerequisites

The data

TRIM5 orthologs from NCBI Gene, as the Datasets v2 API serves them for human Gene ID 85363. Nothing to download: the commands below fetch these.

What it isFile
The orthologs with their genes, coding sequences and proteins, fetched with the NCBI Datasets CLI85363-download

The script takes the domain boundaries on the human protein (UniProtKB Q9C035) from InterPro’s PROSITE profile calls: RING (PS50089, aa 15-59), B-box (PS50119, aa 90-132) and B30.2/SPRY (PS50188, aa 281-493).

The files this page writes, hosted for the figure links
What it isFile
The codon alignment the commands below writetrim5-cds.stock
Its exon structuretrim5-exons.gff

1. Download the orthologs

NCBI Datasets downloads every ortholog of human Gene ID 85363:

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 an outgroup of 5 rodents and squirrels.

2. Pick one coding sequence per species

cds.fna lists every predicted transcript of every gene, 21 for human alone. Keep the curated (NM_) RefSeq transcript where the species has one, and otherwise the predicted (XM_) model closest in length to the human isoform (1482 nt, 493 amino acids plus the stop codon):

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 lengths are multiples of 3, but the black snub-nosed monkey (Rhinopithecus bieti) has one transcript on file, XM_017887309.1, which translates with internal stops at codons 260 and 265. The script drops it, and 32 orthologs go forward.

3. Align the codons

MACSE aligns in amino-acid space and writes back nucleotides, so every codon fills three whole columns:

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

The result has 1800 nucleotide columns (600 codons). Each row ends on its stop codon, and the gaps put each stop in a different column, so the script blanks the last codon of each row separately.

4. Infer a tree

FastTree infers the tree from the protein translation, which holds more signal than the codons at this depth of divergence:

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"

The weakest node has a local support of 0.135 and the second-weakest 0.177, both between Old World monkey lineages.

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.

5. Add the exon structure

genestructure projects the human MANE Select transcript, NM_033034.3, 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 alignment with the 7 exons overlaid, one color per exon, projected through the gaps of each row. RING and B-box fall in the large first exon, and the B30.2/SPRY domain spans exons 5 to 7, most of it in exon 7.

6. Add per-codon selection from HyPhy

HyPhy’s FEL (Fixed Effects Likelihood) fits a synonymous rate (alpha) and a non-synonymous rate (beta) at every codon and tests one against the other:

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 is 0.82, mild purifying selection. At the default p<0.1, FEL calls 33 of 600 codons (5.5%) positively selected and 80 (13.3%) negatively selected.

The track draws beta/alpha per codon, clamped at 5, because the ratio nears 10^8 at human residue 337, the strongest positive site (beta = 23.5, p = 0.0063). The script writes the values to trim5-dnds-values.json, one per nucleotide column, for a columnTracks bar (see data layers):

Red bars are FEL’s per-codon dN/dS, clamped at 5, above the exon-colored alignment. The bars cluster over the SPRY half on the right and thin out over the first three exons on the left.

7. Label the SPRY variable patch

Sawyer et al. (2005) found most of the gene’s positively selected sites in a small SPRY patch that overlaps variable region 1 (V1), human residues 326-341, where the HIV-1 restriction specificity maps (Stremlau et al. 2005; Yap et al. 2005; Ohkura et al. 2006). Two highlights label the patch and R332 at its center:

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

start and end count nucleotides along the human row: residue 326 starts at 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 bars of the dN/dS track. FEL calls 3 of the patch’s 16 codons (18.8%) positively selected and none negatively selected, twice the positive fraction of the whole B30.2/SPRY domain (19 of 213 codons, 8.9%). The 16-codon windows starting at human residues 376-382 hold 4 positive sites each, one more than the patch.
R332 boxed at the center of the V1 patch, under a dN/dS bar shorter than most around it. At R332, alpha = 1.77 and beta = 3.80 give a ratio of about 2.1, but p = 0.52 misses FEL’s threshold in this 32-species set, so R332 is not one of the 3.

8. Check the RING and B-box control

RING and B-box are the zinc fingers that fold the protein and drive its E3 ubiquitin ligase activity, so the same analysis should call them conserved:

RING (human aa 15-59, 45 codons) and B-box (aa 90-132, 43 codons) under the dN/dS track. 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.

9. Open the result

The snapshot combines every layer above:

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.

ggmsaview reads the same hosted alignment, exon GFF and dN/dS values, and takes the three highlights as a data frame:

library(ggmsaview)
data <- "https://gmod.org/JBrowseMSA/demo/data/trim5/"
dnds <- jsonlite::fromJSON(paste0(data, "trim5-dnds-values.json"))
sites <- data.frame(
  row = "human",
  start = c(43, 976, 994), end = c(396, 1023, 996),
  label = c("RING + B-box (purifying control)", "SPRY V1 patch", "R332"),
  color = c("rgba(21,101,192,0.25)", "rgba(255,140,0,0.3)", "#e3242b")
)
w <- msaview(
  msa = paste0(data, "trim5-cds.stock"),
  gff = paste0(data, "trim5-exons.gff"),
  color_scheme = "nucleotide"
) +
  geom_msa_track(dnds, name = "dN/dS (FEL)", domain_max = 5,
    color = "#c0392b") +
  geom_msa_highlight(sites)
w
ggplot2::ggsave("codon-final-view-r.png",
  msaview_plot(w),
  width = 12, height = 8)
The ggplot2 figure of all 1,800 columns, with the exons in alternating shades and the tree on the left. The dN/dS bars above thin out over the RING+B-box control on the left and crowd over the SPRY V1 patch and R332 on the right.

Reproduce it end to end

build_codon_selection.sh runs the steps above in order:

  1. It downloads every TRIM5 ortholog and keeps one coding sequence per species, the curated NM_ transcript where one exists, and drops any that translates with an internal stop codon.
  2. MACSE aligns the coding sequences codon by codon, and FastTree infers a tree from the protein translation. The script writes both into one Stockholm file.
  3. react-msaview-cli genestructure projects the human exons onto every row.
  4. HyPhy FEL tests every codon for selection, and the script writes the clamped dN/dS ratios as the track’s values and prints the numbers quoted above.

It needs the tools under Prerequisites.

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

The script writes every file into the current directory. For another gene, edit GENE_ID and the --transcript argument to genestructure.

See also

References

  • Sawyer SL, Wu LI, Emerman M, Malik HS. Positive selection of primate TRIM5alpha identifies a critical species-specific retroviral restriction domain. Proceedings of the National Academy of Sciences 102:2832-2837 (2005).
  • Stremlau M, Perron M, Welikala S, Sodroski J. Species-specific variation in the B30.2(SPRY) domain of TRIM5alpha determines the potency of human immunodeficiency virus restriction. Journal of Virology 79:3139-3145 (2005).
  • Yap MW, Nisole S, Stoye JP. A single amino acid change in the SPRY domain of human Trim5alpha leads to HIV-1 restriction. Current Biology 15:73-78 (2005).
  • Ohkura S, Yap MW, Sheldon T, Stoye JP. All three variable regions of the TRIM5alpha B30.2 domain can contribute to the specificity of retrovirus restriction. Journal of Virology 80:8554-8565 (2006).
  • Ranwez V, Douzery EJP, Cambon C, Chantret N, Delsuc F. MACSE v2: toolkit for the alignment of coding sequences accounting for frameshifts and stop codons. Molecular Biology and Evolution 35:2582-2584 (2020).
  • Price MN, Dehal PS, Arkin AP. FastTree 2: approximately maximum-likelihood trees for large alignments. PLOS ONE 5:e9490 (2010).
  • Kosakovsky Pond SL, Frost SDW. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Molecular Biology and Evolution 22:1208-1222 (2005).