北京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 病毒傳播方面繼續發揮著橋梁作用.