程金瑾, 金 鑫, 汪 圳, 婁 潔

(上海大學(xué)理學(xué)院, 上海 200444)

在過去的20 年里, 中國(guó)實(shí)施了多項(xiàng)有效的公共衛(wèi)生措施來預(yù)防和治療獲得性免疫缺陷綜合癥(acquired immune deficiency syndrome, AIDS, 也稱艾滋病), 有效減少了艾滋病的傳播[1-4].然而, 男男同性性行為(men who have sex with men, MSM)群體中人類免疫缺陷病毒(human immunodeficiency virus, HIV)感染者人數(shù)的快速增長(zhǎng)與艾滋病的整體流行狀況背道而馳[5-6].在報(bào)道的中國(guó)HIV/AIDS 新增病例中, MSM 人群從2005 年占比0.7%到2014 年占比25.8%, MSM 人群感染HIV 呈現(xiàn)快速增長(zhǎng)的態(tài)勢(shì)[7-8].近10 年來, 中國(guó)MSM 人群主要有2 大特點(diǎn): ①政府的“四免一關(guān)懷”政策;②年輕人對(duì)同性戀行為、MSM 的認(rèn)同度和態(tài)度發(fā)生了很大改變.在1990 年出版的《中國(guó)性文明兩萬例》一書中, 本課題組調(diào)查了人們對(duì)同性戀的認(rèn)知、態(tài)度和行為等問題的看法[9], 并分別于2000 年和2010 年對(duì)新世紀(jì)大學(xué)生進(jìn)行了2 次調(diào)查.2000 年, 本課題組選取了26 所高校的5 070 名大學(xué)生作為樣本[10];2010 年, 本課題組從46 所大學(xué)中選取了7 829 名大學(xué)生作為調(diào)查對(duì)象[11].調(diào)查結(jié)果顯示: 近10 年, 大學(xué)生對(duì)同性戀行為和同性戀人群的正面認(rèn)同和理性態(tài)度都有了顯著提高, 各種錯(cuò)誤認(rèn)知也越來越少.表1 為2000 年和2010 年大學(xué)生對(duì)同性戀行為認(rèn)知程度對(duì)比表[11].研究還發(fā)現(xiàn), 從2000 年到2010 年, 大學(xué)生對(duì)邊緣同性戀行為的態(tài)度也發(fā)生了顯著的變化.大學(xué)生群體中同性間邊緣性行為(如同性間的擁抱、親吻、撫摸行為)的發(fā)生情況也有了明顯的增加[11].

1 研究方法

傳統(tǒng)量化流行病傳播的方法, 依賴于收集到的監(jiān)測(cè)數(shù)據(jù)和確定性模型.將確定性模型和監(jiān)測(cè)數(shù)據(jù)相結(jié)合, 可以幫助研究人員估計(jì)傳染病潛伏期、感染期的分布以及有效再生數(shù)Re.此方法的大部分監(jiān)測(cè)數(shù)據(jù)都來源于案例報(bào)告的發(fā)病率和醫(yī)院的記錄.依靠這種方法獲得的數(shù)據(jù)易受人為影響, 如在公共衛(wèi)生基礎(chǔ)設(shè)施較差地區(qū)的數(shù)據(jù)會(huì)出錯(cuò)和不完整.除了監(jiān)測(cè)數(shù)據(jù), 還可以考慮另一個(gè)重要信息來源: 基因組測(cè)序數(shù)據(jù).許多病原體的復(fù)制周期短, 突變率高, 在傳播過程中會(huì)出現(xiàn)大量遺傳變異, 因此即使在較短的時(shí)間跨度內(nèi), 流行病學(xué)進(jìn)程也會(huì)在取自宿主的病毒序列遺傳結(jié)構(gòu)中留下信號(hào).而貝葉斯系統(tǒng)發(fā)育方法就是從病毒遺傳數(shù)據(jù)中推斷流行病學(xué)過程的常用方法[12-13].貝葉斯系統(tǒng)發(fā)育方法的核心算法是Metropolis-Hastings 馬爾科夫鏈蒙特卡洛方法(Markov chain Monte Carlo, MCMC)[14-15].本工作的目標(biāo)是得到給定基因序列進(jìn)化參數(shù)的后驗(yàn)分布, 將采用出生-死亡天際線(birth-death skyline, BDSKY)方法進(jìn)行推斷, 該方法允許以非參數(shù)的方式從系統(tǒng)發(fā)育樹中提取這些信息.

