Using sequence similarity networks for visualization of relationships across diverse protein superfamilies

Abstract

The dramatic increase in heterogeneous types of biological data--in particular, the abundance of new protein sequences--requires fast and user-friendly methods for organizing this information in a way that enables functional inference. The most widely used strategy to link sequence or structure to function, homology-based function prediction, relies on the fundamental assumption that sequence or structural similarity implies functional similarity. New tools that extend this approach are still urgently needed to associate sequence data with biological information in ways that accommodate the real complexity of the problem, while being accessible to experimental as well as computational biologists. To address this, we have examined the application of sequence similarity networks for visualizing functional trends across protein superfamilies from the context of sequence similarity. Using three large groups of homologous proteins of varying types of structural and functional diversity--GPCRs and kinases from humans, and the crotonase superfamily of enzymes--we show that overlaying networks with orthogonal information is a powerful approach for observing functional themes and revealing outliers. In comparison to other primary methods, networks provide both a good representation of group-wise sequence similarity relationships and a strong visual and quantitative correlation with phylogenetic trees, while enabling analysis and visualization of much larger sets of sequences than trees or multiple sequence alignments can easily accommodate. We also define important limitations and caveats in the application of these networks. As a broadly accessible and effective tool for the exploration of protein superfamilies, sequence similarity networks show great potential for generating testable hypotheses about protein structure-function relationships.

Figure 1. Sequence similarity network topology changes in a predictable way with the stringency of the threshold.

A. Thresholded sequence similarity networks represent sequences as nodes (circles) and all pairwise sequence relationships (alignments) better than a threshold as edges (lines). The same network, depicting three simulated protein classes, is shown here at four different thresholds. At stringent thresholds, the sequences break up into disconnected groups; within each group the sequences are highly similar. The relative positioning of disconnected groups has no meaning, while the lengths of connecting edges tend to correlate with the relative dissimilarities of each pair of sequences. As the threshold is relaxed and edges associated with less significant relationships are added to the network, groups merge together and eventually become completely interconnected. B. Simulated dendrogram for a sequence set that might give rise to the network in A.

Figure 2. Comparison of trees and networks: amine-binding GPCRs.

A. Neighbor-Joining tree describing the interrelationships of 42 amine-binding human GPCR domains. Sequences are labeled according to the common name for their class (e.g., the sequence labeled α1D is adrenoceptor α1D; see additional data file 5 for all sequence database identifiers). B. Sequence similarity network including the same 42 sequences as in (A). This network was thresholded at a BLAST E-value of 1×10⁻³³: only edges associated with E-values more significant than 1×10⁻³³ are included in the network. This network contains 324 edges; the worst edges displayed correspond to a median of 30% identity over an alignment length of 280 amino acids. See Table I for a quantitative comparison of the two representations. The sequences labeled (a) and (b) are discussed in the text.

Figure 3. Sequence similarity networks are useful tools for exploration of the kinase superfamily.

Two ways of coloring the same network of 513 human kinase domains are shown. The network is thresholded at a BLAST E-value of 1×10⁻²⁵. The worst edges displayed correspond to a median of 29% identity over alignments of 260 residues. A. Network colored by kinase class. B. Network colored by the presence of a catalytic Lys in the “VAIK” motif: Each of the 513 sequences was aligned to a sequence model of the kinase domain, and the identity of the residue at the catalytic Lys position is mapped to the network. *Note that MAP2K1 and MAP2K2 registered a Lys to Arg substitution due to a sequence alignment error. The other labeled kinases truly do not contain a homologous catalytic K, but only the WNK kinases have been shown to have kinase activity. See Table II for statistics.

Figure 4. Visualizing very distant interrelationships between GPCRs.

A. 605 human Class A: Rhodopsin-like GPCR domains. This sequence set includes the 42 amine-binding sequences from Table II and Fig. 2. This network was thresholded at a BLAST E-value of 1×10⁻¹¹; the worst edges displayed correspond to a median of 24% identity over an alignment length of 210 amino acids. Sequences colored red for “Class A: Rhodopsin-like” were not classified to a specific subgroup within the class. B. 766 human GPCR domains. This sequence set includes all of the 605 Class A sequences from (A), now colored dark grey. This network was thresholded at an E-value of 1×10⁻², and the worst edges displayed correspond to a median of 22% identity over an alignment length of 120 amino acids. Also included in both networks is a set of 17 sequences in light grey. These sequences were used here as negative controls, and were randomly drawn from the human proteome and not annotated as GPCRs. The red and green clan labels correspond to PFAM clans. The sequences that are associated with or that are extremely similar to high resolution structures are noted [PDB identifiers 1F88, 2VT4, 3EML, 2RH1, and 2R4R[25]]. See Table S1 for network statistics.

Figure 5. Crotonase superfamily: sequence similarity network from full-length sequences and from just the domain in common.

The displayed networks all describe the pairwise relationships between 1,170 sequences from the crotonase superfamily. A. Network colored by family annotation, involving full-length sequences, thresholded at an E-value of 1×10⁻³⁰. The worst edges displayed correspond to a median of 33% identity over alignments of 250 residues. B. The full-length network from A with nodes colored by sequence length and edges colored by alignment length. The same bifunctional enoyl-CoA hydratases (bECH) are marked with a dashed oval in B and C. C. Network colored by family annotation, involving just the crotonase domain, thresholded at 1×10⁻²⁹. The worst edges displayed correspond to a median of 38% identity over alignments of 180 residues. D. 17 selected edges from the network in A and B. In the left panel, for each pair of sequences participating in an alignment, the log E-value versus the HMM used to define the crotonase domain is shown for each sequence—the single domain ECH (sECH) is on the bottom, and the second member of the pair is on the top—and the log BLAST E-value for the alignment between the two is in the middle. Two example bECH and sECH sequences (not alignments) are shown at the bottom of the left and middle panels. In the middle panel, each amino acid in each sequence from the 17 alignments is colored according to whether it was aligned to the crotonase domain defined by the HMM, and/or was paired to the other sequence in the BLAST alignment used to define an edge. Locations of six of these edges are marked in the enlarged view of the network in A in the right panel. The locations of the example bECH and sECH sequences are marked in the right panel using stars. See Tables S1 and S3 for quantitative comparisons.

Figure 6. Domain shuffling in the enoyl-CoA hydratase family.

The displayed networks contain all 410 enoyl-CoA hydratases from the crotonase superfamily network in Fig. 5A. The network is thresholded at a BLAST E-value of 1×10⁻⁵⁰; the worst edges displayed correspond to a median of 40% identity over alignments of 260 amino acids. A. Network nodes colored by sequence length and edges colored by alignment length. B. Network nodes colored by species kingdom (Fungi, Metazoa, Viridplantae) or superkingdom (Bacteria, Eukaryota, Archaea). The same archaebacterial bifunctional enzymes are marked with a dashed oval in both A and B. C. Representative domain structures for the three major classes of enoyl-CoA hydratase-containing sequences, with domains defined using PFAM HMMs.