使用R 分析單核苷酸多態(tài)性位點(diǎn)-環(huán)境交互作用 對(duì)疾病影響的方法

張漪麗，陳 昕

（湖南師范大學(xué)醫(yī)學(xué)院，長(zhǎng)沙 410013）

R 語(yǔ)言是由新西蘭奧克蘭大學(xué)的羅斯·伊哈卡和羅伯特·杰特曼開(kāi)發(fā)的一款免費(fèi)開(kāi)源統(tǒng)計(jì)軟件。相較于其它商業(yè)統(tǒng)計(jì)軟件，R 可以根據(jù)自身需求創(chuàng)建自定義函數(shù)并以包（packages）的形式分享給世界上各地的使用者，達(dá)到各種功能的即時(shí)裝載和演算。目前，R 在世界上已經(jīng)擁有了近千萬(wàn)的用戶，在流行病學(xué)、統(tǒng)計(jì)學(xué)、生命科學(xué)、金融學(xué)、農(nóng)學(xué)等領(lǐng)域廣泛應(yīng)用［1］。單核苷酸多態(tài)性（Single nucleotide polymorphisms，SNPs）指由單個(gè)核苷酸的改變而引起的DNA 序列的改變，造成包括人類在內(nèi)的物種之間染色體基因組的多樣性。自從2008 年全基因組關(guān)聯(lián)研究的方案提出至今，人類已經(jīng)建立了完備的基因數(shù)據(jù)庫(kù)，大量SNPs 位點(diǎn)被測(cè)序研究。［2］目前已經(jīng)證實(shí)，SNPs 位點(diǎn)的突變與疾病發(fā)生發(fā)展呈高度關(guān)聯(lián)，研究SNPs 位點(diǎn)對(duì)疾病的作用，繼而以此作為有效的篩檢指標(biāo)在早期預(yù)防疾病已成為流行病學(xué)研究的熱點(diǎn)。并且在生物醫(yī)學(xué)和臨床醫(yī)學(xué)領(lǐng)域，研究個(gè)案的SNPs 位點(diǎn)對(duì)個(gè)體化醫(yī)療和精準(zhǔn)醫(yī)療提供了理論支撐。交互作用是在各學(xué)科研究中常見(jiàn)的問(wèn)題。指納入研究的兩個(gè)以上的因素之間產(chǎn)生相互效應(yīng)，使因素聯(lián)合作用時(shí)產(chǎn)生的效應(yīng)不等于單獨(dú)使用時(shí)的效應(yīng)之和［3］。交互作用按照方向分為協(xié)同作用和拮抗作用，按照意義分為統(tǒng)計(jì)學(xué)交互作用、生物學(xué)交互作用和公共衛(wèi)生學(xué)交互作用。在流行病學(xué)和生物醫(yī)學(xué)的研究中，常通過(guò)統(tǒng)計(jì)學(xué)交互作用為探討生物學(xué)交互作用提供線索［4］。

1 材料與方法

1.1 材料 本研究所使用的資料為R 包“SNPassoc”中內(nèi)置數(shù)據(jù)集“SNPs”，調(diào)用ID 號(hào)，病例/對(duì)照分組，性別和35 個(gè)SNP 位點(diǎn)，共157 個(gè)案例。建立虛擬因子“risk”作為環(huán)境危險(xiǎn)因素納入數(shù)據(jù)集，“risk”的取值為二分類0，1變量，其值通過(guò)設(shè)置種子數(shù)，用sample 函數(shù)隨機(jī)產(chǎn)生，個(gè)數(shù)與案例數(shù)相同。以此演示兩因素的交互作用分析。

