原始论文：论文地址

摘要：

动机：

HIV-1 antiviral resistance is a major cause of antiviral treatment failure. The in vivo fitness landscape experienced by the virus in presence of treatment could in principle be used to determine both the susceptibility of the virus to the treatment and the genetic barrier to resistance. We propose a method to estimate this fitness landscape from cross-sectional clinical genetic sequence data of different subtypes, by reverse engineering the required selective pressure for HIV-1 sequences obtained from treatment naive patients, to evolve towards sequences obtained from treated patients. The method was evaluated for recovering 10 random fictive selective pressures in simulation experiments, and for modeling the selective pressure under treatment with the protease inhibitor nelfinavir.

HIV-1病毒抗药性是抗病毒治疗失败的主要原因。理论上，病毒在治疗存在的情况下体验到的体内适应性景观可以用来确定病毒对治疗的敏感性和抗药性的遗传障碍。我们提出了一种方法，通过反向工程确定治疗前的HIV-1患者所需的选择性压力，使其发展到得到治疗的患者的序列，从而从不同亚型的临床遗传序列数据中估计这种适应性景观。该方法在模拟实验中被用来恢复10种随机虚构的选择性压力，并用于模拟使用蛋白酶抑制剂奈非那韦治疗下的选择性压力。

结果：

The estimated fitness function under nelfinavir treatment considered fitness contributions of 114 mutations at 48 sites. Estimated fitness correlated significantly with the in vitro resistance phenotype in 519 matched genotype-phenotype pairs ( $R^2 = 0.47(0.41-0.54)$ ) and variation in predicted evolution under nelfinavir selective pressure correlated significantly with observed in vivo evolution during nelfinavir treatment for 39 mutations (with $FDR = 0.05$ ).

在奈非那韦治疗下估计的适应性函数考虑了48个位点上114个突变的适应性贡献。在519个匹配的基因型-表型对中，估计的适应性与体外抗药性表型显著相关（ $R^2 = 0.47（0.41 － 0.54）$ ）。在奈非那韦选择压力下预测的进化变化与在奈非那韦治疗期间观察到的体内进化在39个突变中显著相关（假阳性发现率FDR = 0.05）。

正文部分

1.介绍

HIV antiviral drugs interfere with viral proteins resulting in the inhibition of HIV replication. In many cases, HIV escapes the inhibition of these drugs by selection of drug resistance mutations, leading to treatment failure (Vandamme, 1999). To combine drugs in an effective treatment therefore requires taking into account the presence of resistance mutations, and resistance testing has become a standard of care (Vandamme et al., 2004). Different viral mechanisms may be distinguished that affect short-term versus long-term response to antiviral treatment. In addition to the impact of other factors such as adherence, potency of therapy and pharmacokinetics, the short-term response to treatment is mainly determined by the susceptibility of the virus to the drugs. In the long-term, susceptible virus may evolve to acquire resistance mutations, and the expected time needed for the virus to evolve the necessary resistance mutations is related to the number of nucleotide substitutions required, which can be quantified as the genetic barrier. Several bioinformatics methods have been used successfully in the field of antiviral drug resistance, including methods that predict in vitro phenotypic resistance from the genetic sequence, and methods that describe qualitative relationships between different mutations selected during treatment (Beerenwinkel et al., 2005b). Recently, these techniques were combined to compare the genetic barrier for individual drugs versus drug combinations (Beerenwinkel et al., 2005c).

HIV抗病毒药物通过干扰病毒蛋白质，从而抑制HIV的复制。在许多情况下，HIV通过选择抗药性突变逃避这些药物的抑制，导致治疗失败(Vandamme, 1999)。因此，有效地将药物组合在一起进行治疗需要考虑到抗药性突变的存在，抗药性测试已成为常规护理的一部分(Vandamme等人，2004)。可以区分影响抗病毒治疗短期和长期反应的不同病毒机制。除了其他影响因素，如依从性、治疗的效能和药物代谢动力学，治疗的短期反应主要由病毒对药物的敏感性决定。长期来看，敏感的病毒可能会演变成具有抗药性突变的病毒，病毒演变出必要的抗药性突变所需的预期时间与所需的核苷酸替代数量有关，这可以量化为遗传障碍。几种生物信息学方法在抗病毒药物抗药性领域成功应用，包括从遗传序列预测体外表型抗药性的方法，以及描述在治疗期间选择不同突变之间定性关系的方法(Beerenwinkel等人，2005b)。最近，这些技术被组合起来比较单一药物与药物组合的遗传障碍(Beerenwinkel等人，2005c)。

