一、為什么要校正case和control數量比例不平衡情況
試問作為生信屆人員，最怕的是什么，當然是統計結果不靠譜。統計結果不靠譜包括兩方面：一個是假陰性，一個是假陽性。假陰性可以理解為白天鵝被誤當成丑小鴨了，假陽性可以理解為一大堆青蛙，你不知道哪個才是你的真命天子。假陰性就罷了，最多讓你錯過發現真理的機會，但萬一假陽性呢，你拿着一個看似完美的結果吭哧吭哧做實驗驗證，一年半載的周期下來，什么結果都驗證不出來，豈不是坑了做實驗的人。因此，我們就要在源頭上，把這個不靠譜的統計結果杜絕出去。
上一篇文章什么！GWAS研究中case和control的比例是有講究的？就講到GWAS分析中，如果case和control數量比例失衡的話，會產出非常多的假陽性結果，而用SAIGE模型做GWAS分析可以校正這種數量比例不平衡的情況。下面具體講講怎么應用SAIGE模型。
二、怎么校正：SAIGE的下載和安裝
1、下載SAIGE
此操作在Linux上進行，系統要求R-3.5.1, gcc >= 5.5.0, cmake 3.8.1
wget https://github.com/weizhouUMICH/SAIGE/blob/master/SAIGE_0.35.8.1_R_x86_64-pc-linux-gnu.tar.gz
2、安裝SAIGE
R CMD INSTALL SAIGE_XX_R_x86_64-pc-linux-gnu.tar.gz
3、安裝SAIGE所依賴的其他R包：Rcpp, RcppArmadillo, RcppParallel, data.table, SPAtest, RcppEigen, Matrix, methods, optparse
以下兩個方法二選一：
如果是用conda的話，則用以下命令：
conda install -n r-env r-Rcpp #r-env是指conda下的R環境 conda install -n r-env r-RcppArmadillo conda install -n r-env r-RcppParallel conda install -n r-env r-SPAtest conda install -n r-env r-RcppEigen conda install -n r-env r-optparse
也可以進入R，用install.packages( )安裝：
install.packages("Rcpp")
install.packages("RcppArmadillo")
install.packages("RcppParallel")
install.packages("SPAtest")
install.packages("RcppEigen")
install.packages("optparse")
三、怎么校正：SAIGE的分析、解讀
1、第一步，計算NULL GLMM
1）計算NULL GLMM的命令：
Rscript step1_fitNULLGLMM.R \ --plinkFile=./input/plinkforGRM_1000samples_10kMarkers \ --phenoFile=./input/pheno_1000samples.txt \ --phenoCol=y \ --covarColList=x1,x2 \ --sampleIDColinphenoFile=IID \ --traitType=binary \ --outputPrefix=./output/example \ --nThreads=4 \ --LOCO=TRUE
plinkFile為plink的輸入文件（bed, bim, fam格式）
phenoFile文件格式如下，第一列代表研究的表型，第二列及第N-1列代表協變量，最后一列IID為個體的ID號：
--phenoCol=y # 指定你要研究的表型列名，在本次例子中，指定y的表型分析。
--covarColList=x1,x2 #指定加入的協變量
--sampleIDColinphenoFile=IID #指定樣本的ID
--traitType=binary #指定研究的表型的類型，binary指二分類，即case和control
--outputPrefix #生成文件的輸出路徑
2）輸出文件的結果解讀：
這個步驟會生成三個文件，分別為：example.rda、example.varianceRatio.txt、example_30markers.SAIGE.results.txt
第一個文件：example.rda，是一個model file
可以用R打開：
load("./output/example.rda")
names(modglmm)
modglmm$theta
第二個文件：example_30markers.SAIGE.results.txt，隨意選取位點的關聯分析結果
第三個文件：example.varianceRatio.txt
2、計算每個marker的SPA得分
1）計算每個marker的SNP得分命令：
Rscript step2_SPAtests.R \ --dosageFile=./input/dosage_10markers.txt \ --dosageFileNrowSkip=1 \ --dosageFileNcolSkip=6 \ --dosageFilecolnamesSkip=CHR,SNP,CM,POS,EFFECT_ALLELE,ALT_ALLELE \ --minMAF=0.0001 \ --sampleFile=./input/sampleIDindosage.txt \ --GMMATmodelFile=./output/example.rda \ --varianceRatioFile=./output/example.varianceRatio.txt \ --SAIGEOutputFile=./output/example.plainDosage.SAIGE.txt \ --numLinesOutput=2 \ --IsOutputAFinCaseCtrl=TRUE
dosage_10markers.txt的文件格式如下，類似於plink的ped格式，前六列分別為：CHR, SNP, CM, POS, COUNTED, ALT, 后面為個體的ID號：
sampleIDindosage.txt文件為樣本ID名：
example.rda、example.varianceRatio.txt為第一步生成的兩個文件。
2）輸出文件的結果解讀：
生成example.plainDosage.SAIGE.txt文件，其內容如下：
其中，P值（紅框）即為我們校正case和control數量比例不平衡以后得到的GWAS結果，p.value.NA為不校正case和control數量不平衡的結果。
參數說明：
CHR: chromosome
POS: genome position
SNPID: variant ID
Allele1: Ref allele
Allele2: Alt allele
AC_Allele2: allele count of Alt allele
AF_Allele2: allele frequency of Alt allele
N: sample size
BETA: effect size
SE: standard error of BETA
Tstat: score statistic
p.value: p value with SPA applied
p.value.NA: p value when SPA is not applied
Is.SPA.converge: whether SPA is converged or not
varT: estimated variance of score statistic with sample related incorporated
varTstar: variance of score statistic without sample related incorporated
AF.Cases: allele frequency of allele 2 in cases (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)
AF.Controls: allele frequency of allele 2 in controls (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)
至此，校正GWAS分析中case和control數量比例嚴重失衡的工作就完成了。導師再也不用擔心你給出一堆假陽性結果了。