1.2 方法 使用自建函數(shù)將單水平的SNP 位點(diǎn)篩除。SNPs 位點(diǎn)通過(guò)“SNPassoc”包進(jìn)行模型篩選。模型建立Logistic 回歸模型，調(diào)整變量為性別。顯性模型、隱性模型、超顯性模型作為二分類變量納入研究，共顯性模型作為三分類變量納入研究。最優(yōu)模型選擇根據(jù)赤池信息標(biāo)準(zhǔn)值（Akaike information criterion，AIC）界定。分別建立兩個(gè)Logistic 回歸模型，因變量為是否患病，自變量為risk 和SNPs 候選位點(diǎn)，第一個(gè)模型僅考慮risk 和SNPs的獨(dú)立作用，第二個(gè)模型在此基礎(chǔ)上加入兩個(gè)因素的交互項(xiàng)，全部建立新的啞變量后再納入研究，其中risk 0 值為未暴露，1 為有暴露。似然比檢驗(yàn)被用于兩個(gè)Logistic回歸模型的統(tǒng)計(jì)學(xué)分析，似然比計(jì)算為卡方檢驗(yàn)法，當(dāng)P 值顯著說(shuō)明交互項(xiàng)在模型間有統(tǒng)計(jì)學(xué)意義。使用R“epiR”包對(duì)顯著交互作用的SNP 位點(diǎn)進(jìn)行定量分析。使用超頻相對(duì)危險(xiǎn)度（relative excess risk of interaction，RERI）和協(xié)同指數(shù)（synergy index，S）指標(biāo)進(jìn)行評(píng)價(jià)。本研究中所有分析P<0.05 被認(rèn)為顯著。所有統(tǒng)計(jì)學(xué)分析均在R 軟件（https：//www.r-project.org/，版本：3.6.2）上完成。

2 結(jié)果

2.1 數(shù)據(jù)集建立 根據(jù)內(nèi)置數(shù)據(jù)集建立本研究的演示集，建立集的代碼如下：

#建立數(shù)據(jù)集

library（"SNPassoc"）

data（SNPs）

SNPs<-SNPs［，-c（4，5）］

set.seed（2020）

risk<-sample（0：1，157，replace=T）

SNPs<-cbind（risk，SNPs）

2.2 SNPs 篩選和模型選擇 將候選35 個(gè)SNPs 位點(diǎn)中僅有一個(gè)水平的位點(diǎn)剔除，之后將符合納入標(biāo)準(zhǔn)的位點(diǎn)通過(guò)“SNPassoc”包建立調(diào)整Logistic 回歸模型，將SNPs 變換為最優(yōu)模型納入后續(xù)研究。最后有21 個(gè)位點(diǎn)納入后續(xù)研究。圖1 展示了數(shù)據(jù)集中所有SNPs 位點(diǎn)的顯著性。表1 展示了SNPs 位點(diǎn)模型的選擇結(jié)果和AIC 檢驗(yàn)值。

2.3 SNPs-環(huán)境危險(xiǎn)因素交互作用檢驗(yàn) 建立條件Logistic 回歸模型和似然比檢驗(yàn)對(duì)SNP-環(huán)境危險(xiǎn)因素進(jìn)行交互作用檢驗(yàn)。經(jīng)檢驗(yàn)，共有1 個(gè)SNP 位點(diǎn)snp10008 顯著，表2 展示了顯著位點(diǎn)信息和似然比檢驗(yàn)結(jié)果。可知snp10008（隱性模型 OR=1.51，95%CI=0.30-7.56）和危險(xiǎn)因素risk 在疾病中存在交互作用（似然χ2=4.11，P=0.043）。

2.4 計(jì)算交互作用指標(biāo) 表3 展示了snp10008 位點(diǎn)的Logistic 檢驗(yàn)結(jié)果，使用“epiR”包進(jìn)行交互作用計(jì)算。RERI=-0.0036（95%CI=-0.236-0.2358），S=0.244（95%CI=0.132-0.541）。交互作用的解釋受限。

3 討論

本文通過(guò)分析程序包內(nèi)置的SNPs 位點(diǎn)的模擬數(shù)據(jù)集，演示了如何對(duì)SNPs-環(huán)境的交互作用進(jìn)行統(tǒng)計(jì)學(xué)分析和解釋。由于資料并不是真實(shí)生物醫(yī)學(xué)研究中的SNP，部分位點(diǎn)有邏輯錯(cuò)誤，并沒(méi)有研究的意義。并且由于157 的樣本量相較于真實(shí)研究中過(guò)少，部分位點(diǎn)呈單水平或者不通過(guò)哈溫平衡檢驗(yàn)，損失了部分信息。