Due to technical shortcomings, problems with the interpretation of the results, and the lack of a genetic barrier concept, in vitro phenotypic assays display limitations in their capacity of predicting therapy outcome (Van Laethem and Vandamme, 2006). Therefore, the usefulness of machine-learning approaches to predict resistance phenotype from genotype may be limited. On the other hand, the success of attempts to directly learn genotypic patterns responsible for reduced treatment response from clinical data has been limited by the lack of sufficient data, and the confounding effect of many other factors (DiRienzo and DeGruttola, 2002). The in vivo fitness of the virus in presence of treatment, which reflects both effects of drug resistance and replication capacity,1 determines the immediate treatment response, but cannot be measured directly.

由于技术上的缺陷、对结果解释的问题，以及遗传障碍概念的缺乏，在体外表型测定法在预测治疗效果的能力上表现出限制（Van Laethem和Vandamme，2006年）。因此，使用机器学习方法根据基因型预测抗药性表型的有效性可能有限。另一方面，直接从临床数据中学习导致治疗反应减弱的基因型模式的成功尝试，受到了数据不足和许多其他因素的混淆效应的限制（DiRienzo和DeGruttola，2002年）。病毒在治疗存在的情况下的体内适应性，反映了药物抗性和复制能力的两种影响，决定了即时的治疗反应，但不能直接测量。

注释：

In this article, the universal meaning of fitness as the capacity to replicate in a given environment is used, rather than as a synonym for replication capacity in a drug-free environment as is often the case in the HIV drug resistance community.

在这篇文章中，我们使用适应性有一个普遍的含义，即在给定环境下的复制能力，而不是作为在无药环境下的复制能力的同义词，这在HIV药物抗性社区中常常是这样。

However, HIV tries to recover its ability to replicate efficiently in presence of treatment by accumulating resistance mutations, thus exploring sequence space in the immediate neighborhood of the current sequence. Therefore, observed evolution in clinical sequences at treatment failure provides information about the fitness landscape, but only in the immediate neighborhood of the current sequence. In this article, we present a method to reverse engineer this fitness landscape experienced by HIV-1 in presence of treatment as a function of the genetic sequence, from observed selection during treatment in clinical sequences. The method searches for a fitness landscape, which explains how an observed population of treated sequences could have evolved from a population of untreated sequences under selective pressure. After showing that random, but known fitness functions could be successfully estimated in this way, we applied the method to model the fitness landscape of HIV-1 in presence of the protease inhibitor (PI) nelfinavir (NFV).

然而，HIV试图通过积累抗性突变来恢复其在治疗存在的情况下有效复制的能力，从而在当前序列的紧邻邻域内探索序列空间。因此，治疗失败时临床序列的观察进化提供了有关适应性景观的信息，但只在当前序列的紧邻邻域内。在这篇文章中，我们提出了一种方法，通过观察治疗期间在临床序列中的选择，反向工程这种HIV-1在治疗存在的情况下体验到的适应性景观作为遗传序列的一个函数。这种方法寻找一个适应性景观，解释了一个观察到的治疗序列的群体是如何在选择压力下从一个未治疗序列的群体演变过来的。在显示出以这种方式能够成功估计随机但已知的适应性函数之后，我们应用了这种方法来模拟在蛋白酶抑制剂（PI）奈非那韦（NFV）存在的情况下HIV-1的适应性景观。

2.材料与方法

2.1 临床数据集

