北京MSM 人群HIV 传播的贝叶斯推断

程金瑾, 金 鑫, 汪 圳, 娄 洁

(上海大学理学院, 上海 200444)

在过去的20 年里, 中国实施了多项有效的公共卫生措施来预防和治疗获得性免疫缺陷综合症(acquired immune deficiency syndrome, AIDS, 也称艾滋病), 有效减少了艾滋病的传播[1-4].然而, 男男同性性行为(men who have sex with men, MSM)群体中人类免疫缺陷病毒(human immunodeficiency virus, HIV)感染者人数的快速增长与艾滋病的整体流行状况背道而驰[5-6].在报道的中国HIV/AIDS 新增病例中, MSM 人群从2005 年占比0.7%到2014 年占比25.8%, MSM 人群感染HIV 呈现快速增长的态势[7-8].近10 年来, 中国MSM 人群主要有2 大特点: ①政府的“四免一关怀”政策;②年轻人对同性恋行为、MSM 的认同度和态度发生了很大改变.在1990 年出版的《中国性文明两万例》一书中, 本课题组调查了人们对同性恋的认知、态度和行为等问题的看法[9], 并分别于2000 年和2010 年对新世纪大学生进行了2 次调查.2000 年, 本课题组选取了26 所高校的5 070 名大学生作为样本[10];2010 年, 本课题组从46 所大学中选取了7 829 名大学生作为调查对象[11].调查结果显示: 近10 年, 大学生对同性恋行为和同性恋人群的正面认同和理性态度都有了显著提高, 各种错误认知也越来越少.表1 为2000 年和2010 年大学生对同性恋行为认知程度对比表[11].研究还发现, 从2000 年到2010 年, 大学生对边缘同性恋行为的态度也发生了显著的变化.大学生群体中同性间边缘性行为(如同性间的拥抱、亲吻、抚摸行为)的发生情况也有了明显的增加[11].

1 研究方法

传统量化流行病传播的方法, 依赖于收集到的监测数据和确定性模型.将确定性模型和监测数据相结合, 可以帮助研究人员估计传染病潜伏期、感染期的分布以及有效再生数Re.此方法的大部分监测数据都来源于案例报告的发病率和医院的记录.依靠这种方法获得的数据易受人为影响, 如在公共卫生基础设施较差地区的数据会出错和不完整.除了监测数据, 还可以考虑另一个重要信息来源: 基因组测序数据.许多病原体的复制周期短, 突变率高, 在传播过程中会出现大量遗传变异, 因此即使在较短的时间跨度内, 流行病学进程也会在取自宿主的病毒序列遗传结构中留下信号.而贝叶斯系统发育方法就是从病毒遗传数据中推断流行病学过程的常用方法[12-13].贝叶斯系统发育方法的核心算法是Metropolis-Hastings 马尔科夫链蒙特卡洛方法(Markov chain Monte Carlo, MCMC)[14-15].本工作的目标是得到给定基因序列进化参数的后验分布, 将采用出生-死亡天际线(birth-death skyline, BDSKY)方法进行推断, 该方法允许以非参数的方式从系统发育树中提取这些信息.

本工作的序列抽样来自2014 年北京44 个16∼25 岁新感染HIV 的MSM 人群, 这些HIV病毒都是属于07-BC 亚型, 序列长度为1 201 个核苷酸的pol 基因的核苷酸序列[8].本工作通过分析这些序列, 推断出北京地区青年MSM 人群中艾滋病流行的分子系统动力学.

BDSKY 过程是对种群变化的一种随机描述, 允许个体在任何时间点出生或死亡.BDSKY 主要包括3 个参数: 传播速度λ、移出率γ和抽样率ρ.

用贝叶斯框架描述流行病传播过程.其中“出生事件”对应个体感染;“死亡事件”对应感染个体变为非感染(从感染仓室移出,这可以是由多个原因造成,如个体死亡、成功治疗或个体行为改变).令时间序列为t0

(1)定义λi >0 为每个感染个体在时间[ti−1,ti]内的传播率,同时定义向量λ=(λ1,λ2,··· ,λm)为整个感染过程的传播率向量.

(2) 定义γi >0 为每个感染个体在时间[ti−1,ti]内的死亡(移出)率, 向量µ= (γ1,γ2,··· ,γm)为整个过程的死亡率, 并且死亡率有可能会大于出生率.

(3) 定义ψi >0 为每个感染个体在时间[ti−1,ti]内被抽样的概率(顺序抽样), 向量ψ=(ψ1,ψ2,··· ,ψm)为整个过程的抽样率.

这里, 允许在时刻ti进行特殊的抽样(同时抽样), 即每个感染个体在时刻ti以概率ρi被抽样(或变为非感染者).当这些序列数据是在同一时间点被收集时, 该模型中的参数ρ=0 且ρ1,ρ2,··· ,ρm=0.

