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
outfile <- settings$outfile
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")
# マンハッタンプロットを描く
if(!is.null(outfile)) {
col <- rep("black", nrow(gwa.sim))
col[gwa.sim$chr %% 2 == 1] <- "gray50"
col[fdr < fdr.thresh] <- "green"
res <- manhattan(gwa.sim, pch = 20, col = col)
abline(v = res$coord[qtn.id], col = "pink")
}
is.close2QTN <- rep(F, nrow(bm.sim@snps)) # QTNおよびその周辺であれば真になる
is.detected <- rep(F, n.qtn) # そのQTNが検出されるかどうか
range <- 1e6 # QTNから1Mb以内は、正しく検出したものとする
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.close2QTN <- is.close2QTN | (cond.chr & cond.pos1 & cond.pos2) # QTNおよびその周辺
is.detected[i] <- sum(fdr < fdr.thresh & (cond.chr & cond.pos1 & cond.pos2)) > 0
}
n.positive <- sum(fdr < fdr.thresh) # 陽性であったSNP数
n.negative <- sum(!(is.close2QTN)) # 陰性となるべき（QTLの近くにない）SNP数
n.false.positive <- sum(!(is.close2QTN) & fdr < fdr.thresh)
fdr.obs <- n.false.positive / n.positive # 偽陽性のSNP数 / 陽性であったSNP数
fpr.obs <- n.false.positive / n.negative # 偽陽性のSNP数 / 陰性となるべきSNP数
n.undetected <- sum(!(is.detected)) # 検出されなかったQTN
fnr.obs <- sum(!(is.detected)) / n.qtn # 検出されなかったQTN数 / 全QTN数
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のある系統を抽出
bm.sim <- select.snps(bm.sim, maf > 0.05) # MAF > 5%のSNPのみを使う
bm.sim## A bed.matrix with 500 individuals and 36502 markers.
## snps stats are set
## ped stats are setstandardize(bm.sim) <- "mu_sigma"
grm <- GRM(bm.sim)
eigen.grm <- eigen(grm)以下が、シミュレーションの本体となる。settingsの数値を変えて実行することで、様々な結果が得られる。
settings <- list(
w.qtn = rep(1.0, 6), # 主働遺伝子の効果と数。ここでは、6QTNの効果を全て1とする。
g.var.ratio = 0.5, #ポリジーン遺伝分散対QTN遺伝分散比。ここでは50%とする。
he = 0.5, # 遺伝率。ここでは、表現型分散の50%が遺伝分散とする。
fdr.thresh = 0.05, # 偽発見率のしきい値。ここでは、5%とする。
outfile = "out/sim_1.0x5_1.0_0.7_0.1" # manhattan plotのファイル名
)
n.sim <- 20 # シミュレーションの回数
# 結果を入れるための入れ物を準備
sim.res <- data.frame(fnr = rep(F, n.sim),
fdr = rep(F, n.sim),
fpr = rep(F, n.sim))
# シミュレーションを実行
if(!dir.exists(settings$outfile)) { # 出力フォルダがなければ
dir.create(settings$outfile) # 出力フォルダの作成
}
for(i in 1:n.sim) { # 1〜n.simまでiを変化させて、ループを回す
print(i) # i回目であることを表示
# 上で定義した関数を使って計算。結果をsim.resに代入
jpeg(paste0(settings$outfile, "/", i, ".jpeg"), width = 640) # 図の書き出しの準備
sim.res[i, ] <- gwas.sim(bm.sim, grm, eigen.grm, settings) # シミュレーション
dev.off() # 図の書き出しの終了
}## [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 1.0000000 1.00000000 5.643022e-05
## 2 0.6666667 0.00000000 0.000000e+00
## 3 0.6666667 0.08333333 2.826935e-05
## 4 0.3333333 0.10714286 8.459760e-05
## 5 0.5000000 0.00000000 0.000000e+00
## 6 1.0000000 NaN 0.000000e+00
## 7 0.1666667 0.03448276 2.834387e-05
## 8 0.5000000 0.00000000 0.000000e+00
## 9 0.3333333 0.54320988 1.251992e-03
## 10 0.8333333 0.00000000 0.000000e+00
## 11 0.5000000 0.00000000 0.000000e+00
## 12 0.8333333 0.00000000 0.000000e+00
## 13 1.0000000 NaN 0.000000e+00
## 14 0.8333333 0.00000000 0.000000e+00
## 15 1.0000000 NaN 0.000000e+00
## 16 1.0000000 NaN 0.000000e+00
## 17 0.1666667 0.05555556 2.818807e-05
## 18 1.0000000 NaN 0.000000e+00
## 19 1.0000000 NaN 0.000000e+00
## 20 0.3333333 0.00000000 0.000000e+00# そののsummary
summary(sim.res) ## fnr fdr fpr
## Min. :0.1667 Min. :0.00000 Min. :0.000e+00
## 1st Qu.:0.4583 1st Qu.:0.00000 1st Qu.:0.000e+00
## Median :0.7500 Median :0.00000 Median :0.000e+00
## Mean :0.6833 Mean :0.13027 Mean :7.389e-05
## 3rd Qu.:1.0000 3rd Qu.:0.07639 3rd Qu.:2.821e-05
## Max. :1.0000 Max. :1.00000 Max. :1.252e-03
## NA's :6# 箱ひげ図を描く
title <- paste(sapply(settings, paste, collapse=":"), collapse=" ")
boxplot(sim.res, ylim = c(0,1), main = title) # 箱ひげ図