Identifying HIV sequences that escape antibody neutralization using random forests and collaborative targeted learning

2.1 Counterfactual antibody resistance probability

We assume that we have access to a database like the aforementioned CATNAP. An example illustration of one such data set is given in Table 1. The data consist of n HIV Env sequences. For each sequence, we have a measure of sensitivity to an antibody of interest (e.g., 1 = can effectively be neutralized by the antibody, 0 = cannot) and basic information about the origin of the sequence. The remaining columns represent the AAs that make up the Env proteins of the sequence. Each AA position, or residue, of the Env protein can assume one of the 22 different values – one of the 20 AAs that build proteins, a frameshift mutation (indicating an insertion or deletion of nucleotide bases), or a gap (an artifact of the sequencing technology to maintain alignment with a referent HIV sequence). For our purposes, it suffices to think of each AA residue as a categorical variable that theoretically could assume up to 22 different values, while in practice we generally observe between two and four unique AAs at each residue. We treat the n sequences as independent, which may be reasonable because (i) almost all sequences in CATNAP are isolated from different individuals and (ii) HIV replicates extremely rapidly, such that any two isolated sequences, even if close in geographic proximity, are likely to be distant ancestors of one another.

We denote the sensitivity outcome by Y , where Y = 1 indicates that the virus is effectively neutralized by the antibody, and Y = 0 indicates that the virus is not effectively neutralized by the antibody. We denote the origin of the virus together with the vector of AA information on the Env sequence as W . For simplification, basic information about the origin is encoded as W 0 . Below, we consider methods for analyzing the causal impact of mutations at each of the J AA residues in turn. It is therefore convenient to introduce the notation W − j to denote the origins and the vector of all Env characteristics with the j th residue removed, where j = 1 , … , J . Thus, the data that we use for the analysis of the j th residue can be considered n copies of the triplets O j = ( W j , W − j , Y ) .

For each AA residue, we define a counterfactual resistance outcome Y ( W j = w ) , which is the binary resistance indicator that would be observed under a hypothetical intervention that fixes the AA at residue j to w ∈ W j , the set of possible AAs at residue j that could be observed in the entire population of virus Env sequences. We also define p j ( w ) = P [ Y ( W j = w ) = 1 ] as the proportion of virus Env sequences in this counterfactual world that would be sensitive to neutralization by the bnAb of interest.

To achieve the goals outlined in Figure 2, we may be interested in estimation and inference about p j ( w ) . For example, to identify gaps in one bnAb, we would be interested in identifying ( j , w ) combinations for which p j ( w ) is low. To complement such a bnAb in a multiple bnAb cocktail, we would look for other antibodies such that the same counterfactual sensitivity probability is high. On the other hand, if the scientific goal of the analysis of the observational neutralization data is to reduce the search space for a sieve analysis, we may instead wish to test the null hypothesis that an AA substitution at residue j has no impact on average bnAb resistance, that is, H 0 : p j ( w ) is constant in w . AA residues for which this null hypothesis is rejected may be advanced for further consideration in a sieve analysis using data from a randomized trial. Regardless of the ultimate scientific goal, the problem of identification and estimation of counterfactual sensitivity probabilities is an important problem.

2.2 Causal identification

We require two main conditions for identifiability of p j ( w ) : (i) conditional exchangeability and (ii) positivity. If these conditions hold, then p j ( w ) is identified by the G-functional p j ( w ) = E [ E ( Y ∣ W j = w , W − j ) ] . We discuss these conditions in detail below and conclude with a discussion of their plausibility in the present context.

2.2.1 Exchangeability

A DAG is displayed in Figure 3. On the left, we display a DAG that is useful for discussing the biology of antibody neutralization. On the right, we display a DAG that is useful for describing backdoor pathways and types of unmeasured confounders that may prove problematic for our proposal.

Figure 3

DAG for HIV resistance. The DAG on the left shows a pipeline of gene expression, which defines how gene regulates downstream activities and further affects the antibody sensitivity through Envelope AAs of interest.