令η= (λ,γ,ψ)为发育树生成模型的参数.利用MCMC 方法实现贝叶斯系统发育推理,其目的是从后验分布f[T,η,θ|data]中抽样得到发育树和模型参数, 其中θ为序列演化模型参数, 如核苷酸替换率向量(以每个位点替换率为单位分支长度, 将时间树转换为系统发育树),以及传播树T用来描述所有抽样序列之间的流行病学关系.该后验分布可由如下贝叶斯公式得到:

这里,f[data|T,θ]是系统发育树的概率密度函数, 可以通过Felsenstein 算法得到[16-17].关于BDSKY 模型的详细介绍可以参考文献[18].

核苷酸替换模型在推断演化树和理解基因序列的进化过程中起着很重要的作用.选择合适的核苷酸替模型不仅可以更精确地推断物种的演化历史, 而且还有助于更好地了解影响序列进化的动力和机制.本系统发育分析采用的是严格分子时钟, 将核苷酸替换速率固定为0.002 55[19], 并分别考虑如下4 种不同的核苷酸替换模型: GTR+Γ+I、HKY+Γ+I、JC69+Γ+I和TN93+Γ+I模型.

表2 给出了待测参数的先验分布.这里的树先验是采用同时抽样的BDSKY 模型.

表2 待测参数的先验分布Table 2 Prior distribution for parameters to the bested

这里的树先验是采用同时抽样的BDSKY 模型.本工作将MCMC 链长设置为2 000 万步, 舍弃前10%的样本, 模型参数每隔1 000 步采样一次, 并确保每个参数的有效抽样样本量达到200 及以上.最后通过计算这4 个模型的赤池信息量准则(Akaike information criterion,AICM)值以确定最合适的模型.

2 结果讨论

4 种核苷酸替换模型下参数都达到收敛状态, 且参数的有效样本量(effective sample size,ESS)都达到几千, 其中最少的ESS 也是582.表3 是4 个不同模型下的参数估值表, 主要包括流行病参数和序列的最近源祖先时间(the most recent common ancestors time, tMRCA).结果显示, 这些来自北京的青年MSM 人群所感染HIV 病毒的共同祖先位与2006 年左右(95%最高后验概率密度区间(the high posterior density, HPD)[2005, 2008]).移出率的均值集中在0.48 左右(95%[0.45,0.53]).

表3 利用贝叶斯推断得到的各参数估计及其HPD 区间Table 3 Parameters estimation and HPD interval from Bayesian inference

用Akaike Information Criterion (AIC)法从这4 个模型中选择一个最合适的模型进行进一步分析.其中AIC 定义为AIC=ak −alnL,k为对应参数个数,L为似然值.AICM 采用method-of-moments 进行估计, AICM 值越低的模型是越合适的.表4 为4 种不同模型之间的AICM 估值, 其中正值表示该行所在的模型拟合度比列所在的模型更好.

表4 AICM 模型比较Table 4 Model compared by AICM

从表中可以看出, GTR+Γ+I模型是这4 个模型中最合适的, 而模型JC69+Γ+I则是其中最差的.

近10 年, 中国MSM 人群有2 大特点: 对阳性者进行免费治疗和社会对MSM 人群接受度大幅提高.我们希望利用这些基因序列来研究这2 个特点对HIV 传播动力学的影响.首先, 定义有效再生数Rie=λi/(γi+ψi),i= 1,2,··· ,n, 其中i为速率随时间改变的次数,Re的维数设为10.图1 显示了有效再生数Re的估计中值及其95%HPD(最高后验密度)区间.从图中可以看出, 在最近源祖先(MRCA)处,R1e的中值为1.32(95%HPD 区间: [0.83, 2.50]), 而且在整个的传播过程中Re都以极大的概率大于1.从2006 年开始,Re呈现迅速增大的态势, 并最终趋于某稳态, 这预示HIV 在该青年人群中的传播增势迅猛.这个结果也与已有报道[7-8]相一致.虽然中国在2003 年已经向HIV/AIS 患者提供免费治疗, 但近10 年来北京年轻人群对MSM 人群及行为接受度大幅度提高可能是该传播增势迅猛的重要原因.

图1 GTR+Γ +I 模型下的Re 估计中值及其95%HPDFig.1 Median estimates and 95% HPD intervals for Re of GTR+Γ +I

3 结束语

基因序列中包含了病毒遗传和变异的所有信息, 甚至隐藏着传染病的动力学信息.从某种程度上来说, 基因序列得到的信息比从传统的统计方法推断的信息更可靠.本工作通过贝叶斯推断分析44 个北京年轻男同HIV 患者的基因序列, 得到了该传染病有价值的传播信息,如有效再生数Re.有效再生数对传染病分析意义重大, 可以帮助对传染病进行后续预测和制定防治措施.令人担忧的是, 虽然年轻人对MSM 的接受度已大幅度提高, 但迫于社会和文化舆论压力, 大部分中国MSM 人群会隐藏自己的性取向而正常结婚, 因此男同性恋者可能会在HIV 病毒传播方面继续发挥着桥梁作用.