SKAT

RのSKATパッケージを使いたい。
Z：行にバリアント、列にサンプルを持つ行列とする。値は0/1/2。
X：共変量行列
まず、null modelを作り(obj)これにテストしたいバリアント行列(Z)を投げる。目的変数が連続値の場合out_type="C"を指定し、二値変数の場合out_type="B"を指定する。

obj = SKAT_Null_Model(y.c ~ X, out_type="C")
SKAT(Z, obj)$p.value

obj = SKAT_Null_Model(y.b ~ X, out_type="C")
SKAT(Z, obj)$p.value

上記のモデルはサンプル数が少ない時にやや保守的であるらしい。最近のバージョンではパーミュテーションp値を返すSKATBinaryという関数が実装されている。

out = SKATbinary(Z, obj)
SKAT(Z, obj)$p.value

SNPごとに効果量を変えるカーネルを使いたい場合、kernelで指定する。

SKAT(Z, obj, kernel = "linear.weighted")$p.value

デフォルトではbeta(maf, 1, 25)に従う重みを与えている。もし、beta(maf, 0.5, 0.5)に従う重みを与えたい場合は以下のようにして変更できる。ちなみにこれはかなりゆっくりとした勾配になる。

SKAT(Z, obj, kernel = "linear.weighted", weights.beta = c(0.5, 0.5))$p.value

SKAT-O

SKAT-Oとは、SKATとburden-testをブレンドした検定方法である。ブレンド率$\rho$はr.corrでコントロールする。
r.corr=0のとき: SKATに等しい
r.corr=1のとき; Burden testに等しい
method="SKATO"を実行すると$\rho$を色々に変えて検定を行い、一番小さなp値を返す(とマニュアルにあるが、実際実験してみるとそうならなかった。ソースコード要確認。
PLINK形式のファイルを用いて、複数のSNPセットに対してSKATを適応する。
ひつようなもの
bed/bim/fam/covとSetID
注：SetIDファイルとは

GENE_01 SNP_01
GENE_01 SNP_02
GENE_01 SNP_03
...
GENE_02 SNP_20
GENE_02 SNP_22

のように記載された、variantがどのセットに属するかのファイル。bedtoolsなどを用いて作ればよい。これらをカレントディレクトリに用意した上で以下のコマンドを実行する。

File.Bed<-"./SKATBinary.example.bed"
File.Bim<-"./SKATBinary.example.bim"
File.Fam<-"./SKATBinary.example.fam"
File.Cov<-"./SKATBinary.example.cov"
File.SetID<-"./SKATBinary.example.SetID"
File.SSD<-"./SKATBinary.example.SSD"
File.Info<-"./SKATBinary.example.SSD.info"
Generate_SSD_SetID(File.Bed, File.Bim, File.Fam, File.SetID, File.SSD, File.Info)

これにより.SSDファイルと.SSD.infoファイルが作成される。今後、この二つのファイルを用いて解析を行う。
phenotypeとcovariateが一緒になったファイルをつくっているところ

FAM<-Read_Plink_FAM_Cov(File.Fam, File.Cov, Is.binary=TRUE, cov_header=FALSE)

SSDファイルをオープンする。あとでクローズしなければならない。

SSD.INFO<-Open_SSD(File.SSD, File.Info)

null modelを作って

obj = SKAT_Null_Model(Phenotype ~ COV1 + COV2, out_type="D", data=FAM, Adjustment=FALSE)

SKAT

out.skat = SKATBinary.SSD.All(SSD.INFO, obj, method="SKAT")

SKAT-O

out.skato = SKATBinary.SSD.All(SSD.INFO, obj, method="SKATO")

多分EPACTS推し