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

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

GWASの実行

データの読み込み

　まず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込む。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用いる。

# 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

GWASのシミュレーション

まず、疑似乱数の種をセットしておく。これにより、以下の解析で同じ結果が得られる。

set.seed(20200619)

データの準備。計算負荷を小さくするためにデータ数を減らす。ここでは、客観的な方法（k-medoids法）を用いて予め選んでおいた代表的な500系統を解析に用いる。

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 > 0.05

繰り返し利用することになる血縁行列とその固有値分解を予め行っておく。

standardize(bm.sim) <- "mu_sigma"
grm <- GRM(bm.sim)
grm <- grm + diag(0.0001, nrow(grm))
eigen.grm <- eigen(grm)

全SNPから、QTN (quantitative trait nucleotides)を無作為に選ぶ。

n <- ncol(bm.sim) # 全SNP数
n.qtn <- 5 # QTNの数
qtn.id <- sample(n, n.qtn) # n.qtn個のSNPをQTNとして無作為に選ぶ
X.qtn <- as.matrix(bm.sim[, qtn.id]) # QTNの遺伝子型スコアを抜き出す

次に、シミュ絵r−ションにおけるQTNの効果を決め、それらQTNによる遺伝変異を全系統について計算する。

w.qtn <- rep(1.0, n.qtn) # 主働遺伝子の効果を、ここでは全て1.0とする
g.qtn <- as.vector(X.qtn %*% w.qtn) # QTNによる変異を全系統について計算
var(g.qtn) # QTN効果による遺伝分散
## [1] 4.454911

効果の小さな遺伝子（ポリジーン、微働遺伝子）による効果を擬似的に生成する。生成には、多変量正規分布からの乱数を用いる。QTNによる遺伝変異の大きさに応じて発生させるように、先にQTNによる遺伝変異の分散（遺伝分散）を計算し、その分散に応じてポリジーンによる効果の大きさを調節する。

# QTNの遺伝分散の何割かをポリジーンによる遺伝分散とする
g.var.ratio <- 0.5 # QTNによる遺伝分散の50%をポリジーンによる遺伝分散にする
# 多変量正規乱数としてポリジーン効果を生成
g.poly <- as.vector(rmvnorm(1, sigma = grm * var(g.qtn) * g.var.ratio)) 
var(g.poly) # ポリジーン効果による分散
## [1] 1.693201

QTNの効果と微働遺伝子の効果を足し合わせて全系統の遺伝変異（遺伝子型値）を計算する。

g <- g.qtn + g.poly # 遺伝子型値 = QTNの効果 + ポリジーン効果
var(g) # 遺伝分散
## [1] 4.972652

次に、指定した遺伝率に従って環境変異を発生させる。

he <- 0.7 # 遺伝率。ここでは、表現型分散の70%が遺伝分散とする。
var.e <- (1 - he) / he * var(g) # 環境分散を計算する
e <- rnorm(length(g), sd = sqrt(var.e)) # var.eを分散にもつ環境変異を生成
var(e) # 環境分散
## [1] 2.104551

最後に、遺伝変異（遺伝子型値）に環境変異を加え、表現型変異（表現型値）を計算する。

y <- g + e # 遺伝変異に環境変異を加えて表現型変異を計算する
var(y) # 表現型分散
## [1] 6.971512

擬似的に生成した表現型値\(ymat\)の1列を用いてGWASを行ってみる。

# GWASの実行
gwa.sim <- association.test(bm.sim, Y = y, 
                            method = "lmm", response = "quantitative", 
                            test = "wald", eigenK = eigen.grm, p = 4)

マンハッタンプロットを描く。

fdr.thresh <- 0.05 # 偽発見率のしきい値
fdr <- p.adjust(gwa.sim$p, method = "BH")
col <- rep("black", nrow(gwa.sim))
col[gwa.sim$chr %% 2 == 1] <- "gray50"
col[fdr < fdr.thresh] <- "green"
manhattan(gwa.sim, pch = 20, col = col)

bm.sim@snps[qtn.id,]
##       chr              id dist      pos A1 A2 quality filter  N0 N1  N2 NAs
## 36652  12 SNP-12.4748568.    0  4749584  A  T       0   PASS 480  3  17   0
## 6770    2 SNP-2.12770342.    0 12770347  T  A       0   PASS 372  6 122   0
## 24025   7 SNP-7.16883224.    0 16884218  G  T       0   PASS 427  6  67   0
## 21490   6 SNP-6.22395465.    0 22396463  T  A       0   PASS 311  9 180   0
## 22431   6 SNP-6.29904766.    0 29905765  C  A       0   PASS 433  4  63   0
##       N0.f N1.f N2.f NAs.f callrate   maf    hz
## 36652   NA   NA   NA    NA        1 0.037 0.006
## 6770    NA   NA   NA    NA        1 0.250 0.012
## 24025   NA   NA   NA    NA        1 0.140 0.012
## 21490   NA   NA   NA    NA        1 0.369 0.018
## 22431   NA   NA   NA    NA        1 0.130 0.008

検出できたQTNの数を調べて、偽陰性率を計算する。

# 検出できたQTN数を調べる
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の割合
fnr.obs # 偽陰性率
## [1] 0.2

偽発見率、偽陽性率を計算する。なお、QTNの近くにあるSNPで検出されたものは真の陽性と判断することとする。

is.true <- rep(F, nrow(bm.sim@snps)) # 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.true <- is.true | (cond.chr & cond.pos1 & cond.pos2) # QTNおよびその周辺
}
fdr.obs <- sum(!(is.true) & fdr < fdr.thresh) / sum(fdr < fdr.thresh)
fdr.obs # 偽発見率
## [1] 0.0952381
fpr.obs <- sum(!(is.true) & fdr < fdr.thresh) / sum(!(is.true))
fpr.obs # 偽陽性率
## [1] 0.0001058705