張漪麗，陳 昕

（湖南師范大學醫(yī)學院，長沙 410013）

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

1 材料與方法

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

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

2 結果

2.1 數據集建立 根據內置數據集建立本研究的演示集，建立集的代碼如下：

#建立數據集

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 個SNPs 位點中僅有一個水平的位點剔除，之后將符合納入標準的位點通過“SNPassoc”包建立調整Logistic 回歸模型，將SNPs 變換為最優(yōu)模型納入后續(xù)研究。最后有21 個位點納入后續(xù)研究。圖1 展示了數據集中所有SNPs 位點的顯著性。表1 展示了SNPs 位點模型的選擇結果和AIC 檢驗值。

2.3 SNPs-環(huán)境危險因素交互作用檢驗 建立條件Logistic 回歸模型和似然比檢驗對SNP-環(huán)境危險因素進行交互作用檢驗。經檢驗，共有1 個SNP 位點snp10008 顯著，表2 展示了顯著位點信息和似然比檢驗結果?？芍猻np10008（隱性模型 OR=1.51，95%CI=0.30-7.56）和危險因素risk 在疾病中存在交互作用（似然χ2=4.11，P=0.043）。

2.4 計算交互作用指標 表3 展示了snp10008 位點的Logistic 檢驗結果，使用“epiR”包進行交互作用計算。RERI=-0.0036（95%CI=-0.236-0.2358），S=0.244（95%CI=0.132-0.541）。交互作用的解釋受限。

3 討論

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

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

表1 候選SNPs位點信息及模型選擇

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

圖1 候選SNPs位點模型顯著性分析

表2 Snp10008位點的模型特征與似然比檢驗結果

表3 Logistic回歸分析結果

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

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