To estimate the fitness function under NFV selective pressure, clinical data was pooled from the Stanford HIV Drug Resistance Database (Kantor et al., 2001), from the University Hospitals, Leuven, Belgium, and from Hospital Egas Monis, Lisbon, Portugal, to create a treated population PT of 1026 sequences from patients with experience to NFV as sole PI, and a naive population PN of 7774 sequences from PI naive patients. At most one treated and one naive sequence per patient was used, and duplicate sequences (that were present in the hospital database but also published in the Stanford database) were identified and removed. The treated population consisted mostly of HIV-1 subtype B sequences, but included also a large number of subtype G and subtype C sequences (Fig. 1), as determined from the protease and partial reverse transcriptase sequences using the REGA HIV-1 Subtype tool v2.0 (de Oliveira et al., 2005). Nucleotide ambiguities that occur commonly in the population sequences were resolved by randomly substituting the mixture with a suitable pure nucleotide. Using a threshold of 0.5% prevalence at variable sites, 114 mutations at 48 protease positions were included in the fitness function models, and are listed in Supplementary Material A. To remove redundancy, the most prevalent mutation at each position was considered the ‘wild type’ and was omitted from the fitness function model, as its presence was implied in the absence of any mutation.

为了估计在奈非那韦选择压力下的适应性函数，我们将斯坦福大学HIV药物抗性数据库（Kantor等人，2001）的临床数据、比利时鲁汶大学医院的数据以及葡萄牙里斯本埃加斯莫尼兹医院的数据汇集在一起，创建了一个经历过以奈非那韦为唯一蛋白酶抑制剂治疗的患者序列的治疗群体PT，人数为1026，并且创建了一个由7774个PI初级患者序列组成的初级群体PN。每个患者最多使用一个治疗序列和一个初级序列，且重复序列（存在于医院数据库但也在斯坦福数据库中公布的序列）已被识别并移除。治疗群体主要由HIV-1亚型B序列组成，但也包括大量的亚型G和亚型C序列（图1），这是通过使用REGA HIV-1亚型工具v2.0（de Oliveira等人，2005）从蛋白酶和部分逆转录酶序列中确定的。在群体序列中常见的核苷酸不确定性通过随机用合适的纯核苷酸代替混合物来解决。使用变化位点0.5%的患病率阈值，在适应性函数模型中包括了48个蛋白酶位置的114个突变，这些突变在补充材料A中列出。为了消除冗余，在每个位置上最常见的突变被认为是“野生型”，并且在适应性函数模型中被省略，因为在没有任何突变的情况下，其存在被暗示。

2.2估计适应度函数的方法

We present a method with objective to learn a functio $F(A_1,...,A_n)$ , where $A_i$ presents presence or absence of a mutation, that represents the fitness landscape of HIV under drug selective pressure. To learn F, we find a function that fits with the evolution of the virus in a naive population of patients $P^N$ to a treated population $P^T$ , and is closest to neutrality (minimizing $| F-1|$ ). The fitness function F incorporates interactions indicated using Bayesian Network (BN) learning, and its parameters are estimated using an iterative procedure where evolution for $P^N$ over the current fitness function estimate is simulated, and compared to $P^T$ .

我们提出了一种方法，其目标是学习一个函数 $F(A_1,...,A_n)$ ，其中 $A_i$ 表示突变的存在或缺失，代表HIV在药物选择压力下的适应性景观。为了找到F，我们找到一个函数，适合从处于初始阶段的病人群体 $P^N$ 到接受治疗的病人群体 $P^T$ 的病毒进化，并且与中性（最小化|F-1|）最接近。适应性函数F采用贝叶斯网络（BN）学习来指示交互作用，其参数使用迭代过程进行估计，在当前适应性函数估计下，模拟 $P^N$ 的演化并与 $P^T$ 进行比较。

2.2.1适应度函数的结构

The protease amino acid sequences from the treated population PT were used to learn interactions between mutations as described before (Deforche et al., 2006). Briefly, a data set was created where a boolean variable indicated the presence of each included mutation. BN structure learning (Myllymaki et al., 2002) on this boolean data was used to discover relationships between these mutations that may indicate epistatic fitness effects. By assuming conditional independencies, the Bayesian Network refactors the Joint Probability Distribution (JPD) in a product of Conditional Probability Distributions (CPD), leading to a reduction in number of parameters to model the JPD. Formally, for n variables $A_1,..,A_n$ (representing amino acid mutations), we would write:

$P(A_1,...,A_n) = \prod_i^n P(A_i | parents(A_i))$

我们使用接受治疗的患者群体PT中的蛋白酶氨基酸序列来学习突变之间的交互，就像之前描述的那样(Deforche等人，2006）。简单地说，创建了一个数据集，其中一个布尔变量指示了每个包含的突变的存在。对这个布尔数据使用贝叶斯网络结构学习(Myllymaki等人，2002)来发现这些突变之间的关系，这可能表明表现性适应性效应。通过假设条件独立性，贝叶斯网络将联合概率分布(JPD)重构为条件概率分布(CPD)的乘积，导致模拟JPD参数数量的减少。正式地，对于n个变量A_1...,A_n(表示氨基酸突变），我们可以写作：

$P(A_1,...,A_n) = \prod_i^n P(A_i | parents(A_i))$

在这里，P(A | B)表示已知B条件下A的概率，parents(Ai)表示在贝叶斯网络结构中的变量Ai的父节点。我们称治疗群体PT中的氨基酸序列的最可能网络，其结构为ST，CPD参数为� T，为BNT(� T, ST)。

我们以与BNT重构JPD相同的方式对相对适应度函数F（A1，…，An）进行建模：

with parents (Ai) the parents in ST, and F(A |B) the Conditional Fitness Contribution (CFC) of the presence of A, depending on the presence of B. The assumption here is that if two mutations are synergistic for example, they would occur more often together than not, and a dependency should be visible in the JPD too. See Supplementary Material B for an example.

其中，parents(Ai)指的是ST中的父节点，F(A |B)表示A的存在依赖于B的存在的条件适应性贡献（CFC）。此处的假设是，如果两个突变是协同的，例如，它们会比不在一起更常见，因此在JPD中也应该可以看到依赖关系。请参阅补充材料B以获取示例。

The CPDs are modeled by specifying the probability for a mutation Ai given any pattern of parent mutations k, in Conditional Probability Tables (CPTs): � i,k ¼ P(Ai ¼ 1| parents(Ai) ¼ k). Similarly, we used Conditional Fitness Tables (CFTs) to model the CFCs for each mutation Ai, which specify a different fitness contribution of the presence of a mutation Ai for every pattern of parent mutations: � i,k ¼ F(Ai ¼ 1| parents(Ai) ¼ k).

CPD通过指定在任何父突变模式k下突变Ai的概率，在条件概率表（CPTs）中：� i,k ¼ P(Ai ¼ 1| parents(Ai) ¼ k)。类似地，我们使用条件适应性表格（CFTs）来模拟每个突变Ai的CFCs，这些表格为每种父突变模式指定突变Ai的存在的适应性贡献：� i，k ¼ F(Ai ¼ 1| parents(Ai) ¼ k)。

2.2.2 对进化的建模

Both for estimating the fitness function parameters, and for prediction of sequence evolution during treatment, the same stochastic model of HIV evolution was used. Evolution was considered as an accumulation of fixations of nucleotide mutations in the HIV intra-host population, as reflected in the consensus sequence, under the selective pressure of an arbitrarily complex fitness function. This corresponds roughly to how HIV resistance evolution is observed in population sequences obtained by genotypic resistance tests (Van Laethem and Vandamme, 2006).

无论是在估计适应性函数参数还是在预测治疗期间的序列演化时，都使用了相同的HIV演化随机模型。演化被认为是在复杂性任意的适应性函数的选择压力下，在HIV宿主内部群体中固定核苷酸突变的累积，这反映在一致性序列中。这大致对应于通过基因型抗性测试获得的群体序列中观察到的HIV抗性演化(Van Laethem和Vandamme，2006年)。

The HIV intra-host population was modeled by a finite ideal Wright– Fisher population with selection and mutation, using empirical estimates of the HIV intra-host effective population size, mutation rate and mutation rate biases derived from literature, and selection coefficients derived from the fitness function F. Analytical results lack for fixation time distributions of mutations in the Wright–Fisher model for all but the simplest cases (Ewens, 1979). Therefore, to sample from these distributions, for a Wright–Fisher model with multiple loci and a complex fitness function with epistatic interactions, an approximate simulation of this model was implemented. For a detailed description of the implemented model and approximations see Supplementary Material C.

HIV内宿主种群被模拟为一个有限的理想的Wright–Fisher种群，其中包含选择和突变，使用的是从文献中得到的HIV内宿主有效种群大小、突变率和突变率偏差的实证估计，以及从适应性函数F派生的选择系数。对于除了最简单的情况以外的Wright–Fisher模型中的突变固定时间分布，缺乏分析结果(Ewens，1979)。因此，为了从这些分布中抽样，对于具有多个位点和复杂适应性函数（其中包含表现性交互作用）的Wright–Fisher模型，实现了这个模型的近似模拟。关于已实施模型和近似的详细描述，请参阅补充材料C。

2.2.3 适应度函数参数 Fitness function parameters

The parameters i, k of the function F are estimated so that evolution over the fitness landscape of a naive population PN resembles the treated population PT. Therefore, evolution is simulated for sequences sampled from the naive population PN using the fitness function, to obtain an evolved population PE. The difference between the sequence populations PE and PT, which must thus be minimized, is measured by comparing the parameters of BNT(� T, ST) of the treated data set, with BNE (� E, ST), a BN estimated from the simulated population using the structure that was learned from the treated data set. Thus, we measure and minimize the difference in prevalence of each mutational pattern that is modeled by the BN, and for which the fitness function specifies a separate fitness contribution.

函数F的参数i，k被估计为使得初级群体PN在适应性景观上的演变类似于治疗群体PT。因此，模拟使用适应性函数进行初级群体PN采样的序列的演变，以获得演化群体PE。通过将治疗数据集的BNT(� T, ST)的参数与从模拟群体估计出来的、使用从治疗数据集学习到的结构的BN BNE(� E, ST)，从而度量并最小化序列群体PE和PT之间的差异，该差异因此必须被最小化。因此，我们度量并最小化了BN模拟的每种突变模式的患病率的差异，并且适应性函数为其指定了一个单独的适应性贡献。

Given this minimization objective, the parameters � i, k are not necessarily unique, since we cannot quantify how unfit unobserved mutational patterns are. Indeed, we can only determine how unfit these patterns must be at least to explain their lack of evolution. We constrain the search to a unique solution by minimizing | � i, k � 1 | for each parameter � i, k, as a secondary objective (with a low weight compared to the first objective).

考虑到这个最小化目标，参数� i, k并不一定是唯一的，因为我们无法量化未观察到的突变模式有多不适应。实际上，我们只能确定这些模式至少要多么不适应才能解释它们缺乏演化。我们通过最小化每个参数� i，k的| � i, k � 1 |，作为二级目标（与第一目标相比权重低），将搜索约束到唯一解。

An iterative algorithm was used to estimate the parameters. Thealgorithm starts initially from a flat fitness function (i.e. F(Ai, ..., A,)1),and in each iteration this function is updated in a step-wisefashion. A population pE is computed using the current estimate of thefitness function F. The fitness function parameters p,k are subsequentlyadjusted based on the difference in the sufficient statistics (which reflectthe counts) related to the BN parameters 0,k and 0,k: p,k is increasedwith a small multiplicative factor (l+,k) when there is too few ofmutation A,for parent combination k in the evolved versus treatedpopulation, or vice versa if there is too much of mutation A,. By usingthe sufficient statistics instead of the actual CPD parameters, uncertainty on these parameters is taken into account. The step sizes d, k aredynamically adjusted depending on the convergence of the corresponding pik. Details, pseudo-code, and convergence properties of thealgorithm are presented in Supplementary Material D.

使用迭代算法来估计参数。该算法最初从一个平坦的适应性函数开始（即 F(Ai，..., A,) = 1），在每次迭代中，这个函数以逐步的方式进行更新。使用当前的适应性函数F的估计计算出一个种群pE。接着，基于与BN参数 0,k 和 0,k 相关的足够统计量（反映了计数）的差异，随后调整适应性函数参数 p,k：当在演化的种群与处理的种群中，对于父组合k，突变 A, 的数量过少，或者相反，如果突变 A, 的数量过多，p,k 就以一个小的乘法因子（l +,k）增加。通过使用足够的统计量代替实际的CPD参数，考虑了这些参数上的不确定性。步长 d, k 根据相应的 pik 的收敛情况动态调整。算法的细节、伪代码和收敛性质在补充材料D中进行了介绍。

The amount of evolution experienced by HIV under treatment depends on many factors such as the baseline viral load, the potency of the combination therapy to suppress residual replication (which depends on the amount of resistance), patient adherence and duration

of the therapy. The probability distribution for the effective number of generations under drug pressure between sequences from drug naive and treated patients P(GT), which is used to create PE, is in general unknown and was assumed to be uniform in the interval [0, Gmax], with Gmax a maximum limit for the number of generations. Variation of Gmax does not affect the shape of the landscape, but its steepness (lower values result in a steeper landscape).

HIV在治疗下经历的演变程度取决于许多因素，如基线病毒载量，组合疗法抑制残余复制的效能（取决于抗药性的程度），患者的依从性和疗程的持续时间。在药物初级和治疗患者的序列之间，在药物压力下的有效代数的概率分布P(GT)，用于创建PE，通常是未知的，并且被假定为在区间[0, Gmax]内是均匀的，其中Gmax是代数的最大限制。Gmax的变化不会影响景观的形状，但会影响其陡度（较低的值会导致更陡峭的景观）。

2.2.4 Phylogenetic guide tree系统发育导向树

The sampling of sequences fromDNwas guided by a phylogenetic tree, and more weight was given tosequences from the naive population that were epidemiologically linkedto the treated population. This assures that the sampled population hada similar epidemiological background to the treated populationavoiding that mutations linked to epidemiologies with a differentdistribution among these population were assigned as arising duringtreatment fan improvement compared to stratifying according toepidemiology (Deforche et al., 2006; Kantor et al., 2006)]. The proteaseand partial reverse transcriptase nucleotide sequences were used to reconstruct a neighbor-joining phylogenetic tree including all naive andtreated isolates used in the training data. The tree was built using PAUP(Swofford, 2000)using the HKY-y substitution model, and codonsrepresenting IAS resistance associated positions (Johnson et al., 2005)were excluded to avoid problems of convergent evolution.
从DN抽取序列是由一个系统发育树来指导的，对与治疗的种群有流行病学联系的初始种群的序列赋予了更多的权重。这保证了被抽样的种群具有与治疗种群相似的流行病学背景，避免了与这些种群在流行病学上有不同分布的突变被标记为在治疗期间产生[相比按照流行病学进行分层（Deforche et al., 2006；Kantor et al., 2006）这是一种改进]。使用蛋白酶和部分逆转录酶核苷酸序列来重建包括所有初级和治疗分离物在内的邻接连接系统发育树，这些分离物在训练数据中使用。利用PAUP（Swofford，2000）和HKY-y替代模型构建树，并排除表示IAS抗性相关位置（Johnson et al., 2005)的密码子，以避免收敛演化的问题。