Beginning with the DAG on the left, the nodes labeled G represent genes on the viral genome. Certain genes on the genome correspond with AA residues on the Env protein (nodes labeled AA). These AA configurations result in certain physiochemical properties of the Env protein. The DAG lists several Env characteristic that are putatively important for antibody binding and neutralization. These characteristics act as mediators of AA configuration on antibody sensitivity. The structure (or shape) of the Env protein is determined by the AAs present; this in turn determines which portions of the protein are “exposed” to antibodies. There are other physiochemical properties that are also impacted by the Env protein AA sequence that could impact affinity of antibodies for binding the protein. One particular process that may be important for HIV infectivity and antibody sensitivity is the presence of glycosylated regions of the Env protein. Such regions can play a key role in binding, both of the virus to the host CD4 T-cells, as well as antibodies to the virus. Therefore, these regions are expected to play an especially important causal role in neutralization of viruses via antibodies.

A simplified DAG is shown on the right with nodes collapsed by the role they play in opening/closing backdoor paths between the outcome and AA residues in the protein. As mentioned above, our approach considers each AA in turn. So for a particular residue j , W j is acting as the “exposure” of interest, while all other AA residues are acting as potential “confounders” W k . The genes are unmeasured and play the role of potential unmeasured confounders, while the origin of the virus can be considered another potential confounder W 0 . The mediators in this problem, here collapsed to a single node M , are the aforementioned structural or physiochemical properties of the Env protein. The DAG on the left suffices to discuss the biological plausibility of causal identification in this problem. In particular, there are three types of pathways arising to be relevant to identification. The first are pathways of the form W j ← G 2 → W k → M → Y . In such pathways, there is a gene responsible for coding both the AA residue of interest and other residues that are associated with mediators of antibody sensitivity. In this case, controlling for W k suffices to block the pathway. The second type of pathway is one of the form W j ← G 2 ← W 0 → G 1 → M → Y . Here, the geographic region ( W 0 ) is associated with two genes that code separate residues on the Env protein. This pathway can be blocked by conditioning on W 0 . Finally, there are pathways of the form W j ← G 3 → M → Y . Here, there is a gene that codes the AA residue of interest, but is also associated with some other mediating pathways of antibody sensitivity. For example, it is possible that having more Env proteins expressed reduces antibody sensitivity of viruses. If a particular gene was associated with the number of Env proteins that are expressed and also associated with the AA residue of interest, then there would be no way to block this pathway, since the gene information is unavailable. Thus, in order to formally establish causal identification, we would need to assume that no such genes exist. If not, then the DAG implies exchangeability, that is Y ( W j = w ) ⊥ W j ∣ W − j .

2.3 Motivation for novel method

There are many methods available that could in principle be applied in a non-causal context in this setting. For example, we could first conduct an analysis using LASSO on half of samples. The best value of the tuning parameter λ can be determined via cross-validation to arrive at a final model. The AA residues included in the final model can be included in a logistic regression on the withheld half of samples for formal inference. Another approach would be to use random forest “variable importance” measures, for which we may obtain p -values using a fast variable importance test [16]. However, we argue that methods that have a causal interpretation, even if under limited circumstances, are appealing in this setting since causal inference is fundamentally at the core of the scientific question. Thus, we are motivated to explore causal inference-inspired methodology for estimating the G-functional E [ E ( Y ∣ W j = w , W − j ) ] and associated significance tests for AA residues.

There are two related challenges associated with causal inference methods in this context: (i) the dimensionality of the covariates and (ii) practical violations of the positivity assumption. The dimensionality of the covariates presents challenges to estimation. The discussion of exchangeability above highlights the need to potentially adjust for a high-dimensional W j , with a limited number of sequences available for a given analysis. Thus, we may be motivated to consider flexible machine learning algorithms that are built specifically for modeling high-dimensional covariates. It would be natural to couple these modeling strategies with doubly robust methods for inference, as these methods typically allow for regular, parametric-rate inference, under certain statistical assumptions. However, these methods can be susceptible to erratic behavior in the presence of practical positivity violations. One approach to improving behavior of doubly robust estimators in the presence of positivity violations involves a careful pre-selection of adjustment covariates, informed by background knowledge [17]. However, this would appear difficult in the present setting due to imperfect understanding of the relationships between AA residues on the Env protein and their relationship with antibody neutralization. Moreover, for a method to be successfully employed to scan the entire Env protein for significant residues, it should be a high throughput method that is capable of delivering stable inference on hundreds of AA residues. An analysis that involves deliberate pre-selection of covariates AA residue-by-residue is infeasible. Therefore, we are motivated to develop an automated procedure for covariate selection in this context. In the subsequent sections, we propose a solution for this problem using CTMLE.

