GWASSimMinimum

SNPデータの読み込み、GWA解析にgastonを用いる。Rのパッケージインストーラで簡単にダウンロード・インストールできる。

require(gaston) # for reading a vcf file and performing GWAS
require(mvtnorm) # 多変量正規分布に従う乱数の発生

シミュレーションの一連の解析を行う関数を作る（GWASSimのパーツの寄せ集め）

gwas.sim <- function(bm.sim, grm, eigen.grm, settings) {
  w.qtn <- settings$w.qtn
  n.qtn <- length(w.qtn)
  g.var.ratio <- settings$g.var.ratio
  he <- settings$he
  fdr.thresh <- settings$fdr.thresh

  n <- ncol(bm.sim)
  qtn.id <- sample(n, n.qtn) # n.qtn個のSNPをQTNとして無作為に選ぶ
  X.qtn <- as.matrix(bm.sim[, qtn.id]) # QTNの遺伝子型スコアを抜き出す
  g.qtn <- as.vector(X.qtn %*% w.qtn) # QTNによる変異を全系統について計算
  var_g.poly <- var(g.qtn)  # ここでは、50%に設定
  g.poly <-
    as.vector(rmvnorm(1, sigma = grm * var(g.qtn) * g.var.ratio))
  g <- g.qtn + g.poly # 遺伝子型値 = QTNの効果 + ポリジーン効果
  var.e <- (1 - he) / he * var(g) # 環境分散を計算する
  e <- rnorm(length(g), sd = sqrt(var.e)) # var.eを分散にもつ環境変異を生成
  y <- g + e # 遺伝変異に環境変異を加えて表現型変異を計算する
  gwa.sim <- association.test(
    bm.sim,
    Y = y,
    method = "lmm",
    response = "quantitative",
    test = "wald",
    eigenK = eigen.grm,
    p = 4
  )
  fdr <- p.adjust(gwa.sim$p, method = "BH")
  n.detected <-
    sum(bm.sim@snps$id[qtn.id] %in% gwa.sim$id[fdr < fdr.thresh])
  fnr.obs <- 1 - n.detected / n.qtn #　1 - 全QTNのうちの検出されたQTNの割合

  is.true <- rep(F, nrow(bm.sim@snps)) # QTNおよびその周辺であれば真になる
  range <- 200000 # QTNから200kb以内は、正しく検出したものとする
  for (i in 1:n.qtn) {
    qtn <- bm.sim@snps[qtn.id[i], ] # 1つのQTNずつ調べる
    cond.chr <- bm.sim@snps$chr == qtn$chr # 染色体番号が一致するか?
    cond.pos1 <- bm.sim@snps$pos > qtn$pos - range # QTNの位置より-range〜
    cond.pos2 <- bm.sim@snps$pos < qtn$pos + range # +rangeに入っているか?
    is.true <-
      is.true | (cond.chr & cond.pos1 & cond.pos2) # QTNおよびその周辺
  }
  fdr.obs <-
    sum(!(is.true) & fdr < fdr.thresh) / sum(fdr < fdr.thresh)
  fpr.obs <-
    sum(!(is.true) & fdr < fdr.thresh) / sum(!(is.true))
  c(fnr.obs, fdr.obs, fpr.obs)
}

データを読み込み、一部を利用。

# vcfファイルの読み込み
bm <- read.vcf("data/HDRA-G6-4-RDP1-RDP2-NIAS.AGCT_MAF005MIS001NR_imputed.vcf.gz") 

## ped stats and snps stats have been set. 
## 'p' has been set. 
## 'mu' and 'sigma' have been set.

bm # データのsummaryが表示される

## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are set

解析のためのデータの準備

id.med <- read.csv("data/id_med.csv", row.names = 1)$id.med  # 代表500系統のid
bm.sim <- bm[id.med, ] # idのある系統を抽出
standardize(bm.sim) <- "mu_sigma"
grm <- GRM(bm.sim)
grm <- grm + diag(0.0001, nrow(grm))
eigen.grm <- eigen(grm)

以下が、シミュレーションの本体となる。settingsの数値を変えて実行することで、様々な結果が得られる。

settings <- list(
  w.qtn = rep(1.0, 5), # 主働遺伝子の効果と数。ここでは、5QTNの効果を全て1とする。
  g.var.ratio = 1.0, #ポリジーン遺伝分散対QTN遺伝分散比。ここでは100%とする。
  he = 0.7, # 遺伝率。ここでは、表現型分散の70%が遺伝分散とする。
  fdr.thresh = 0.1 # 偽発見率のしきい値。ここでは、10%とする。
)

n.sim <- 20 # シミュレーションの回数

# 結果を入れるための入れ物を準備
sim.res <- data.frame(fnr = rep(F, n.sim), 
                      fdr = rep(F, n.sim), 
                      fpr = rep(F, n.sim))

# シミュレーションを実行
for(i in 1:n.sim) { # 1〜n.simまでiを変化させて、ループを回す
  print(i) # i回目であることを表示
  # 上で定義した関数を使って計算。結果をsim.resに代入
  sim.res[i, ] <- gwas.sim(bm.sim, grm, eigen.grm, settings) 
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20

# シミュレーション回数分のFNR, FDR, FPR
sim.res

##    fnr        fdr          fpr
## 1  0.4 0.38709677 3.121748e-04
## 2  0.4 0.28571429 3.635703e-04
## 3  0.2 0.52830189 7.265556e-04
## 4  0.6 0.25000000 2.595178e-05
## 5  0.0 0.30769231 3.114618e-04
## 6  0.4 0.68000000 4.409400e-04
## 7  0.4 0.29032258 2.336570e-04
## 8  0.2 0.15789474 1.556743e-04
## 9  0.0 0.07692308 2.597605e-05
## 10 0.2 0.25000000 2.337237e-04
## 11 0.0 0.22222222 1.038152e-04
## 12 0.0 0.30000000 1.559454e-04
## 13 0.0 0.49333333 9.597427e-04
## 14 0.2 0.12903226 1.037829e-04
## 15 0.4 0.00000000 0.000000e+00
## 16 0.2 0.44444444 2.076196e-04
## 17 0.0 0.47887324 8.834611e-04
## 18 0.0 0.22727273 1.299072e-04
## 19 0.2 0.04545455 2.599293e-05
## 20 0.0 0.54285714 4.926875e-04

# そののsummary
summary(sim.res) 

##       fnr            fdr              fpr           
##  Min.   :0.00   Min.   :0.0000   Min.   :0.0000000  
##  1st Qu.:0.00   1st Qu.:0.2061   1st Qu.:0.0001038  
##  Median :0.20   Median :0.2880   Median :0.0002206  
##  Mean   :0.19   Mean   :0.3049   Mean   :0.0002946  
##  3rd Qu.:0.40   3rd Qu.:0.4531   3rd Qu.:0.0003829  
##  Max.   :0.60   Max.   :0.6800   Max.   :0.0009597

# 箱ひげ図を描く
title <- paste(sapply(settings, paste, collapse=":"), collapse=" ")
boxplot(sim.res, ylim = c(0,1), main = title) # 箱ひげ図