2.2.5 Constants 常数

2.3 验证实验

2.3.1 与奈非那韦耐药性表型的相关性

Using apublic data set of matched genotypephenotype pairs for subtype Bsequences ffrom the Stanford HIV Drug Database (Kantor et al.2001)], estimated fitness was compared to in vitro resistance fold changephenotype. Sequences with unknown amino acid mutations ('Z’or '’)were removed, as were sequences with a fold change at the upperdetection range of the assay. For each of the remaining 5l9 amino acidsequences j, fitness f, was estimated from resistance fold change Rusing (Holford and Sheiner, 1982):

使用来自斯坦福HIV药物数据库（Kantor等人，2001年）的匹配基因型和表型配对的公开数据集针对B型序列，将估计的适应性与体外抗性倍数变化表型进行比较。移除了具有未知氨基酸突变（'Z’或'’）的序列，以及折叠变化在测定的上限检测范围内的序列。对于剩下的519个氨基酸序列j，根据抗性倍数变化R使用(Holford和Sheiner，1982)的公式估计适应性f：

with e,the replication capacity of the virus and D the effective drugconcentration.The values e,are generally unknown, and e;= 1 wasassumed for all strains when computing f. The fitness estimated fromthe phenotypes was then compared with the fitness computed using theestimated fitness landscape by computing the correlation coefficient,which was indifferent to the valule of D.