2.4 TMLE

Our CTMLE builds on targeted minimum loss-based estimators (TMLE) of the average treatment effects [18]. TMLEs are constructed using estimates of two key quantities: the OR and the GPS. The OR is defined as the conditional mean of the outcome given the treatment and confounding factors, which we denote by Q ¯ ( W j , W − j ) = P ( Y = 1 ∣ W j , W − j ) . The GPS is the conditional distribution of AAs at residue j given all other AA residues and sequence information, which we denote by g j ( w ∣ W − j ) = P ( W j = w ∣ W − j ) for w ∈ W j . A TMLE for estimating E [ Y ( W j = w ) ] for a particular j and w can be constructed in two steps. In step 1, initial estimator Q ¯ n of the OR can be obtained using any learning technique for classification of a binary outcome. Similarly, the estimator g n , j of GPS could be estimated using any multi-class classification technique. In step 2, a single iteration of boosting is used to update the initial OR. For a given OR estimate, Q ¯ , we can compute its empirical risk,

Next, we update Q ¯ n by minimizing empirical loss (1) along a particular univariate regression model indexed by Q ¯ n . For each w , we define the logistic regression model

We can obtain a maximum likelihood estimate ε n ( w ) = arg min ε ( w ) L n , j ( Q ¯ n , ε ( w ) ) and define the updated OR estimate

The final estimate is computed as

The key idea motivating the second step of the TMLE is that the bias of the revised OR estimator Q ¯ n ∗ is generally smaller than that of Q ¯ n with respect to estimating p j ( w ) .

This procedure can be repeated for each w ∈ W j . With a minor abuse of notation, we define the implied updated estimate of E { Y ∣ W j , W − j } as

Recalling our null hypothesis of interest that E [ Y ( W j = w ) ] is constant in w , a TMLE-based test based on these estimates can be derived as follows. For w ∈ W j , we define

and let D j , i be a ∣ W j ∣ -length row vector consisting of D j , i ( w ) for each w ∈ W j . Let D j be a n × ∣ W j ∣ matrix formed by stacking the row vectors D j , i , i = 1 , … , n . The estimates p ˆ j = { p ˆ j ( w ) : w ∈ W j } of p j = { p j ( w ) : w ∈ W j } have asymptotic covariance that is consistently estimated by Σ ˆ = n − 1 D j ⊤ D j .

We propose to test H 0 using a ∣ W j ∣ 2 − 1 degree-of-freedom Wald test. Specifically, we define a contrast matrix A , where each row defines a contrast between p ˆ j ( w ) and p ˆ j ( w ′ ) for w , w ′ ∈ W j . For example, if ∣ W j ∣ = 4 , then

and our null hypothesis can be written as H 0 : A p j = 0 .

If instead inference on a single p j ( w ) is desired, the standard error of p ˆ j ( w ) is given by the square root of the ( j , j ) element of Σ ˆ . Assuming n is large enough, a two-sided Wald-style hypothesis test that rejects H 0 : p ˆ j ( w ) = β whenever ∣ p ˆ j ( w ) − β ∣ / σ ˆ n is greater than the 1 − α / 2 quantile of a standard normal distribution will have type I error no larger than α .

2.5 CTMLE

A challenge in utilizing TMLE in the present setting is that the estimated GPS may be extremely small for some virus sequences due to high correlations between residues in the Env protein. This can lead to erratic estimates ε n ( w ) of ε ( w ) , which can degrade the OR estimate Q ¯ n ∗ in (2) considerably. The resulting ATE estimate could be highly biased, with correspondingly inflated type I errors and poor confidence interval coverage. There are fixes available [17]; however, the approaches require manual manipulation from the analysts (e.g., to identify overlapping covariate regions and redefine the causal parameter accordingly). Unfortunately, in our motivating example we wish to consider hundreds of AA residues. This motivates the pursuit of a more high-throughput analysis approach that maintains stability even in rather extreme scenarios.

Recently, data-driven PS model-building strategies have emerged [14,13]. A key insight of these approaches is that PS models should adjust only for variables that are related to the outcome. Thus, estimated ORs can be used to inform variable inclusion in propensity score models. By appropriately screening so-called instrumental variables (those that impact only the GPS but not the OR), we may generate less extreme estimates of the GPS and thereby attain more stable behavior of TMLE.

