Abstract

INTRODUCTION

During the course of an immune response, B cells that initially bind antigen with low affinity through their immunoglobulin (Ig) receptor are modified through cycles of somatic hypermutation and affinity-dependent selection to produce high-affinity memory and plasma cells. This affinity maturation is a critical component of T cell-dependent adaptive immune responses. It helps guard against rapidly mutating pathogens and underlies the basis for many vaccines. Somatic hypermutation is a process unique to B cells responding to antigen that results in a mutation rate that is 7–8 orders of magnitude above normal background, thus introducing about one point mutation per cell per division in the Ig receptor (1,2). Understanding the somatic hypermutation process also has applications far beyond pathogen responses. It has been found to occur in autoimmune responses, and in several proto-oncogenes (3). We have also demonstrated that somatic hypermutation can act genome-wide and thus represents a risk for genomic instability (4).

The ability to detect selection from mutated Ig sequences is a critical part of many studies. Current methods are based on comparing the observed frequency of replacement (i.e. non-synonymous) mutations to their expected frequency under the null hypothesis of no selection (5–8). Elevated levels indicate positive selection, while decreased levels indicate negative selection with significance determined by a binomial test. It is common to look for negative selection in the framework regions (FWRs), which provide the structural backbone of the receptor, and positive selection in the complementary determining regions (CDRs), where most contact residues for antigen binding are found. As the intrinsic biases of somatic hypermutation can give the appearance of selection (9), a significant challenge for these methods is calculating the expected frequency of replacements under the null hypothesis of no selection. We previously developed the Focused binomial test for detecting selection that improved upon existing methods by fully accounting for microsequence specificity and base substitution bias in somatic hypermutation (10). The Focused test also corrects for the decrease in specificity due to cross-talk in other methods by using a carefully derived null model of mutation.

Our web site implementing the Focused test has been increasingly used by a large number of groups since the initial publication of the method in 2008 (10). Here, we present an improved web server that includes many of the most requested features from users, such as V(D)J germline segment identification with IMGT alignment and batch submission of sequences. For comparison, we also integrate results from the previously proposed Local binomial (5) and multinomial (7) tests. Note that we have shown that the multinomial test is mathematically equivalent to a much simpler binomial test, referred to in this article as the Global binomial test. Along with these features, we have implemented an improved Z-score-based statistic that increases sensitivity and allows for the combined analysis of multiple-independent Ig sequences or clones.

METHODS

Analyzing sequences for detecting selection

Calculating the observed number of mutations

By comparing the input sequence to its associated germline, our program identifies the mutations and determines the number of replacement (R) and silent (S) mutations. Each mutation is considered independently in its germline context when determining whether it is an R or S. R and S mutations that fall into CDR and FWR are tabulated individually using the boundaries indicated by the user. A checkbox is provided to indicate that sequences are clonally related. In this case, each group of sequences associated with a germline is analyzed as a single unit and only unique mutations are used in the analysis. That is, the same base substitution occurring at the same position in multiple sequences is counted only once.

Estimating the expected frequencies of mutations

Having computed the observed number of mutations, the next step is to compute the expected number of mutations under the null hypothesis of no selection. Expectations are computed independently for each germline sequence in the input as previously described (10). A significant advantage of the Focused test over previous methods is that our null model fully accounts for the effects of microsequence specificity (16) and also introduces the well-characterized transition bias of somatic hypermutation (17,18). Briefly, the expected number of R mutations in the CDR () is the sum of the product of an individual point mutation falling in the CDR and the probability the base substitution results in an R mutation [Equation (1)]:

(1)

where i is summed over all positions in the region (i.e. CDR or FWR) and b over all possible nucleotides ({A,C,T,G}). In this equation, is a vector containing the nucleic content of each position in the germline sequence, is the mutability index for position i in germline [as explained in (10)], M_a→b is the relative rate in which nucleotide a mutates to b (while M_a→a=0) and is an indicator function that is 1 in cases where a mutation in position i from a to b results in a replacement mutation and 0 otherwise.

The binomial framework for detecting selection

The Local binomial, multinomial and Focused binomial tests for selection all determine whether the observed number of R mutations (x), in either the CDR or FWR, is significantly different than the expected number () out of n observed mutations. The expected frequency (x/n) under the null hypothesis of no selection is denoted by p. For x/n≤p (an indication of negative selection), the significance of the test is calculated as the probability of observing x or fewer R mutations by adding half the probability density function (P) at x to the cumulative distribution function (F) at (x−1). The P-value of x/n>p (an indication of positive selection) is one minus that of negative selection.

(2)

where , and Θ(x) is 1 for x>0 and 0 elsewhere. Table 1 defines these parameters for each implemented test.