其中e表示病毒的复制能力，D表示有效的药物浓度。e的值通常是未知的，当计算f时，假定所有菌株的e都等于1。然后，通过计算相关系数，将从表型估计出的适应性与使用估计的适应性景观计算出的适应性进行比较，这是不受D值影响的。

2.3.2 与观察到的奈非那韦抗性演化的相关性

In 404 patients for which a baseline and consecutivefollow-up sequence during NFV treatment was available, the accuracyof the model to predict observed resistance evolution was evaluatedThese pairs were independent from the cross-sectional training data. Foreach wild type and mutation included in the fitness function, correlationof observed evolution (0 or 1) with predicted evolution (0<p<1) wasevaluated.Correlation with observed evolution was analyzed with alinear model, which included next to the predicted evolution a non-linearcorrection for the number of observed substitutions for each sequenceCorrection for multiple testing was done using the Benjamini andHochberg method with FDR = 0.05. An observed mixture (such asL63LPA) was evaluated against predicting evolution of the mutations(for this example L63P or L63A). Prediction of loss of mutations inmixtures (such as LPA63P) was not considered.

在404名有基线和连续随访序列在奈非那韦治疗期间可用的患者中，评估了该模型预测观察到的抗性演化的准确性。这些配对与横断面训练数据无关。对于适应性函数中包含的每个野生型和突变，评估了观察到的演化（0或1）与预测的演化（0 < p < 1）的相关性。使用线性模型分析与观察到的演化的相关性，该模型除了预测的演化外，还包括对每个序列的观察到的替代数的非线性校正。使用Benjamini和Hochberg方法（FDR = 0.05）进行多重测试的校正。对观察到的混合物（如L63LPA）进行评估，以预测突变的演化（对于这个例子，L63P或L63A）。没有考虑预测混合物中突变的丧失（如LPA63P）。