本工作的序列抽樣來自2014 年北京44 個(gè)16～25 歲新感染HIV 的MSM 人群, 這些HIV病毒都是屬于07-BC 亞型, 序列長(zhǎng)度為1 201 個(gè)核苷酸的pol 基因的核苷酸序列[8].本工作通過分析這些序列, 推斷出北京地區(qū)青年MSM 人群中艾滋病流行的分子系統(tǒng)動(dòng)力學(xué).

BDSKY 過程是對(duì)種群變化的一種隨機(jī)描述, 允許個(gè)體在任何時(shí)間點(diǎn)出生或死亡.BDSKY 主要包括3 個(gè)參數(shù): 傳播速度λ、移出率γ和抽樣率ρ.

用貝葉斯框架描述流行病傳播過程.其中“出生事件”對(duì)應(yīng)個(gè)體感染;“死亡事件”對(duì)應(yīng)感染個(gè)體變?yōu)榉歉腥?從感染倉(cāng)室移出,這可以是由多個(gè)原因造成,如個(gè)體死亡、成功治療或個(gè)體行為改變).令時(shí)間序列為t0

(1)定義λi >0 為每個(gè)感染個(gè)體在時(shí)間[ti?1,ti]內(nèi)的傳播率,同時(shí)定義向量λ=(λ1,λ2,··· ,λm)為整個(gè)感染過程的傳播率向量.

(2) 定義γi >0 為每個(gè)感染個(gè)體在時(shí)間[ti?1,ti]內(nèi)的死亡(移出)率, 向量μ= (γ1,γ2,··· ,γm)為整個(gè)過程的死亡率, 并且死亡率有可能會(huì)大于出生率.

(3) 定義ψi >0 為每個(gè)感染個(gè)體在時(shí)間[ti?1,ti]內(nèi)被抽樣的概率(順序抽樣), 向量ψ=(ψ1,ψ2,··· ,ψm)為整個(gè)過程的抽樣率.

這里, 允許在時(shí)刻ti進(jìn)行特殊的抽樣(同時(shí)抽樣), 即每個(gè)感染個(gè)體在時(shí)刻ti以概率ρi被抽樣(或變?yōu)榉歉腥菊?.當(dāng)這些序列數(shù)據(jù)是在同一時(shí)間點(diǎn)被收集時(shí), 該模型中的參數(shù)ρ=0 且ρ1,ρ2,··· ,ρm=0.

令η= (λ,γ,ψ)為發(fā)育樹生成模型的參數(shù).利用MCMC 方法實(shí)現(xiàn)貝葉斯系統(tǒng)發(fā)育推理,其目的是從后驗(yàn)分布f[T,η,θ|data]中抽樣得到發(fā)育樹和模型參數(shù), 其中θ為序列演化模型參數(shù), 如核苷酸替換率向量(以每個(gè)位點(diǎn)替換率為單位分支長(zhǎng)度, 將時(shí)間樹轉(zhuǎn)換為系統(tǒng)發(fā)育樹),以及傳播樹T用來描述所有抽樣序列之間的流行病學(xué)關(guān)系.該后驗(yàn)分布可由如下貝葉斯公式得到:

這里,f[data|T,θ]是系統(tǒng)發(fā)育樹的概率密度函數(shù), 可以通過Felsenstein 算法得到[16-17].關(guān)于BDSKY 模型的詳細(xì)介紹可以參考文獻(xiàn)[18].

核苷酸替換模型在推斷演化樹和理解基因序列的進(jìn)化過程中起著很重要的作用.選擇合適的核苷酸替模型不僅可以更精確地推斷物種的演化歷史, 而且還有助于更好地了解影響序列進(jìn)化的動(dòng)力和機(jī)制.本系統(tǒng)發(fā)育分析采用的是嚴(yán)格分子時(shí)鐘, 將核苷酸替換速率固定為0.002 55[19], 并分別考慮如下4 種不同的核苷酸替換模型: GTR+Γ+I、HKY+Γ+I、JC69+Γ+I和TN93+Γ+I模型.

表2 給出了待測(cè)參數(shù)的先驗(yàn)分布.這里的樹先驗(yàn)是采用同時(shí)抽樣的BDSKY 模型.

表2 待測(cè)參數(shù)的先驗(yàn)分布Table 2 Prior distribution for parameters to the bested

這里的樹先驗(yàn)是采用同時(shí)抽樣的BDSKY 模型.本工作將MCMC 鏈長(zhǎng)設(shè)置為2 000 萬步, 舍棄前10%的樣本, 模型參數(shù)每隔1 000 步采樣一次, 并確保每個(gè)參數(shù)的有效抽樣樣本量達(dá)到200 及以上.最后通過計(jì)算這4 個(gè)模型的赤池信息量準(zhǔn)則(Akaike information criterion,AICM)值以確定最合適的模型.

