單核苷酸多態性與疾病的相關性是遺傳關聯性研究的典型代表。相對于傳統的二分類數據，單核苷酸多態性的數據具有其自身的特點，在開展 Meta 分析時常使用 5 種遺傳模型。本文通過實例展示如何使用 R 軟件 Meta 程序包實現單核苷酸多態性數據的 Meta 分析。
引用本文： 李柄輝, 王朝陽, 翁鴻, 鄭忠立, 任學群, 曾憲濤. 應用 R 軟件 meta 程序包實現遺傳關聯性研究的 Meta 分析. 中國循證醫學雜志, 2017, 17(12): 1471-1477. doi: 10.7507/1672-2531.201704062 復制
單核苷酸多態性（single nucleotide polymorphism，SNP）作為一種常見的可遺傳變異，廣泛存在于人類基因組中。隨著基因檢測技術的發展，SNP 的篩查和檢測更加快速和便捷，對 SNP 和疾病相關性的研究越來越多。針對于同一疾病相同位點的不同原始研究可能會得出不一致性的結果，或結果較為一致但因樣本量小而導致可信度低。因此，需要基于原始研究進行 Meta 分析[1]。R 軟件作為一種免費的開源性的統計軟件，已廣泛運用于 Meta 分析中[2]。盡管已有研究[3]介紹如何使用 R 軟件 meta 程序包實現二分類數據的 Meta 分析，但與傳統的二分類數據不同，SNP 的數據具有自身特點，在開展 Meta 分析時常使用 5 種遺傳模型。本文結合實例介紹如何使用 R 軟件 meta 程序包進行 SNP 數據的 Meta 分析。
1 資料與方法
1.1 示例數據
本文以 Zeng 等[4]的研究為例，該 Meta 分析探討了白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性，共計納入了 19 個病例-對照研究。本文根據實際需要，僅提取了作者、發表時間、種族、樣本量、各基因表型頻數（表 1）。在表 1 中，CC 代表野生純合型、TT 代表突變純合型、CT 代表雜合型，總樣本量用 Total 代表。
表1
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的示例數據
1.2 數據處理
開展 Meta 分析時，需要完成 5 種遺傳模型的比較：Tvs. C、TT vs. CC、CT vs. CC、CT+TT vs. CC、TT vs. CT+CC。本文以 T vs. C 為例，進行 Meta 分析，由于納入研究均為基因型的頻數，在進行 Meta 分析前，需把等位基因 T 和 C 的頻數分別計算出來，再轉化成時間發生數與總樣本量數，即轉化為經典的二分類數據。
采用 CC1、CT1、TT1 代表病例組的基因型，采用 CC0、CT0、TT0 代表對照組基因型，則病例組 C 等位基因的頻數：C1=2×CC1+1×CT1；T 等位基因的頻數：T1=2×TT1+1×CT1；總等位基因頻數：Total 1=C1+T1。以表 1 中 Gore 1998 的數據為例進行計算，病例組 T1=2×4+1×15=23，C1=2×13+1×15=41，Total1=T1+C1=41+23=64，對照組 T0=2×3+1×16=22，C0=2×13=1×16=42，Total0=T0+C0=64。其他研究數據均以此方式計算，計算結果見表 2。表 2 中 T1 和 C1 代表病例組等位基因 T 和 C，T0 和 C0 代表對照組等位基因 T 和 C，Total 1 和Total 0 分別代表病例組和對照組等位基因 T、C 出現的頻數總和。
以此類推，另外四種模型轉化后的數據分別為：TT vs. CC+TT、CT vs. CC+CT、CT+TT vs. Total、TT vs. Total。
1.3 軟件和程序包安裝
R 軟件的安裝程序包可從官方網站（http://www.r-roject.org/）下載，使用者可根據自身情況選擇鏡像下載和軟件安裝。R 軟件安裝完成后，進入 R 軟件主界面，安裝 meta 程序包，安裝命令如下：
>install.packages（“meta”）
選擇合適的鏡像安裝（CRAN）。加載 meta 程序包，命令如下：
>library（meta）
注：R 語言中區分字母大小寫，因此在安裝、加載及命令輸入時應特別注意。
2 軟件操作過程
2.1 數據輸入
R 軟件可手工輸入數據，也可從不同的數據文件中讀取數據。R 軟件可讀取文本數據、SPSS 數據、Stata 數據、SAS 數據等[5]，常用的方法是將數據記錄在 Excel 中，保存成后綴名為“csv”的文件，然后利用 read.csv（）命令讀取。本文將數據保存在計算機 D 盤下 example.csv（文件路徑為 D:/example.csv）。
注：R 語言中采用反斜杠“/”來指明文件路徑。
讀取文件，將數據定義為 data1，其命令如下：
>data1<-read.csv（“D:/example.csv”）
>data1
讀取結果見表 2。
表2
白細胞介素-1β C511T 等位基因 T 和 C 頻數計算結果
2.2 執行 Meta 分析
R 軟件 meta 程序包采用 metabin（）命令進行二分類資料的 Meta 分析，其基本命令如下：
>data2<-metabin（event.e，n.e，event.c，n.c，sm，studlab，method，data）
上述命令中，event.e 指病例組事件發生的例數，n.e 指病例組總例數；event.c 指對照組事件發生的例數，n.c 指對照組總例數，sm 為合并效應值的類型（如 RR、OR 等），method 為合并效應值的算法（如 M-H、Peto、Inverse等），data 為指定分析的數據。
根據讀取的數據，有兩種合并效應值的命令。選擇效應量為 OR，使用 MH 模型的命令分別如下：
>data2<-metabin(event.e=T1,n.e=T1+C1,event. c=T0,n.c=T0+C0,studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR")
>data2
或
>data2<-metabin(event.e=T1,n.e=total1,event. c=T0,n.c=total0,studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR")
>data2
命令執行后，計算結果見圖 1。
圖1
白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性效應值的計算結果
表3
R 軟件 meta 程序包敏感性分析結果
從圖 1 可看出，各研究間存在的異質性較大，所以選擇隨機效應模型的結果[OR=1.029 3，95%CI（0.850 2，1.246 1），P=0.767 3]。
2.3 繪制森林圖和漏斗圖
該程序包提供了繪制森林圖的功能，運行命令如下：
>forest（data2）
結果見圖 2。
圖2
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的Meta分析
R 軟件 meta 程序包采用 funnel（）命令進行漏斗圖的繪制，運行命令如下：
>funnel（data2）
漏斗圖結果見圖 3，從漏斗圖中未觀察到明顯的不對稱。
圖3
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性發表偏倚檢驗的漏斗圖
2.4 發表偏倚檢測
在 meta 程序包中，提供了多種檢測發表偏倚的方法。本文展示常用的漏斗圖法和非參數剪補法。
常用漏斗圖法為 metabias（）命令。Metabias 基于秩相關法和線性回歸法進行漏斗圖不對稱檢驗。常見方法有：Egger’s 檢驗、Begg’s 檢驗、Harbord 檢驗、Peter 檢驗、Rücker 檢驗等。對于連續型資料，一般選用 Egger’s 檢驗；對于二分類資料，如 RR、OR 等，當異質性較小的時候一般選用 Harbord 檢驗、Peter 檢驗；當異質性較大的時候，一般選用反正弦變化后的檢驗（AS-Thompson 檢驗）[6]。本例中異質性較大，因此選用 AS-Thompson 檢驗。在計算反正弦變化前，需將 metabin 命令中 sm 選擇“ASD”（即 Arcsine difference）。命令如下：
>meta1.as<-metabin(event.e=T1,n.e=T1+C1, event.c=T0,n.c=T0+C0,studlab=paste(authors,year, sep=","),method="MH",data=data1,sm="ASD")
>metabias(meta1.as,method.bias="mm")
運行結果如下：
Linear regression test of funnel plot asymmetry (methods of moment)
data: data2
t=0.816 98, df=17, p-value=0.425 2
alternative hypothesis: asymmetry in funnel plot
sample estimates:
bias se.bias slope
1.18001002 1.44435532 -0.05323743
非參數剪補法檢驗發表偏倚時采用 meta 程序包中 trimfill（）命令，命令如下：
>tf1<-trimfill(data2)
>summary(tf1,comb.random=T)
結果如下：
Number of studies combined: k=24
Random effects model: OR=0.886 3, 95%CI(0.726 6, 1.081 0), P=0.233 4
Quantifying heterogeneity: χ2=0.191 9, H=2.47, I2=83.7%
Test of heterogeneity: Q=140.75, d.f.=23, P<0.000 1
繼續使用下述命令繪制剪補后的漏斗圖（圖 4）：
>funnel(tf1,pch=ifelse(tf1$trimfill,1,23),level=0.95,comb.random=T)
圖4
白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性發表偏倚檢驗剪補后漏斗圖* *圓圈為剪補的研究
漏斗圖不對稱檢驗結果提示 P>0.05，非參數剪補法結果顯示：剪補后與剪補前無明顯改變，因此，該項指標結果不存在明顯的發表偏倚。
2.5 亞組分析
根據人種（ethnicity）進行亞組分析，命令如下：
>data3<-metabin(event.e=T1,n.e=T1+T0,event. c=C1,n.c=C1+C0,measure="OR",studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR",byvar=ethnicity)
>(data3)
結果見圖 5。
圖5
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的Meta分析的亞組分析
亞組分析結果顯示，該基因多態性和慢性牙周病發病率二者之間的相關性不受種族因素影響。
2.6 敏感性分析
在 meta 程序包中，采用 metainf（）命令通過依次剔除單項納入研究的方法進行敏感性分析。命令如下：
>metainf(data2,pooled="random")
結果見表 3。
還可繪制敏感性分析的森林圖，命令如下：
>forest(metainf(data2,pooled="random"),comb.random=TRUE)
結果見圖 6，表明本 Meta 分析結果穩定性良好。
圖6
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的 Meta 分析的敏感性分析
3 結語
R 軟件具有多個實現 Meta 分析的程序包，但針對的數據類型各異[7-14]，其中最為經典的程序包為 meta 程序包。在使用 meta 程序包進行 SNP 數據的 Meta 分析時，要注意結合 SNP 數據的特點：① 有些 SNP 的位點可能會出現三種等位基因類型，如 ABO 血型的等位基因；② SNP 數據 Meta 分析一般有 5 種基因模型進行比較，其基本原理仍為兩兩比較，與二分類數據相同；③ 在分析等位基因模型時，需先計算每個等位基因出現的頻數，然后按照傳統二分類數據進行 Meta 分析；④ 在使用 meta 程序包進行分析時，要注意參數代表的意義，以免因參數設置差錯導致結果錯誤。
從上述介紹來看，meta 程序包在進行 SNP 數據 Meta 分析時，操作簡單，用法靈活，更易接受和學習；但也存在一定的局限性[2，15]：繪制的森林圖不清晰，不能進行多變量的混合效應模型的運算，主要采用 DL 法對異質性進行估計，只能對一個分類變量進行回歸分析，不能進行假設檢驗等。但對于初學者而言，該程序包可作為 R 軟件進行 SNP 數據 Meta 分析的入門程序包。
單核苷酸多態性（single nucleotide polymorphism，SNP）作為一種常見的可遺傳變異，廣泛存在于人類基因組中。隨著基因檢測技術的發展，SNP 的篩查和檢測更加快速和便捷，對 SNP 和疾病相關性的研究越來越多。針對于同一疾病相同位點的不同原始研究可能會得出不一致性的結果，或結果較為一致但因樣本量小而導致可信度低。因此，需要基于原始研究進行 Meta 分析[1]。R 軟件作為一種免費的開源性的統計軟件，已廣泛運用于 Meta 分析中[2]。盡管已有研究[3]介紹如何使用 R 軟件 meta 程序包實現二分類數據的 Meta 分析，但與傳統的二分類數據不同，SNP 的數據具有自身特點，在開展 Meta 分析時常使用 5 種遺傳模型。本文結合實例介紹如何使用 R 軟件 meta 程序包進行 SNP 數據的 Meta 分析。
1 資料與方法
1.1 示例數據
本文以 Zeng 等[4]的研究為例，該 Meta 分析探討了白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性，共計納入了 19 個病例-對照研究。本文根據實際需要，僅提取了作者、發表時間、種族、樣本量、各基因表型頻數（表 1）。在表 1 中，CC 代表野生純合型、TT 代表突變純合型、CT 代表雜合型，總樣本量用 Total 代表。
表1
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的示例數據
1.2 數據處理
開展 Meta 分析時，需要完成 5 種遺傳模型的比較：Tvs. C、TT vs. CC、CT vs. CC、CT+TT vs. CC、TT vs. CT+CC。本文以 T vs. C 為例，進行 Meta 分析，由于納入研究均為基因型的頻數，在進行 Meta 分析前，需把等位基因 T 和 C 的頻數分別計算出來，再轉化成時間發生數與總樣本量數，即轉化為經典的二分類數據。
采用 CC1、CT1、TT1 代表病例組的基因型，采用 CC0、CT0、TT0 代表對照組基因型，則病例組 C 等位基因的頻數：C1=2×CC1+1×CT1；T 等位基因的頻數：T1=2×TT1+1×CT1；總等位基因頻數：Total 1=C1+T1。以表 1 中 Gore 1998 的數據為例進行計算，病例組 T1=2×4+1×15=23，C1=2×13+1×15=41，Total1=T1+C1=41+23=64，對照組 T0=2×3+1×16=22，C0=2×13=1×16=42，Total0=T0+C0=64。其他研究數據均以此方式計算，計算結果見表 2。表 2 中 T1 和 C1 代表病例組等位基因 T 和 C，T0 和 C0 代表對照組等位基因 T 和 C，Total 1 和Total 0 分別代表病例組和對照組等位基因 T、C 出現的頻數總和。
以此類推，另外四種模型轉化后的數據分別為：TT vs. CC+TT、CT vs. CC+CT、CT+TT vs. Total、TT vs. Total。
1.3 軟件和程序包安裝
R 軟件的安裝程序包可從官方網站（http://www.r-roject.org/）下載，使用者可根據自身情況選擇鏡像下載和軟件安裝。R 軟件安裝完成后，進入 R 軟件主界面，安裝 meta 程序包，安裝命令如下：
>install.packages（“meta”）
選擇合適的鏡像安裝（CRAN）。加載 meta 程序包，命令如下：
>library（meta）
注：R 語言中區分字母大小寫，因此在安裝、加載及命令輸入時應特別注意。
2 軟件操作過程
2.1 數據輸入
R 軟件可手工輸入數據，也可從不同的數據文件中讀取數據。R 軟件可讀取文本數據、SPSS 數據、Stata 數據、SAS 數據等[5]，常用的方法是將數據記錄在 Excel 中，保存成后綴名為“csv”的文件，然后利用 read.csv（）命令讀取。本文將數據保存在計算機 D 盤下 example.csv（文件路徑為 D:/example.csv）。
注：R 語言中采用反斜杠“/”來指明文件路徑。
讀取文件，將數據定義為 data1，其命令如下：
>data1<-read.csv（“D:/example.csv”）
>data1
讀取結果見表 2。
表2
白細胞介素-1β C511T 等位基因 T 和 C 頻數計算結果
2.2 執行 Meta 分析
R 軟件 meta 程序包采用 metabin（）命令進行二分類資料的 Meta 分析，其基本命令如下：
>data2<-metabin（event.e，n.e，event.c，n.c，sm，studlab，method，data）
上述命令中，event.e 指病例組事件發生的例數，n.e 指病例組總例數；event.c 指對照組事件發生的例數，n.c 指對照組總例數，sm 為合并效應值的類型（如 RR、OR 等），method 為合并效應值的算法（如 M-H、Peto、Inverse等），data 為指定分析的數據。
根據讀取的數據，有兩種合并效應值的命令。選擇效應量為 OR，使用 MH 模型的命令分別如下：
>data2<-metabin(event.e=T1,n.e=T1+C1,event. c=T0,n.c=T0+C0,studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR")
>data2
或
>data2<-metabin(event.e=T1,n.e=total1,event. c=T0,n.c=total0,studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR")
>data2
命令執行后，計算結果見圖 1。
圖1
白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性效應值的計算結果
表3
R 軟件 meta 程序包敏感性分析結果
從圖 1 可看出，各研究間存在的異質性較大，所以選擇隨機效應模型的結果[OR=1.029 3，95%CI（0.850 2，1.246 1），P=0.767 3]。
2.3 繪制森林圖和漏斗圖
該程序包提供了繪制森林圖的功能，運行命令如下：
>forest（data2）
結果見圖 2。
圖2
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的Meta分析
R 軟件 meta 程序包采用 funnel（）命令進行漏斗圖的繪制，運行命令如下：
>funnel（data2）
漏斗圖結果見圖 3，從漏斗圖中未觀察到明顯的不對稱。
圖3
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性發表偏倚檢驗的漏斗圖
2.4 發表偏倚檢測
在 meta 程序包中，提供了多種檢測發表偏倚的方法。本文展示常用的漏斗圖法和非參數剪補法。
常用漏斗圖法為 metabias（）命令。Metabias 基于秩相關法和線性回歸法進行漏斗圖不對稱檢驗。常見方法有：Egger’s 檢驗、Begg’s 檢驗、Harbord 檢驗、Peter 檢驗、Rücker 檢驗等。對于連續型資料，一般選用 Egger’s 檢驗；對于二分類資料，如 RR、OR 等，當異質性較小的時候一般選用 Harbord 檢驗、Peter 檢驗；當異質性較大的時候，一般選用反正弦變化后的檢驗（AS-Thompson 檢驗）[6]。本例中異質性較大，因此選用 AS-Thompson 檢驗。在計算反正弦變化前，需將 metabin 命令中 sm 選擇“ASD”（即 Arcsine difference）。命令如下：
>meta1.as<-metabin(event.e=T1,n.e=T1+C1, event.c=T0,n.c=T0+C0,studlab=paste(authors,year, sep=","),method="MH",data=data1,sm="ASD")
>metabias(meta1.as,method.bias="mm")
運行結果如下：
Linear regression test of funnel plot asymmetry (methods of moment)
data: data2
t=0.816 98, df=17, p-value=0.425 2
alternative hypothesis: asymmetry in funnel plot
sample estimates:
bias se.bias slope
1.18001002 1.44435532 -0.05323743
非參數剪補法檢驗發表偏倚時采用 meta 程序包中 trimfill（）命令，命令如下：
>tf1<-trimfill(data2)
>summary(tf1,comb.random=T)
結果如下：
Number of studies combined: k=24
Random effects model: OR=0.886 3, 95%CI(0.726 6, 1.081 0), P=0.233 4
Quantifying heterogeneity: χ2=0.191 9, H=2.47, I2=83.7%
Test of heterogeneity: Q=140.75, d.f.=23, P<0.000 1
繼續使用下述命令繪制剪補后的漏斗圖（圖 4）：
>funnel(tf1,pch=ifelse(tf1$trimfill,1,23),level=0.95,comb.random=T)
圖4
白細胞介素-1β C511T 多態性與慢性牙周病發病率的相關性發表偏倚檢驗剪補后漏斗圖* *圓圈為剪補的研究
漏斗圖不對稱檢驗結果提示 P>0.05，非參數剪補法結果顯示：剪補后與剪補前無明顯改變，因此，該項指標結果不存在明顯的發表偏倚。
2.5 亞組分析
根據人種（ethnicity）進行亞組分析，命令如下：
>data3<-metabin(event.e=T1,n.e=T1+T0,event. c=C1,n.c=C1+C0,measure="OR",studlab=paste(authors,year,sep=","),method="MH",data=data1,sm="OR",byvar=ethnicity)
>(data3)
結果見圖 5。
圖5
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的Meta分析的亞組分析
亞組分析結果顯示，該基因多態性和慢性牙周病發病率二者之間的相關性不受種族因素影響。
2.6 敏感性分析
在 meta 程序包中，采用 metainf（）命令通過依次剔除單項納入研究的方法進行敏感性分析。命令如下：
>metainf(data2,pooled="random")
結果見表 3。
還可繪制敏感性分析的森林圖，命令如下：
>forest(metainf(data2,pooled="random"),comb.random=TRUE)
結果見圖 6，表明本 Meta 分析結果穩定性良好。
圖6
白細胞介素-1β C511T 多態性與慢性牙周病發病率相關性的 Meta 分析的敏感性分析
3 結語
R 軟件具有多個實現 Meta 分析的程序包，但針對的數據類型各異[7-14]，其中最為經典的程序包為 meta 程序包。在使用 meta 程序包進行 SNP 數據的 Meta 分析時，要注意結合 SNP 數據的特點：① 有些 SNP 的位點可能會出現三種等位基因類型，如 ABO 血型的等位基因；② SNP 數據 Meta 分析一般有 5 種基因模型進行比較，其基本原理仍為兩兩比較，與二分類數據相同；③ 在分析等位基因模型時，需先計算每個等位基因出現的頻數，然后按照傳統二分類數據進行 Meta 分析；④ 在使用 meta 程序包進行分析時，要注意參數代表的意義，以免因參數設置差錯導致結果錯誤。
從上述介紹來看，meta 程序包在進行 SNP 數據 Meta 分析時，操作簡單，用法靈活，更易接受和學習；但也存在一定的局限性[2，15]：繪制的森林圖不清晰，不能進行多變量的混合效應模型的運算，主要采用 DL 法對異質性進行估計，只能對一個分類變量進行回歸分析，不能進行假設檢驗等。但對于初學者而言，該程序包可作為 R 軟件進行 SNP 數據 Meta 分析的入門程序包。