Table 1.

Definition of x, n and p for each of the implemented tests to detect selection in the CDR

Test	x	n	p
Focused	RCDR	RCDR+SCDR+SFWR
Local	RCDR	RCDR+SCDR
Global	RCDR	RCDR+SCDR+RFWR+SFWR

Open in a separate window

Equivalent tests for detecting selection in the FWR are obtained by swapping CDR and FWR in the definitions.

Using Z-scores to compute P-values

The standard way of calculating P-value (P_Binom) described in the previous section is conservative in the sense that the resulting specificity is often greater than the cutoff used to obtain it. This limits sensitivity and is particularly a problem when the total number of mutations per sequence is relatively low (~10), as is common for many Ig data sets, but is still significant for dozens of mutations. To correct this problem, we have implemented a Z-score-based method for computing the P-value. The Z score is defined as follows: for a random variable x_i corresponding to the x for sequence or clone i as defined in Table 1, the associated z_i score (which is itself a random variable) is defined as:

(3)

whereμ_i=p_in_i and are the mean and standard deviation of x_i, respectively. In order to obtain a P-value from the z_i, we use the normal approximation and get:

(4)

where erf is the Gauss error function; the reader is referred to Figure 1 to see graphically the different P-values definitions [normal approximation and the traditional way of defining it for a binomial distribution, Equation (2)]. We call this the Focused Z test (P_z) to distinguish it from the Focused binomial test (P_Binom).

Open in a separate window

Figure 1.

Graphic representation of the calculation for P_Binom, the binomial-based test (projected upwards) and, P_Z, the Z-score based test (projected downwards). The binomial distribution depicted by the bars corresponds to x=11, n=20 and p=0.389 (typical parameters from the simulations presented in Figures 2 and ​and3).3). P_Binom [Equation (2)] is the sum of the yellow bars, while P_Z is calculated by approximating the binomial distribution with a normal distribution with the same mean (dotted line) and variance (1 SD interval around the mean is shaded in darker gray) and taking the integral (area under the curve) from x=11 to ∞ (shaded red).

Detecting selection in groups of sequences

Our web server also implements a test for detecting selection in a group of independent sequences. This can be helpful to improve sensitivity if the sequences do not exhibit statistically significant selection when analyzed individually. Previous tools did not allow for such analysis since grouping P-values computed by Equation (2) from sequences with different p's and n's cannot be done for a fixed P-value cutoff in a simple way. However, this can be done using the Z-score approach, through the application of Stouffer's Z-score method (19), which has proven to be superior to alternatives (20). The resulting P value tests whether the collection of observed mutations came from binomial distributions with associated probabilities p_i. To analyze G-independent sequences, Z is defined as the mean of the previously computed Z-scores z_i, i1…G from Equation (3). Since the x_i's in this equation are independent, the expected mean and variance of Z are 0 and , respectively. Using the normal approximation once again, the resulting P value for detecting selection is:

(5)

By default, this test is applied to compute P values for each group of sequences associated with a single germline, as well as for the set of all sequences provided in the input. In order to compare with the original Focused binomial test, a checkbox allows the use of P_Binom for individual sequences. If the user indicates that sequences are clonally related then the set of clones is analyzed as a group (in which case, G will be the number of germlines provided in the input).

PERFORMANCE RESULTS

We have previously shown, using a simulation-based validation strategy, that the Focused binomial test provides the best trade-off between sensitivity and specificity compared with other available methods (10). Furthermore, we found that the Focused binomial test is able to detect selection in experimentally derived Ig sequences that have undergone affinity maturation, while maintaining good specificity on non-functional Ig sequences where selection is not a force (10). An independent study has corroborated these findings (21). As described below, we use a simulation-based validation strategy to evaluate the performance characteristics of our updated Z-score-based methods.

Simulation of somatic hypermutation and selection

Simulation enables us to produce synthetic sequence data with a prescribed number of somatic mutations subject to varying amounts of positive and/or negative selection pressures. The simulation is initiated with a single IMGT formatted Ig V germline sequence. Mutations are introduced one-by-one along the entire length of the sequence (excluding gaps) in two steps. First, the position is chosen stochastically based on the microsequence specificity of each nucleotide. Second, the particular substitution is determined accounting for transition bias. Selection (ξ) is implemented separately in the CDR and FWR as uniformly increasing or decreasing by a log factor the expected probability of the expected frequency of R mutations in each region:

(6)

where p_region is the expected frequency of R mutations defined by the Local test (Table 1) and π_region is the actual probability applied to the simulated sequences in the region of interest (i.e. CDR or FWR). For example, ξ_CDR values of −1, 0 and 1 yield synthetic data with negative, neutral and positive selection in the CDR, respectively.

