疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R package 'SKAT'を用いて、SKATを実行する

今回はSNP-set Kernel Association Test(SKAT)と呼ばれるrare-variantの解析手法についてRのパッケージ'SKAT'を用いて、解説する。

追記

20160908に新しいSKATの記事を公開しています。基本的にはこのページの内容を参照していただき、さらに詳しいところは下記のリンクへ
jojoshin.hatenablog.com


はじめに

SKATについての詳細は参考文献 1 2を参照していただきたいが、従来のGenome Wide Association Studyで大きな問題点とされる『多重比較』や『結果の再現性』などの問題とともに扱われる『レアバリアントの検出』にアプローチする方法として画期的な方法である。

実践

パッケージのインストールからサンプルデータを使用できる状態にする。

install.packages("SKAT")
library(SKAT)
data(SKAT.example)
names(SKAT.example)
## [1] "Z"   "X"   "y.c" "y.b"
attach(SKAT.example)
1. Association test

アウトカムが連続変数の場合
SKATの帰無仮説のモデル(SNP setの影響は0である)での残差を算出する。
Xというのは共変量であり、それだけで説明するモデルを帰無仮説下でのモデルとする(odj)。

obj <- SKAT_Null_Model(y.c ~ X, out_type = "C")
str(obj)
## List of 9
##  $ res            : Named num [1:2000] 0.545 -0.717 0.229 -0.469 0.383 ...
##   ..- attr(*, "names")= chr [1:2000] "1" "2" "3" "4" ...
##  $ X1             : num [1:2000, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2000] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:3] "(Intercept)" "X1" "X2"
##   ..- attr(*, "assign")= int [1:3] 0 1 1
##  $ res.out        : NULL
##  $ out_type       : chr "C"
##  $ n.Resampling   : num 0
##  $ type.Resampling: chr "bootstrap"
##  $ id_include     : int [1:2000] 1 2 3 4 5 6 7 8 9 10 ...
##  $ s2             : num 0.976
##  $ n.all          : int 2000
##  - attr(*, "class")= chr "SKAT_NULL_Model"
SKAT(Z, obj)$p.value
## [1] 0.002877041

アウトカムが二値変数の場合

obj <- SKAT_Null_Model(y.b ~ X, out_type = "D")
SKAT(Z, obj)$p.value
## [1] 0.1401991

SKATはアウトカムが二値変数で、サンプルサイズが小さい場合には保守的な結果になる。そのため、サンプルが200以下の場合にはデフォルトで調整をする設定になっている。
調整した場合とそうでない場合を実際に比較してみる。
調整あり(With adjustment: default)

IDX <- c(1:100, 1001:1100)
obj.s <- SKAT_Null_Model(y.b[IDX] ~ X[IDX,], out_type = "D")
## Sample size (non-missing y and X) = 200, which is < 2000. The small sample adjustment is applied!
SKAT(Z[IDX,], obj.s, kernel="linear.weighted")$p.value
## [1] 0.1328266

調整なし(Without adjustment)

obj.s <- SKAT_Null_Model(y.b[IDX] ~ X[IDX,], out_type = "D", Adjustment = FALSE)
SKAT(Z[IDX,], obj.s, kernel="linear.weighted")$p.value
## [1] 0.147093
2. Assgin weights for each SNP

rare variantsの方がより大きなeffect sizeを有することが一般的に知られており、SNPの頻度によって重み付けをすることが考案されている(SNP-weighted kernel)。

#ちなみに、MadsenとBrowningらの重み付けは、w=1/[MAF(1-MAF)]^1/2で定義されている 
SKAT(Z, obj, kernel="linear.weighted", weights.beta=c(0.5, 0.5))$p.value
## [1] 0.4931639

自分で重み付けを作成する。

MAF <- 1:1000/1000
W <- Get_Logistic_Weights_MAF(MAF, par1=0.07, par2=150)
par(mfrow=c(1,2))
#MAFとweightの関係(全体)を図示する。
plot(MAF, W, xlab="MAF", ylab="Weights", type="l")
#MAFとweightの関係(1から100すなわちMAF=0.001〜0.1までの部分)を図示する。
plot(MAF[1:100], W[1:100], xlab="MAF", ylab="Weights", type="l")

f:id:ryosukefujii0320:20160212230447p:plain


Z(遺伝子多型)のデータを用いて、重み付けをして解析する。

weights <- Get_Logistic_Weights(Z, par1=0.07, par2=150)
SKAT(Z, obj, kernel="linear.weighted", weights=weights)$p.value
## [1] 0.3293643
3. PLINKファイルからの解析