本研究中評(píng)價(jià)交互作用的方法為似然比檢驗(yàn)法，在SNP 相關(guān)研究中，相較于Wald 系數(shù)檢驗(yàn)法，似然比檢驗(yàn)由于其特征，規(guī)避了可能的混雜效應(yīng)，更加精確可靠，受到廣泛的認(rèn)可［5］。

表1 候選SNPs位點(diǎn)信息及模型選擇

根據(jù)國(guó)內(nèi)外學(xué)者的研究報(bào)告，在生物學(xué)交互作用的研究中，相加模型的統(tǒng)計(jì)指標(biāo)相較于相乘模型更能反映真實(shí)的效應(yīng)［6］。然而，本研究選擇Logistic 回歸法，本質(zhì)仍是相乘模型。原因是目前在研究SNPs 交互作用中，Logistic 回歸仍為最廣泛而有效的手段。為了探索生物學(xué)意義，我們?cè)谙喑四Ｐ偷慕Y(jié)果上納入了RERI 和S 等相加模型指標(biāo)進(jìn)行分析，以期綜合評(píng)價(jià)。值得注意的是，相乘模型指標(biāo)的比值大于1 時(shí)，相加模型檢驗(yàn)指標(biāo)也有協(xié)同的顯著性意義。但是相乘模型小于1 時(shí)，相加模型檢驗(yàn)指標(biāo)會(huì)出現(xiàn)協(xié)同、拮抗、無(wú)統(tǒng)計(jì)學(xué)意義三種可能情況［7］。其原因是統(tǒng)計(jì)學(xué)的交互作用與生物學(xué)的交互作用是復(fù)雜的，并沒(méi)有絕對(duì)的對(duì)應(yīng)關(guān)系。因此，在遇到這種情況時(shí)，應(yīng)該結(jié)合自身的專業(yè)知識(shí)和既往研究報(bào)告，做出合理的判斷。比如在SNP 研究中可以查詢?cè)撐稽c(diǎn)的文獻(xiàn)報(bào)告情況和位置基因，判斷基因是否有與待研究因素的特征。

圖1 候選SNPs位點(diǎn)模型顯著性分析

表2 Snp10008位點(diǎn)的模型特征與似然比檢驗(yàn)結(jié)果

表3 Logistic回歸分析結(jié)果

本研究中的顯著位點(diǎn)snp10008 由于樣本量有限導(dǎo)致其在計(jì)算參數(shù)時(shí)標(biāo)準(zhǔn)誤過(guò)大，在研究中通常要考慮采用其他指標(biāo)評(píng)估是否SNP 位點(diǎn)真正顯著，而不能僅考慮logistic 回歸值。由于數(shù)據(jù)不夠理想，未顯示出交互作用。 另外，本文使用的是參數(shù)檢驗(yàn)法，研究SNPs 交互作用還有一種非參數(shù)檢驗(yàn)法，廣義多因子降維法，該方法采用數(shù)學(xué)模型對(duì)研究因素進(jìn)行降維，彌補(bǔ)了Logistic回歸在處理高階交互作用時(shí)的局限性［8］。在常見(jiàn)的慢性病如高血壓糖尿病中已有廣泛而成功應(yīng)用。但是值得注意的是，GMDR 對(duì)數(shù)據(jù)要求較高，并且采用無(wú)模型方法，不可完全替代傳統(tǒng)回歸分析法。

隨著計(jì)算機(jī)技術(shù)和人工智能的發(fā)展，高通量基因組學(xué)的研究成為現(xiàn)實(shí)，SNPs 位點(diǎn)已經(jīng)成為生物醫(yī)學(xué)研究的熱點(diǎn)。而R 語(yǔ)言擁有良好的生態(tài)系統(tǒng)，軟件包多為專業(yè)的統(tǒng)計(jì)學(xué)家或計(jì)算機(jī)專業(yè)人員制作，并在公布前進(jìn)行內(nèi)部同行評(píng)審，保證了該統(tǒng)計(jì)軟件的精確性和可靠性。本文旨在提供新的研究方法，其中使用到的所有R 包已在基因組學(xué)研究中廣泛應(yīng)用，為流行病學(xué)研究人員帶來(lái)了便利。