In particular, we propose using variable importance measures from an OR model to select variables to include in the GPS model. We focus on random forest models explicitly, for estimation of both the OR and GPS, though the methods theoretically extend to any machine learning framework that has associated variable importance measures. We focus on random forests in particular for two reasons: (i) robust and fast software implementations with variable importance measures are readily available; (ii) past work has shown that random forests perform extremely well in predicting antibody neutralization, even when compared with other state-ofthe-art approaches like deep learning [11,19].

Our approach entails first using random forests to estimate the OR. During the random forest construction, variable importance for each AA residue is summarized by the mean decrease in Gini index. The Gini index quantifies how much each feature contributes to the decline of node impurities on average. We select a fixed number of the top-ranked features to include in the GPS model. The OR and GPS are then combined into estimates of E [ Y ( W j = w ) ] for each j and w and a test of the significance of each AA residue. We propose a CTMLE algorithm for selecting the optimal number of features to include in the GPS model.

3.1 Design

We conducted Monte-Carlo simulation studies to compare the performance of the proposed CTMLE to several implementations of TMLE. The first was a standard implementation of TMLE, where all features were included in the GPS model, while the other TMLEs used reduced numbers of features (either 5, 10, 50, or 100) based on their estimated importance in the OR model. For all estimators, GPS estimates were truncated at 0.01. We considered two sample sizes n ∈ { 500 , 1,000 } , which are similar to the sample sizes of the CATNAP data. For each sample size, 1,000 data sets were generated as follows. To represent a sequence of AA residues, we simulated 200 moderately correlated categorical variables, W 1 , W 2 , … , W 200 , each containing four levels, W j ∈ { 1 , 2 , 3 , 4 } . To generate a realization of W = ( W 1 , … , W 200 ) ⊤ , we drew a random variable X = ( X 1 , … , X 200 ) ⊤ from a multivariate normal distribution N ( μ 0 , Σ 0 ) with μ 0 = ( 0 , 0 , 0 , 0 ) ⊤ and an autoregressive-1 covariance structure, Σ 0 = σ 2 Σ ( ρ ) , where Σ ( ρ ) is a matrix with ( j , k ) th entry equal to ρ ∣ j − k ∣ . In our example, σ and ρ were fixed at 0.9 and 0.75, indicating a moderate correlation between adjacent features. Next, we set W j = 1 if X j < q 0.25 , W j = 2 if q 0.25 < X j < q 0.50 , W j = 3 if q 0.50 < X j < q 0.75 , and W j = 4 if X j > q 0.75 , where q p is the p th quantile of a standard normal random variable. Five true signals were randomly selected among all features ( W j : j ∈ J , where J = { 37 , 87 , 94 , 135 , 151 } ). Given W = W i , the outcome Y i was generated from a Bernoulli distribution with success probability = expit ∑ j ∈ J ∑ ℓ = 1 4 β j l I ( W j , i = ℓ ) (Table 2).

We focused on evaluating operating characteristics for three specific AA residues: (i) position 10, a position unrelated to the outcome and essentially uncorrelated with any true signals; (ii) position 85, a position unrelated to the outcome, but highly correlated with a true signal; and (iii) position 87, a position that is truly related to the outcome ( p 87 ( 1 ) = 0.514 , p 87 ( 2 ) = 0.588 , p 87 ( 3 ) = 0.342 , p 87 ( 4 ) = 0.545 ). For each of these sites, the C/TMLE point estimates were evaluated by their bias and mean-squared error (MSE) and we compared the power of the Wald test for rejecting the null hypothesis of no relationship between AAs and outcome.

3.2 Results

We found that the bias and MSE of the TMLE estimators tended to increase as more variables were included in the GPS model (Figure 1, left/middle column), with the largest bias seen for the true signal (position 87, blue plus). The CTMLE had considerably reduced bias and MSE relative to the standard TMLE that included all 200 features in the GPS model. Moreover, in the two null scenarios (positions 10 and 85), the standard TMLE had drastically inflated type I error (Figure 4, right column). For example, the tests incorrectly rejected the null hypothesis more than 60% of the time when n = 1,000 . On the other hand, CTMLE-based tests appropriately rejected the null hypothesis about 5% of the time for the null positions. For the position exhibiting a true signal (position 87), CTMLE rejected the null hypothesis 57% of the time at n = 500 and 82% of the time at n = 1,000 .