Using Z-scores improves sensitivity

To evaluate performance, a synthetic data set of mutated sequences was generated. This data set included 10000 sequences with strong positive selection in the CDR (ξ_CDR=1), and 10000 sequences with neutral selection in the CDR (ξ_CDR=0), allowing us to evaluate sensitivity and specificity, respectively. All sequences included strong negative selection in the FWR (ξ_FWR=−1). Each simulation starts with the IgHV7–3 germline sequence and introduces between 5 and 35 mutations per sequence, stochastically determined to reflect the numbers seen in experimental data. Results starting with other germline segments are similar. Plotting the fraction of sequences predicted as being positively selected for varying alpha cut-off values in the first data set against the second produces a receiver operating characteristic (ROC) curve. The ROC curves (Figure 2) for both the Focused binomial test (yellow) and the Focused Z test (red) are comparable, confirming the validity of the Z-score-based method. However, the position on the ROC curve where the α cutoff is 0.05 for the Focused binomial test (yellow X) falls short of the expected 1−specificity of 0.05, indicated by the dashed line. Applying the Z-score-based method corrects for this discrepancy leading to an improvement in sensitivity of ~10%, as shown by the red X marking the same position on the ROC curve for the Focused Z test. The accuracy of the Focused Z test at α=0.05 is 0.72 while it is 0.68 for the Focused binomial test. Along with the Focused test, a checkbox on our web site also allows users to include results from the Local and Global (also known as the multinomial) tests (Table 1). However, it is important to note that these tests, as originally proposed, did not fully include for microsequence specificity or substitution biases when computing the expected frequencies of mutations, whereas our implementation accounts for these intrinsic biases. In addition, P-values for all the approaches are calculated using the Z-score-based method. Only results from the Focused test are output by default since this method exhibits better sensitivity than the Local test. Furthermore, we strongly caution against using the Global test because it fails to maintain the defined specificity [Figure 3 and (10)], but nevertheless make its output available for comparative studies.

Open in a separate window

Figure 2.

ROC curves comparing P_Binom and P_Z on simulated data. Sensitivity is based on detecting selection in simulated sequences with ξ_CDR=1 and ξ_FWR=−1. Specificity is based on detecting selection in simulated sequences with ξ_CDR=0 and ξ_FWR=−1. The dotted line indicates the expected specificity at α=0.05. The position of this cutoff on the P_Binom and P_Z curves is indicated by yellow and red Xs, respectively. The inset shows the ROC curves from P_Z,G calculated for independent sequences grouped into sizes (G) of: 1,2,4 and 8.

Open in a separate window

Figure 3.

P_Z ROC curves comparing Focused (red triangle), Local (blue square) and Global (green circle) tests for detecting selection. The respective colored X's indicates the positions on the ROC curves that corresponds to the α=0.05 cutoff. The sensitivity and specificity are computed on the same simulated data set from Figure 2.

DISCUSSION

The ability to analyze mutation patterns in Ig sequences and detect antigen-driven selection is critical to understanding adaptive immunity. We have presented a web site that makes available our previously published Focused test, along with the Local and Global tests based on statistics proposed by other groups (5–7). Consistent with previous results from our own group and others (10,21), a simulation-based validation found that the Focused test exhibited the best performance. The web site offers an integrated pipeline where users can carry out V(D)J identification with IMGT alignment using SoDA (13), quantify the mutational load in their sequences and analyze the mutations for evidence of positive and/or negative selection (Figure 4). The ability to carry out batch processing and analyze related sequences as a single group will be critical to gain insights from large-scale data sets generated by emerging techniques such as expression cloning (23) and deep sequencing of B cell repertoires (24).

Open in a separate window

Figure 4.

Example output from the website. (A) When sequences are independent, P_Z is output for each individual sequence and P_Z,G is output for the set of sequences associated with each germline. (B) If the sequences are clonally related, then P_Z is output for each clone. Results from each germline (A) or clone (B) are further combined to produce a final P_Z,G for the entire data set (labeled as and , respectively). The resulting P values for the different tests in each region (CDR and FWR) are color coded: red for positive selection and green for negative selection, with light and dark colors indicating non-significant (|P|>0.05) and significant (|P|≤0.05) selection, respectively. Note that P values less than zero are indicative of negative selection. P_Z and P_Z,G are indicated on the side of the plot for selected cases.

FUNDING

National Institutes of Health (R03AI092379-01) (in part). Funding for open access charge: National Institutes of Health (R03AI092379-01).

Conflict of interest statement. None declared.

ACKNOWLEDGEMENTS

REFERENCES