SKATではPLINKの二値データを使用することができる。PLINKのファイルを使用するために、bed/bim/famとsetid(Generate SSD SetIDで生成可能)のファイルが必要になる。ここで使用するデータはSKAT/Meta SKAT google groupのページでダウンロードできる。PLINKで扱うbed/bim/famファイルの詳細はこちらのウェブサイトを参照してほしい。

データの読み込み

#まずはじめに、SSDファイルを作成する必要がある。
File.Bed <- "./Example1.bed"
File.Bim <- "./Example1.bim"
File.Fam <- "./Example1.fam"
#SetIDのファイルは一つ目の列にSet IDを、二つ目の列にSNP IDとする。(headerはなし)
File.SetID <- "./Example1.SetID"
#SSD fileは遺伝子多型のbinary data
File.SSD <- "./Example1.SSD"
#SSD info fileはデータの一般的な情報(最初の6列)とSNPsetの情報(残りの8列)
File.Info <- "./Example1.SSD.info"
#Set IDを作る。
Generate_SSD_SetID(File.Bed, File.Bim, File.Fam, File.SetID, File.SSD, File.Info)
## 
## Check duplicated SNPs in each SNP set
## No duplicate
## 1000 Samples, 10 Sets, 984 Total SNPs
## [1] "SSD and Info files are created!"
#FAMを読み込む。
FAM <- Read_Plink_FAM(File.Fam, Is.binary=FALSE)
y<-FAM$Phenotype

#SSDファイルを使用するときははじめにopenし、終わる前に必ず閉じること。
SSD.INFO<-Open_SSD(File.SSD, File.Info)
SSD.INFO$nSample
## [1] 1000
SSD.INFO$nSets
## [1] 10

実際に解析する

obj<-SKAT_Null_Model(y ~ 1, out_type="C")
out<-SKAT.SSD.All(SSD.INFO, obj)
out
## $results
##      SetID    P.value N.Marker.All N.Marker.Test
## 1  GENE_01 0.77747880           94            94
## 2  GENE_02 0.06245208           84            84
## 3  GENE_03 0.38416582          108           108
## 4  GENE_04 0.46179268          101           101
## 5  GENE_05 0.18548863          103           103
## 6  GENE_06 0.93255760           94            94
## 7  GENE_07 0.18897220          104           104
## 8  GENE_08 0.73081683           96            96
## 9  GENE_09 0.67366458          100           100
## 10 GENE_10 0.40310682          100           100
## 
## $P.value.Resampling
## NULL
## 
## attr(,"class")
## [1] "SKAT_SSD_ALL"

PLINK形式の共変量のファイルを持っている場合、famとcovを結合した形としてファイルを整形し、その後共変量を入れた形でSKATを用いる

File.Cov <- "./Example1.Cov"
FAM_Cov <- Read_Plink_FAM_Cov(File.Fam, File.Cov, Is.binary=FALSE)
FAM_Cov[1:5,]
#X1とX2は共変量で、それ以外はfam file形式のデータである。
##
##      FID IID PID MID Sex Phenotype         X1 X2
## 1 FID454   1   0   0   1  0.679793  1.0297614  1
## 2 FID977   1   0   0   1  0.836566  0.1846235  1
## 3 FID462   1   0   0   1 -0.408388 -0.6141158  1
## 4 FID958   1   0   0   1 -0.522305 -2.0226759  0
## 5 FID668   1   0   0   1 -0.328300 -0.8213776  0
X1 <- FAM_Cov$X1
X2 <- FAM_Cov$X2
y <- FAM_Cov$Phenotype
obj <- SKAT_Null_Model(y ~ X1 + X2, out_type="C")
out <- SKAT.SSD.All(SSD.INFO, obj)
out
## $results
##      SetID    P.value N.Marker.All N.Marker.Test
## 1  GENE_01 0.77771227           94            94
## 2  GENE_02 0.06157071           84            84
## 3  GENE_03 0.39818504          108           108
## 4  GENE_04 0.46548442          101           101
## 5  GENE_05 0.18981516          103           103
## 6  GENE_06 0.94073952           94            94
## 7  GENE_07 0.18779019          104           104
## 8  GENE_08 0.74559501           96            96
## 9  GENE_09 0.66573796          100           100
## 10 GENE_10 0.40204308          100           100
## 
## $P.value.Resampling
## NULL
## 
## attr(,"class")
## [1] "SKAT_SSD_ALL"

余談

Michael Wuの研究グループは現在、SeattleにあるFred Hutchinson Cancer Research Center (FHCRC)で研究をしている。昨年、幸運なことにお会いできる機会があったのだが、非常に気さくで人格者であった。

参考文献

  1. Michael Wu, et al. (2010) Powerful SNP-set analysis for case-control genome-wide association studies. Am J Hum Genet
  2. Michael Wu, et al. (2011) Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test. Am J Hum Genet

20160212
RF