Families, evolution and variants

A protein is understood through its relatives. GL1X gathers a chain's homologs from the archive and aligns them, measures how conserved each position is and paints it on the structure, builds trees, finds coevolving residue pairs that predict contacts, clusters sequences, searches the archive with profile HMMs, annotates sequences, carries annotations over from structural neighbours, and places variants on structures with the features models of mutational effects need — with GL1F models on chain predicting them. Every method follows its reference implementation and is checked against it (validation).

Alignments · Conservation · Trees · Couplings · Clustering · Profile HMMs · Sequence annotation · Annotation transfer · Variants · API and MCP · Validation

Multiple sequence alignment #

Homologs. Given a record's chain, its whole sequence is searched with BLAST against every protein of the archive (E ≤ 10⁻³); hits covering at least half of it are kept, the query's part of each, dropping any sequence more than 95% identical to one already kept, up to 150 (all adjustable). Records and FASTA sequences can be added or aligned alone.

Alignment follows MUSCLE's progressive strategy (R C Edgar. Nucleic Acids Res 32, 1792–1797, 2004):

1. distances: k-mer similarity for more than 30 sequences, else identity of pairwise alignments (Kimura's correction);
2. a UPGMA guide tree;
3. sequence weights (CLUSTAL W's, from the tree);
4. profile–profile alignment up the tree, columns scored by weighted sum of pairs with BLOSUM62 (+1.5 offset), gap opening 8 scaled by each column's occupancy and by the share of gaps already opening there, extension 1, terminal gaps at half cost;
5. the tree estimated again from the draft alignment and the sequences realigned along it;
6. refinement: the alignment split at each edge of the tree, the two profiles realigned, the change kept when the sum-of-pairs score improves (two passes; one above 150 sequences, none above 400).

Up to 500 sequences of up to 3,000 residues; an alignment stops when its time limit (150 s) is reached or its requester leaves. Alignments made elsewhere — up to 5,000 sequences of 10,000 columns — are read as aligned FASTA, A2M/A3M, Clustal or Stockholm and scored the same way; a tree is drawn for 1,000 sequences at most (and not when the pairs of sequences times the columns exceed 2.5 × 10⁹).

Conservation and profiles #

Jensen–Shannon divergence (J A Capra, M Singh. Bioinformatics 23, 1875–1882, 2007), exactly as their score_conservation.py: each column's amino-acid distribution (Henikoff sequence weights, pseudocount 10⁻⁷) against the BLOSUM62 background distribution, multiplied by the weighted share of non-gap residues; columns with more than 30% gaps are not scored; the score is smoothed over a window of 3 columns on each side with weight λ = 0.5. 0 is random, 1 is invariant.

The reference chain's residues sit on the columns, so conservation is painted on its structure (white to deep purple) and conserved surface patches — functional sites — stand out.

PSSM. Position-specific log-odds (bits) with PSI-BLAST's BLOSUM62 conditional pseudocounts (β = 10, α = the number of distinct residues in the column − 1, at least 1: PSI-BLAST's Nc − 1 taken column by column), the substitution scores of the variant tools.

Sequence weights, the identity of every pair and of each sequence to the reference come with each alignment.

Trees #

Neighbour joining (N Saitou, M Nei 1987; Studier–Keppler's formulation) on Kimura distances for proteins (−ln(1 − p − 0.2p²)), or p-distances, or Poisson-corrected; written as Newick, drawn in the workbench.

Direct coupling analysis #

Residues in contact coevolve. Mean-field DCA (F Morcos et al. PNAS 108, E1293–E1301, 2011), as the authors' dca.m: 21 states (gap and 20 amino acids), sequences reweighted at 80% identity, pseudocount 0.5, the coupling matrix from the inverse of the correlation matrix (Cholesky), direct information for each pair of columns, and the APC-corrected Frobenius norm of the couplings in the zero-sum gauge without gaps (S Ekeberg et al. Phys Rev E 87, 012707, 2013), the better contact predictor. Columns with more than half gaps are left out; up to 300 columns. Coupled pairs come with their distance in the reference structure, and the precision of the top L/5, L/2 and L pairs as contacts (Cβ within 8 Å, six or more residues apart). A family with fewer effective sequences than half its length gives a warning: its couplings are noise.

DCA runs as a job (POST /api/bio/msa/dca, then GET /api/bio/msa/dca/{id}): one runs at a time and four may be queued or running (two of one client), an alignment given is read when the job is asked for (one that cannot run is refused at once), results kept for an hour.

Clustering and non-redundant sets #

CD-HIT's greedy incremental clustering (W Li, A Godzik. Bioinformatics 22, 1658–1659, 2006): sequences longest first, each joining the first representative it matches at the identity threshold (identical residues of a global alignment over the shorter sequence), else founding a cluster; a short-word filter (5-mers at 70% identity and above, 4-mers at 60%, 3 at 50%, 2 below) skips alignments that cannot reach the threshold — at 90% identity most of them, below about 80% few, so that every sequence is then aligned with every representative. Up to 5,000 sequences of 5,000 residues, within two minutes. Datasets (POST /api/bio/dataset) take non_redundant to keep one representative per cluster — the standard guard against overfitting.

Profile HMM search #

A profile hidden Markov model (a Pfam family, or one built with hmmbuild) is searched against every protein of the archive with hmmsearch's pipeline (S R Eddy. PLoS Comput Biol 7, e1002195, 2011), in GL1X's own code:

• the model configured as hmmsearch configures it: local, multihit, the length model set for each target, HMMER's background and null model, degenerate residues scored by their average;
• filters: MSV (ungapped segments, P ≤ 0.02, in HMMER's integer arithmetic), the bias filter (a two-state composition model), Viterbi (P ≤ 10⁻³) and Forward (P ≤ 10⁻⁵, raised to E/Z so a reportable hit is never dropped); P-values from the model's own calibration (STATS LOCAL lines);
• domain definition by posterior decoding: regions where residues are likely in the model, split into domains by clustering 200 sampled paths when more than one is expected, each domain rescored alone;
• the null2 composition correction: sequence and domain scores lose the part a biased composition explains — without it, low-complexity proteins score tens of bits too high;
• E-values = P × the number of sequences searched; a domain's conditional E-value (c-Evalue) = P × the number of sequences reported and its independent E-value (i-Evalue) = P × the number searched, as hmmsearch gives them; a sequence is included at E ≤ 0.01, a domain of it at c-Evalue ≤ 0.01. A model without the MSV calibration line is not filtered by MSV.

Each hit lists its domains with envelope and alignment coordinates on the sequence and the model, scores, biases and the alignment. Searches of 2,000 sequences or fewer (a collection) score every sequence fully, like hmmsearch --max.

Sequence annotation #

For a record's chain the structure's view comes beside the prediction: B-factors or pLDDT, the residues it lacks, DSSP and accessibility per residue — a predicted disordered stretch missing from a crystal structure, or a low-pLDDT loop of a model, confirms each other.

Annotation transfer #

What a chain's structural neighbours (structure search, TM-score ≥ 0.5 by the query, one per record) say about it, each at the confidence their similarity supports:

• EC numbers, truncated to the digits the neighbour's identity in the superposition carries (W Tian, J Skolnick. J Mol Biol 333, 863–882, 2003): four at ≥ 60%, three at ≥ 40%, two at ≥ 30% with TM-score ≥ 0.7, the class alone from a TM-score of 0.5; each with how many neighbours support it;
• ligand sites: the neighbours' ligands (not additives or ions) with their binding residues mapped onto the query through the structural alignment — where a cofactor or a drug would sit in a structure that has none;
• pocket sites: the neighbours' pockets carried onto the query the same way, when most of a pocket's lining on the aligned chain aligns (three residues at least); pockets of several neighbours sharing 30% of their carried residues (Jaccard) are one site, ranked by the TM-scores of the neighbours that support it, with the best and mean druggability, how many held a ligand and which, and each residue with how many pockets line it — a binding site an apo structure or a model does not show by itself;
• names the neighbours give their molecule.

Gene Ontology terms are not in structure files, so GL1X does not invent them; EC numbers, ligands and names are what the archive holds.

Variants and mutational effects #

Variants are written as A123V, p.Ala123Val, 123V or with a chain A:K72E, numbered along the chain's whole sequence or by the author; the variants of one request are of one chain. Each is placed on the structure and described by 61 features (mut_), the inputs of models of mutational effects:

With measured effects (labels, e.g. ΔΔG), the variants become a training dataset (CSV) for a GL1F model. With model, a GL1F model on chain — a free or tips forest model with public internals, a regression (ΔΔG) or a binary classification (damaging), whose inputs are mut_ features — predicts each variant; GL1X replays the model reconstructed from the chain and checks its first prediction against the chain's own answer.

A scan scores every substitution of a stretch (up to 1,500 positions): the family log-odds change when an alignment or homologs are given, else BLOSUM62, drawn as a residues × amino acids heatmap; a model's predictions replace them when given.

In the explorer #

The menu's Bioinformatics onchain section opens each tool on its own page, and the research workbench runs them on one structure:

• Families & evolution (/bio/families): Alignment & conservation aligns a record's homologs, sequences or records — or scores an alignment made elsewhere — and shows the alignment with its conservation bar, the structure painted in ConSurf's nine grades (a column click shows its residue), the neighbour-joining tree, and downloads (aligned FASTA, Newick, conservation table); Couplings (DCA) runs a job with its progress and draws the strongest couplings on a contact map against the structure's contacts; Clustering, HMM search (paste or load a HMMER3 model; hits with their domain alignments) and Sequence annotation (tracks of disorder, coiled coils, hydropathy and low complexity) complete it.
• In the workbench (/bio/workbench/{contract}/{token}): Conservation paints the family's conservation and draws the couplings as rods between residues (green when the structure has them in contact); Variants · ΔΔG scans every substitution of a stretch as a heatmap — family log-odds, BLOSUM62, or a GL1F model's ΔΔG or damage probability — with the mean per residue painted on the structure, and places named variants; Annotation paints disorder, coiled coils, transmembrane helices or B-factors and pLDDT; Neighbours lists the EC numbers, ligand sites and pocket sites carried over and paints the chosen site.

API and MCP #

Validation #

• Alignment, on 8 curated reference alignments (Pfam seeds of 12TM_1, 3Beta_HSD, ArsP_1, COX2_TM, Cyclin_N with 95 sequences, Phage_Coat_Gp8, Cys_knot, and HMMER's 45 globins), mean sum-of-pairs / total-column scores: GL1X 0.900 / 0.754; MUSCLE 5.1 0.901 / 0.756; MAFFT 7.505 0.889 / 0.733; Clustal Omega 1.2.4 0.868 / 0.680. On simulated families (40–200 sequences, 30–60% identity): GL1X 0.839, MAFFT 0.880, Clustal Omega 0.793. 500 sequences of ~300 residues align in 2.8 s.
• Conservation follows score_conservation.py term by term (Henikoff weights, pseudocount, background, gap penalty, window), checked against a computation by hand; DCA's direct information equals dca.m (ProDy's reference output) to 10⁻⁶; a family of 1,000 sequences × 200 columns takes 11.7 s on two CPUs.
• Trees: splits and branch lengths identical to Biopython 1.88's neighbour joining on 15 matrices of 3–25 taxa.
• Clustering against CD-HIT 4.8.1 on four simulated families at 90%, 70% and 50% identity: cluster counts equal in 7 of 12 cases and within 10% in all (450 vs 451, 131 vs 133, 38 vs 42).
• Profile HMM search against hmmsearch 3.4 with default settings on 22,338 sequences (2,005 NCBI proteins, 215 globins and cyclins, 79 Pfam family members, 20,050 shuffled decoys) and 8 models of 52–449 match states: the same hits at E ≤ 0.01 for every model (none missed, none extra); sequence scores within 0.05 bits (hmmsearch prints tenths) for single-domain regions and within 0.33 bits where paths are sampled; composition corrections within 0.38 bits; MSV filter decisions identical. The whole pipeline runs at about 1.5 ns per residue and model state on one core.
• Coiled coils equal ncoils (26,128 residues, within 0.0006).