2 結(jié)果討論

4 種核苷酸替換模型下參數(shù)都達(dá)到收斂狀態(tài), 且參數(shù)的有效樣本量(effective sample size,ESS)都達(dá)到幾千, 其中最少的ESS 也是582.表3 是4 個(gè)不同模型下的參數(shù)估值表, 主要包括流行病參數(shù)和序列的最近源祖先時(shí)間(the most recent common ancestors time, tMRCA).結(jié)果顯示, 這些來自北京的青年MSM 人群所感染HIV 病毒的共同祖先位與2006 年左右(95%最高后驗(yàn)概率密度區(qū)間(the high posterior density, HPD)[2005, 2008]).移出率的均值集中在0.48 左右(95%[0.45,0.53]).

表3 利用貝葉斯推斷得到的各參數(shù)估計(jì)及其HPD 區(qū)間Table 3 Parameters estimation and HPD interval from Bayesian inference

用Akaike Information Criterion (AIC)法從這4 個(gè)模型中選擇一個(gè)最合適的模型進(jìn)行進(jìn)一步分析.其中AIC 定義為AIC=ak ?alnL,k為對(duì)應(yīng)參數(shù)個(gè)數(shù),L為似然值.AICM 采用method-of-moments 進(jìn)行估計(jì), AICM 值越低的模型是越合適的.表4 為4 種不同模型之間的AICM 估值, 其中正值表示該行所在的模型擬合度比列所在的模型更好.

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

從表中可以看出, GTR+Γ+I模型是這4 個(gè)模型中最合適的, 而模型JC69+Γ+I則是其中最差的.

近10 年, 中國(guó)MSM 人群有2 大特點(diǎn): 對(duì)陽性者進(jìn)行免費(fèi)治療和社會(huì)對(duì)MSM 人群接受度大幅提高.我們希望利用這些基因序列來研究這2 個(gè)特點(diǎn)對(duì)HIV 傳播動(dòng)力學(xué)的影響.首先, 定義有效再生數(shù)Rie=λi/(γi+ψi),i= 1,2,··· ,n, 其中i為速率隨時(shí)間改變的次數(shù),Re的維數(shù)設(shè)為10.圖1 顯示了有效再生數(shù)Re的估計(jì)中值及其95%HPD(最高后驗(yàn)密度)區(qū)間.從圖中可以看出, 在最近源祖先(MRCA)處,R1e的中值為1.32(95%HPD 區(qū)間: [0.83, 2.50]), 而且在整個(gè)的傳播過程中Re都以極大的概率大于1.從2006 年開始,Re呈現(xiàn)迅速增大的態(tài)勢(shì), 并最終趨于某穩(wěn)態(tài), 這預(yù)示HIV 在該青年人群中的傳播增勢(shì)迅猛.這個(gè)結(jié)果也與已有報(bào)道[7-8]相一致.雖然中國(guó)在2003 年已經(jīng)向HIV/AIS 患者提供免費(fèi)治療, 但近10 年來北京年輕人群對(duì)MSM 人群及行為接受度大幅度提高可能是該傳播增勢(shì)迅猛的重要原因.

圖1 GTR+Γ +I 模型下的Re 估計(jì)中值及其95%HPDFig.1 Median estimates and 95% HPD intervals for Re of GTR+Γ +I

3 結(jié)束語

基因序列中包含了病毒遺傳和變異的所有信息, 甚至隱藏著傳染病的動(dòng)力學(xué)信息.從某種程度上來說, 基因序列得到的信息比從傳統(tǒng)的統(tǒng)計(jì)方法推斷的信息更可靠.本工作通過貝葉斯推斷分析44 個(gè)北京年輕男同HIV 患者的基因序列, 得到了該傳染病有價(jià)值的傳播信息,如有效再生數(shù)Re.有效再生數(shù)對(duì)傳染病分析意義重大, 可以幫助對(duì)傳染病進(jìn)行后續(xù)預(yù)測(cè)和制定防治措施.令人擔(dān)憂的是, 雖然年輕人對(duì)MSM 的接受度已大幅度提高, 但迫于社會(huì)和文化輿論壓力, 大部分中國(guó)MSM 人群會(huì)隱藏自己的性取向而正常結(jié)婚, 因此男同性戀者可能會(huì)在HIV 病毒傳播方面繼續(xù)發(fā)揮著橋梁作用.