Figure 4

Bias, MSE, type I error rate for null cases and power for true signals. Simulation results with sample size of 500 (first row) and 1,000 (second row) were visualized for two null cases and one true signal. Compared with standard TMLE approach using all 200 features, CTMLE largely reduced bias and MSE, especially for the non-signals. it also maintained a relatively high power with better controlled type I error as good as the best model among the five proposed TMLE scenarios. As sample size increases, the power for CTMLE increased from 0.57 to 0.82.

3.3 Comparison with non-causal approach

We also compared our CTMLE-based test for identifying residues causally related to antibody neutralization to a test developed for random forest-based importance measures. Using the simulation design above, we compared the tests’ type I error rates (proportion of tests in which the null hypothesis was rejected for W 10 and W 85 ) and power under the alternative (proportion of tests in which the null hypothesis was rejected for W 87 ). We found that both tests controlled the type I error rate for the uncorrelated residue ( W 10 ); however, for the correlated residue ( W 85 ), the type I error for the random forest-based test was inflated relative to 0.05. The inflation remained at larger sample sizes. On the other hand, the random forest-based test exhibited higher power under the alternative ( W 87 ) (Table 3).

4.1 CATNAP datasets

The proposed methods were applied to the Compile, Analyze and Tally NAb Panels (CATNAP) database [7]. This database contains information on HIV virus sequences including the HIV envelope AA sequence (the features of interest), other key variables describing structural and biological information (e.g., geographic origin), and a binary measure of whether the virus can be neutralized by a particular antibody (our outcome of interest). We analyzed five antibodies: VRC01, VRC26.08, PGT145, PGT121, and 10-1074. For each antibody, a pre-screening step was applied to exclude AA positions that were nearly constant across all virus sequences. In practice, some residues in Env may exhibit little or no variation. For these residues, we wish to avoid hypothesis testing since there is no hope of garnering sufficient statistical evidence for signal detection. In particular, at each position, AAs that were observed in fewer than 30 total virus sequences were combined into a single AA category. After the combination, only residues still with smallest level count greater than 30 were retained for analysis. Since the number of remaining AA residues that pass this variable screen procedure is still large, we used a Bonferroni correction to control the chance of falsely rejecting null hypotheses. Namely, we tested each individual residue at a significance level of α / N A A ( post-screening ) , where N A A ( post-screening ) is the number of AA residues passing the screening above.

4.2 Results

The numbers of AA residues that passed the screening process are summarized in Table 4 and the plots of − log 10 ( p -value) for each residue are shown in Figure 5. For each antibody, the figure highlights functional regions of the gp120 protein, which are putatively associated with antibody activity. Our analysis revealed that many significant residues fell within known regions – the V1/V2 loop in particular contained many significant residues – however, there were also many significant residues in gp41, particularly for VRC01.

Figure 5

Significant AA residues for five antibodies. Orange highlights denote V1/V2/V3 loop regions that are putatively associated with these antibodies; purple regions highlight additional functional regions of the gp120 protein. Gray points denote residues that did not pass variability screening. (a) VRC01, (b) VRC26.08, (c) 10-1074, (d) PGT121 and (e) PGT145.

4.3 Comparison with other approaches

We compared our results for the VRC01 antibody to the LASSO and random forest approach by feature selection described above (Figure 6). For LASSO-based approach, there were more than 15 multi-level features selected in the CV-selected model, so the inference from the GLM in the held-out data was highly unstable. This illustrates a potential difficulty of this approach, which requires the analyst to make arbitrary, post-hoc modeling decisions to stabilize inference. On the other hand, our proposed CTMLE provides a fully automated procedure that balances type I and II errors.

The random forest-based approach provided greater stability and identified 76 significant residues. However, less than half (16/76) were sites in gp120, which has been identified through other experiments as the likely region of importance for neutralization by VRC01. This result raises the possibility of a high false positive rate of the random forest procedure in this context. Additionally, the associated variable importance measures have dubious relevance in terms of their scientific interpretation. In contrast, our CTMLE is able to provide biologically meaningful interpretations of estimated parameters.

Figure 6

Results for the VRC01 antibody based on three approaches: LASSO, random forests (RF), and CTMLE. The dashed lines indicate thresholds from a Bonferroni-correction to compensate multiple comparisons.