3. 结果

3.1 模拟实验

To illustrate convergence properties of the method, the method was tested in two series of simulation experiments: (i) to recover 10 random but known fitness functions for HIV, each corresponding to a random fictive selective pressure on protease, and (ii) to re-estimate the fitness function that was estimated for HIV in presence of NFV, from four training data sets of varying size.

为了说明该方法的收敛性质，对该方法进行了两组模拟实验的测试：(i)恢复10个随机但已知的HIV适应性函数，每个函数都对应于蛋白酶上的一个随机的虚构选择压力，和 (ii)从四个不同大小的训练数据集中重新估计HIV在NFV存在下的适应性函数。

Each of the 10 random fitness functions used the same 114 mutations as those considered for the NFV fitness function (listed in Supplementary Material A), and was generated with 100 random interactions (BN arcs), and random fitness contributions for presence of different mutations and patterns of mutations sampled from the distribution 1 þ U(0, 1)6. For each random fitness function, a training data set of 1000 sequences was generated by evolving PI naive sequences using the fitness function for g generations, sampled from distribution P(GT). Similarly, the estimated NFV fitness function was used to create training data sets of size 513, 1026, 2052 and 10 260. The generated training sequences together with the PI naive sequences were then used to estimate the original (random or estimated NFV) fitness function.

这10个随机适应性函数每个都使用了与NFV适应性函数相同的114个突变（在补充材料A中列出），并生成了100个随机交互（BN弧线），以及不同突变和突变模式存在的随机适应性贡献，这些贡献从分布1 + U(0, 1)^6中抽样。对于每个随机适应性函数，通过使用适应性函数对PI初级序列进行g代的演化，从分布P(GT)中抽样，生成了1000个序列的训练数据集。同样，使用估计的NFV适应性函数创建了大小为513、1026、2052和10260的训练数据集。然后，将生成的训练序列和PI初级序列一起使用，以估计原始（随机或估计的NFV）适应性函数。

