GWAS Simulation
SNPデータの読み込み、GWA解析にgastonを用いる。Rのパッケージインストーラで簡単にダウンロード・インストールできる。また、多変量正規分布にしたがう乱数を発生させるためにmvtnormを用いる。こちらも読み込んでおく。
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 setGWASのシミュレーション
まず、疑似乱数の種をセットしておく。これにより、以下の解析で同じ結果が得られる。
set.seed(20230525)データの準備。計算負荷を小さくするためにデータ数を減らす。ここでは、客観的な方法（k-medoids法）を用いて予め選んでおいた代表的な500系統を解析に用いる。なお、k-medioids法での代表選びのためのRコードはhttps://rpubs.com/hiroiwata/627490にある。
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%以下のものを除いておく
bm.sim## A bed.matrix with 500 individuals and 36502 markers.
## snps stats are set
## ped stats are set500系統の36,502マーカーのデータとなっている。
マーカースコアを平均0、分散1に基準化しておく。
standardize(bm.sim) <- "mu_sigma" # 基準化繰り返し利用することになる血縁行列とその固有値分解を予め行っておく。
grm <- tcrossprod(as.matrix(bm.sim)) / nrow(bm.sim@snps)
eigen.grm <- eigen(grm)全SNPから、QTN (quantitative trait nucleotides)を無作為に選ぶ。
n <- ncol(bm.sim) # 全SNP数
n.qtn <- 6 # QTNの数
qtn.id <- sample(n, n.qtn) # n.qtn個のSNPをQTNとして無作為に選ぶ
X.qtn <- as.matrix(bm.sim[, qtn.id]) # QTNの遺伝子型スコアを抜き出す次に、シミュレーションにおける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] 13.18241効果の小さな遺伝子（ポリジーン、微働遺伝子）による効果を擬似的に生成する。生成には、多変量正規分布からの乱数を用いる。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] 7.190041QTNの効果と微働遺伝子の効果を足し合わせて全系統の遺伝変異（遺伝子型値）を計算する。
g <- g.qtn + g.poly # 遺伝子型値 = QTNの効果 + ポリジーン効果
var(g) # 遺伝分散## [1] 11.1873次に、指定した遺伝率に従って環境変異を発生させる。
he <- 0.5 # 遺伝率。ここでは、表現型分散の50%が遺伝分散とする。
var.e <- (1 - he) / he * var(g) # 環境分散を計算する
e <- rnorm(length(g), sd = sqrt(var.e)) # var.eを分散にもつ環境変異を生成
var(e) # 環境分散## [1] 11.08311最後に、遺伝変異（遺伝子型値）に環境変異を加え、表現型変異（表現型値）を計算する。
y <- g + e # 遺伝変異に環境変異を加えて表現型変異を計算する
var(y) # 表現型分散## [1] 21.38792擬似的に生成した表現型値\(y\)を用いて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"
res <- manhattan(gwa.sim, pch = 20, col = col)
QTNとして検出されたSNPを示す。
res[fdr < fdr.thresh,]## chr pos id A1 A2 freqA2 h2 beta sd
## 12789 4 1074839 SNP-4.1070389. T C 0.334 0.4724265 1.225496 0.2414696
## 21953 7 7091066 SNP-7.7090070. T A 0.141 0.4724131 1.550123 0.3121123
## 21960 7 7167883 SNP-7.7166887. A G 0.143 0.4589641 1.494236 0.3077564
## 21965 7 7216442 SNP-7.7215446. T C 0.135 0.4682119 1.598557 0.3118301
## 21969 7 7239623 SNP-7.7238627. G A 0.126 0.4662441 2.004560 0.3211789
## p coord
## 12789 3.871776e-07 0.3129213
## 21953 6.815229e-07 0.5878182
## 21960 1.202320e-06 0.5880242
## 21965 2.953815e-07 0.5881545
## 21969 4.340735e-10 0.5882167次に、QTNとしたSNPの検定の結果、および、マンハッタンプロットのx軸上での位置を確認する。
res[qtn.id,]## chr pos id A1 A2 freqA2 h2 beta sd
## 12789 4 1074839 SNP-4.1070389. T C 0.334 0.4724265 1.2254961 0.2414696
## 14209 4 20465566 SNP-4.20293609. T C 0.703 0.4655281 1.1702477 0.5042255
## 3472 1 31596132 SNP-1.31595087. C T 0.738 0.4835603 1.1185380 0.4758560
## 25902 8 19423857 SNP-8.19421143. C A 0.682 0.4960053 0.7317241 0.2897714
## 21969 7 7239623 SNP-7.7238627. G A 0.126 0.4662441 2.0045597 0.3211789
## 23346 7 23926645 SNP-7.23925650. G A 0.590 0.4532506 1.4282800 0.3416838
## p coord
## 12789 3.871776e-07 0.31292131
## 14209 2.029323e-02 0.36493828
## 3472 1.874414e-02 0.08475882
## 25902 1.156401e-02 0.70049712
## 21969 4.340735e-10 0.58821669
## 23346 2.913528e-05 0.63298079QTNとしたSNPの位置をピンク色の垂線で表す。
res <- manhattan(gwa.sim, pch = 20, col = col)
abline(v = res$coord[qtn.id], col = "pink") # 垂線を引く
偽発見率、偽陽性率を計算する。なお、QTNの近くにあるSNPで検出されたものは真の陽性と判断することとする。
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.positive## [1] 5n.negative <- sum(!(is.close2QTN)) # 陰性となるべき（QTLの近くにない）SNP数
n.negative## [1] 35182n.false.positive <- sum(!(is.close2QTN) & fdr < fdr.thresh)
# 偽陽性のSNP数（本当は陰性なのに陽性になったSNP数）
n.false.positive## [1] 0fdr.obs <- n.false.positive / n.positive # 偽陽性のSNP数 / 陽性であったSNP数
fdr.obs # 偽発見率## [1] 0fpr.obs <- n.false.positive / n.negative # 偽陽性のSNP数 / 陰性となるべきSNP数
fpr.obs # 偽陽性率## [1] 0n.undetected <- sum(!(is.detected))
# 検出されなかったQTN数（近くにSNPが検出されなかったQTN）
n.undetected## [1] 4fnr.obs <- sum(!(is.detected)) / n.qtn # 検出されなかったQTN数 / 全QTN数
fnr.obs # 偽陰性率## [1] 0.6666667