A weight w was defined for each arc in a random BN as the change in likelihood of that BN in the corresponding generated data set, when removing the arc, and reflects the strength of the interaction in the generated sequences. For around 60% of arcs, a value of w51 indicated that the implied interaction was not considered during evolution. These arcs were in general not present in the learned BN. Of arcs with w41, over 60% were present in the learned BN, with higher probability for arcs with higher w. Around 12% of arcs in the learned BNs indicated the obvious antagonism between different mutations at a single position (a consequence of the chosen data representation), and other artefact arcs mostly indicated associations between polymorphisms. The accuracy of the estimated fitness functions was evaluated by correlation of known and estimated fitness for 1000 independent validation sequences that were generated in the same way as the training data. Results were similar for the random and estimated NFV fitness functions. Correlation coefficients showed considerable variation with values for R2 between 0.47 and 0.84 (see Supplementary Material E). Low correlation values were associated with a banding pattern in the scatter plots. When adding subtype as an explanatory variable to the linear regression, these banding patterns disappeared, and correlation improved with R2 ranging between 0.79 and 0.94. This indicates that the estimated fitness landscapes performed well to explain intra-subtype variation of fitness under treatment, but not inter-subtype variation, which is indeed not related to fitness changes under treatment.

在随机BN中，为每个弧定义了一个权重w，作为在相应生成的数据集中该BN的似然变化，当移除弧时，反映了在生成的序列中的交互强度。对于约60%的弧，w的值小于或等于1，表示在演化过程中没有考虑到暗示的交互作用。这些弧线通常不出现在学习到的BN中。对于w超过1的弧线，超过60%出现在学习到的BN中，对于w值更高的弧线，其出现的概率更大。学习到的BN中约12%的弧线表明了同一位置上不同突变之间的明显的对抗性（这是所选择的数据表示的结果），其他的人工弧线主要表明了多态性之间的关联。通过将已知和估计的适应性对产生的1000个独立验证序列进行相关性评价，以评估估计的适应性函数的准确性。随机和估计的NFV适应性函数的结果是相似的。相关系数显示出相当大的变化，R^2的值在0.47和0.84之间（见补充材料E）。低相关性值与散点图中的带状模式相关联。当将亚型添加为线性回归的解释变量时，这些带状模式消失了，相关性得到了改善，R^2的范围在0.79和0.94之间。这表明，估计的适应性景观能够很好地解释在治疗下的适应性的亚型内变异，但不能解释亚型间变异，这实际上与治疗下的适应性变化无关。

3.2.奈非那韦适应性函数

3.2.1 Nelfinavir Bayesian network

A BN was estimatedfrom the treated sequences to estimate epistatic fitnessinteractions between the included mutations. The networkwith highest a posteriori probability, that served as a blueprintfor the fitness function, included 271 arcs (of which 33 indicatedantagonisms between different mutations at a single position.which is an artifact caused by the boolean data representation).and the corresponding fitness interactions were included in thefitness model. The network was similar to the one describedpreviously (Deforche et al., 2006), but here we allowed for more

3.2.1奈非那韦贝叶斯网络

在处理过的序列中估计出一个BN，用于估计包括的突变之间的表现性适应性交互作用。具有最高后验概率的网络作为适应性函数的蓝图，包含了271个弧（其中33个弧表示了单个位置上不同突变之间的对抗性，这是由布尔数据表示引起的人为因素），并且相应的适应性交互作用被包含在适应性模型中。网络与之前描述的一样（Deforche等人，2006年），但在这里我们允许了更多的假定交互作用，因为没有运行bootstrap程序来减小整个模型。弧的Bootstrap支持反映了弧对数据集中采样效应的鲁棒性，并且只有在从弧的存在或缺乏中推出结论时才重要。

论文研读：估计HIV-1在药物选择压力下体验的体内适应性景观对于预测治疗期间的药物抗性进化有用。