バイオメトリックス第6回 ゲノムワイドアソシエーション研究(GWAS）
例として利用するデータ
ここでは、Rを用いてGWASおよびGSモデル構築を行う方法を解説する。 なお、Zhaoら（2011; Nature Communications 2:467）がイネ遺伝資源のGWASに用いた表現型データと、 McCouchら（2016; Nature Communication 7:10532）が用いたHigh Density Rice Array (HDRA)を用いてジェノタイピングされたSNP遺伝子型データ（いずれも、Rice Diversityからダウンロードできる）を、例として用いる。 後者のデータは、70万SNPsのデータのうち、他のマーカーとの冗長性が無く、遺伝子型データの欠測率が低く（< 1%）、かつ、マイナーアリル頻度（minor allele frequency: MAF）が高い（> 5%）、38,769 SNPsのデータを用いる。なお、欠測値については、Beagle 5.1（Blowning et al. 2018）を用いて補完をしてある。
利用するRのパッケージ
SNPデータの読み込み・解析にgastonとpcaMethodsを、GWASにはgastonパッケージを用いる。
require(gaston) # for reading a vcf file and performing GWASGWASの実行
データの読み込み
まず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込む。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用いる。
# read a vcf file
bm <- read.vcf("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 # show the summary of the file## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are setイネ遺伝資源アクセッション1568系統についての、38,769 SNPsの遺伝子型データが含まれている。
bmには、SNP遺伝子型が数値としてコードされたものが含まれている。
# show the first 10 rows and 10 columns of marker score data
as.matrix(bm[1:10, 1:10])## SNP-1.8562. SNP-1.18594. SNP-1.24921. SNP-1.25253.
## IRGC121316@c88dcbba.0 0 0 0 0
## IRGC121578@950ba1f2.0 0 0 0 0
## IRGC121491@0b7d031b.0 0 0 0 2
## IRGC121440@d02dfe4e.0 0 0 0 2
## IRGC121896@89688646.0 0 0 0 0
## IRGC121323@57c71a3b.0 0 0 0 0
## IRGC121251@392016d6.0 0 0 0 0
## IRGC121902@f1be6ab1.0 0 0 0 1
## IRGC124483@9a9187e6.0 0 0 0 0
## IRGC121314@6c2404e4.0 0 0 0 0
## SNP-1.29213. SNP-1.30477. SNP-1.31732. SNP-1.32666.
## IRGC121316@c88dcbba.0 0 2 0 2
## IRGC121578@950ba1f2.0 0 0 0 2
## IRGC121491@0b7d031b.0 0 0 2 2
## IRGC121440@d02dfe4e.0 0 0 2 2
## IRGC121896@89688646.0 0 0 0 2
## IRGC121323@57c71a3b.0 0 0 0 2
## IRGC121251@392016d6.0 0 0 0 2
## IRGC121902@f1be6ab1.0 0 0 1 2
## IRGC124483@9a9187e6.0 0 2 0 2
## IRGC121314@6c2404e4.0 0 0 0 2
## SNP-1.40457. SNP-1.75851.
## IRGC121316@c88dcbba.0 0 0
## IRGC121578@950ba1f2.0 0 0
## IRGC121491@0b7d031b.0 0 0
## IRGC121440@d02dfe4e.0 0 0
## IRGC121896@89688646.0 0 0
## IRGC121323@57c71a3b.0 0 0
## IRGC121251@392016d6.0 0 0
## IRGC121902@f1be6ab1.0 0 0
## IRGC124483@9a9187e6.0 0 0
## IRGC121314@6c2404e4.0 0 0SNP遺伝子型が、“0”と”2”でコードされていることが分かる。また、4列目の上から８番目のように、“1”でコードされているSNP遺伝子型もある。“0”, “1”, “2”は、それぞれ、A1/A1、A1/A2、A2/A2の遺伝子型を表し、“0”と”2”はホモ接合体、“1”はヘテロ接合体である。なお、イネは自殖であるため、ヘテロ接合の頻度は低い。
なお、A1やA2の対立遺伝子が、それぞれ、どの塩基に対応するのかは、SNP情報で確認できる。
# show the first 10 SNPs
head(bm@snps, 10)bm@snpsには、染色体番号（chr）、マーカー名（id）、物理位置（pos）、A1とA2の対立遺伝子（A1, A2)、3つのSNP遺伝子型（A1/A1, A1/A2, A2/A2）の観察数（N0, N1, N2）、非欠測率（callrate）、マイナーアリル頻度（maf）などの情報が含まれている。
系統/品種の属性情報と表現型データの読み込み
次に、おなじ遺伝資源アクセッションについて得られた表現型データを読み込む。
# read phenotypic data
pheno <- read.csv("RiceDiversityPheno.csv", row.names = 1)
str(pheno) # show the data summary## 'data.frame': 388 obs. of 36 variables:
## $ Flowering.time.at.Arkansas : num 75.1 89.5 94.5 87.5 89.1 ...
## $ Flowering.time.at.Faridpur : int 64 66 67 70 73 NA 74 75 55 NA ...
## $ Flowering.time.at.Aberdeen : int 81 83 93 108 101 158 122 81 74 NA ...
## $ FT.ratio.of.Arkansas.Aberdeen : num 0.927 1.078 1.016 0.81 0.882 ...
## $ FT.ratio.of.Faridpur.Aberdeen : num 0.79 0.795 0.72 0.648 0.723 ...
## $ Culm.habit : num 4 7.5 6 3.5 6 3 2.5 6.5 3 NA ...
## $ Leaf.pubescence : int 1 0 1 1 1 1 NA 1 NA NA ...
## $ Flag.leaf.length : num 28.4 39 27.7 30.4 36.9 ...
## $ Flag.leaf.width : num 1.283 1 1.517 0.892 1.75 ...
## $ Awn.presence : int 0 0 0 0 0 1 0 1 1 NA ...
## $ Panicle.number.per.plant : num 3.07 4.05 3.12 3.7 2.86 ...
## $ Plant.height : num 111 144 128 154 148 ...
## $ Panicle.length : num 20.5 26.8 23.5 29 30.9 ...
## $ Primary.panicle.branch.number : num 9.27 9.17 8.67 6.36 11.17 ...
## $ Seed.number.per.panicle : num 4.79 4.44 5.08 4.52 5.54 ...
## $ Florets.per.panicle : num 4.91 4.73 5.15 4.73 5.68 ...
## $ Panicle.fertility : num 0.879 0.75 0.935 0.813 0.871 0.907 0.781 0.887 0.537 NA ...
## $ Seed.length : num 8.06 7.71 8.24 9.71 7.12 ...
## $ Seed.width : num 3.69 2.95 2.93 2.38 3.28 ...
## $ Seed.volume : num 2.59 2.11 2.15 1.94 2.2 ...
## $ Seed.surface.area : num 3.91 3.66 3.71 3.7 3.66 ...
## $ Brown.rice.seed.length : num 5.79 5.6 6.12 7.09 5.14 ...
## $ Brown.rice.seed.width : num 3.11 2.57 2.47 2.04 2.8 ...
## $ Brown.rice.surface.area : num 3.51 3.29 3.33 3.3 3.28 ...
## $ Brown.rice.volume : num 7.74 5.11 5.13 4.09 5.56 ...
## $ Seed.length.width.ratio : num 2.19 2.61 2.81 4.08 2.17 ...
## $ Brown.rice.length.width.ratio : num 1.86 2.18 2.48 3.47 1.83 ...
## $ Seed.color : int 0 0 0 0 0 0 0 0 0 NA ...
## $ Pericarp.color : int 0 0 0 0 0 0 0 0 0 NA ...
## $ Straighthead.suseptability : num 4.83 7.83 NA 8.33 8.17 ...
## $ Blast.resistance : int 8 4 3 5 3 2 2 9 2 1 ...
## $ Amylose.content : num 15.6 23.3 23.1 19.3 23.2 ...
## $ Alkali.spreading.value : num 6.08 5.64 5.53 6.03 5.44 ...
## $ Protein.content : num 8.45 9.2 8 9.6 8.5 ...
## $ Year07Flowering.time.at.Arkansas: num 73 88 105.5 86.5 85.5 ...
## $ Year06Flowering.time.at.Arkansas: num 77.2 91 83.5 88.5 92.7 ...phenoには、36形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
ここでは、形質例として種子の長さを用いる。種子の長さのデータを抜き出した後、各要素（各系統に対応）に名前を付ける。その後、na.omit関数を用いて欠測値を除く。
y <- pheno[, "Seed.length"] # select seed length as the target trait
names(y) <- rownames(pheno) # naming the elements of y
y <- na.omit(y) # remove missing elements解析に必要なマーカー遺伝子型データの抽出
実は、全ての系統について表現型データが得られているわけではない。そこで、bmから、表現型データをもつ系統のSNPデータを抜き出す。
# select samples with ids contained in rownames(pheno)
bm.wp <- bm[bm@ped$id %in% names(y),]
bm.wp # show the summary of the selected data## A bed.matrix with 361 individuals and 38769 markers.
## snps stats are set
## There is one monomorphic SNP
## ped stats are set361品種/系統のマーカー遺伝子型が抜き出される。
最後に、表現型データ\(y\)を、 SNPデータbm.wpと同じ順番で並び替える。
y <- y[bm.wp@ped$id] # reorder the elements of y as same as the bm@ped$id分集団構造解析
次に、356系統の遺伝的背景を調べるために、マーカー遺伝子型データに基づく主成分分析を行う。マーカー遺伝子型データは、上述したように”0”, “1”, “2”でスコア化され、行列のようなかたちでbm内に保存されている。
このデータを取り出し、その主成分分析を行う。ここでは、pcaMethodsパッケージのpca関数を利用する予定であったが、RStudio Cloudでは計算に長い時間を要するため、prcompを用いることとした。なお、pca関数では、欠測値がある場合でも、主成分分析を行うことができる。
res <- prcomp(as.matrix(bm.wp)) # use prcomp for pca
pcs <- res$x[,1:4]
rm(res) # remove the PCA result because the object is large in size.なお、主成分分析の結果は、サイズが大きいので、主成分得点を抽出した後に削除しておく。
主成分分析の結果を図示する。
op <- par(mfrow = c(1,2))
plot(pcs[,1:2]) # a scatter plot graph between pc1 and pc2
plot(pcs[,3:4]) # a scatter plot graph between pc3 and pc4
par(op)遺伝的背景と形質の表現型の間に関連があると、GWASで偽陽性を生じる（実際、GWASでは、遺伝的背景の違いによる偽陽性が大きな問題なる）。ここでは、それを視覚化して確認する。
op <- par(mfrow = c(1,2))
symbols(pcs[,1:2], circle = y - min(y), inches = 0.1)
symbols(pcs[,3:4], circle = y - min(y), inches = 0.1)
par(op)遺伝的背景との関係は少なからずありそうである。
GWASの準備
ここからは、GWASを行う。まずは、そのための準備を行う。
まず、マイナーアリル頻度（MAF）が5%以上のSNPsだけを選ぶ。
bm.wp <- select.snps(bm.wp, maf > 0.05) # filter with MAF
bm.wp # show the data summary## A bed.matrix with 361 individuals and 36112 markers.
## snps stats are set
## ped stats are setこれにより、36,112のSNPマーカーがGWAS解析用に残る。
次に、遺伝的背景を説明するために、ゲノム関係行列を計算しておく。この行列を含めて（具体的には、この行列に規定される共分散構造をもつ多変量正規分布に従う効果を含める、詳細は省く）解析することで、偽陽性を抑え込むことができる。ゲノム関係行列を計算する前に、SNPスコア行列であるbm.wpを標準化（standardize）しておくしておく必要がある。ここでは、平均0、標準偏差1となるように標準化する（なお、自分で標準化を行わない場合は、平均\(2p\)、標準偏差はHardy-Weinberg平衡に従うときの期待値\(2p(1-p)\)で標準化される。ここで、\(p\)はA2対立遺伝子の頻度である）。
standardize(bm.wp) <- "mu_sigma"
#grm <- GRM(bm.wp)
#head(grm, 3)[,1:3]
grm <- tcrossprod(as.matrix(bm.wp)) / ncol(bm.wp)
head(grm, 3)[,1:3]## NSFTV80@02d095ba.0 NSFTV87@64fa0112.0 NSFTV96@32a6808e.0
## NSFTV80@02d095ba.0 1.08064106 -0.03870567 -0.04754667
## NSFTV87@64fa0112.0 -0.03870567 0.63496932 0.34106057
## NSFTV96@32a6808e.0 -0.04754667 0.34106057 0.63313979GWASの実行
準備しておいたデータを用いて、GWASを実行する。
gwa <- association.test(bm.wp, Y = y,
method = "lmm",
response = "quantitative",
test = "lrt",
eigenK = eigen(grm),
p = 4)gasonのassociation.test関数は計算が速く、短い時間で計算が終わる。
ここで、簡単に関数の引数について説明しておく。
bm.wpは解析に用いるSNPデータ、Y = yは、解析に用いる表現型データを指定している。method = “lmm”は、線形混合モデルをあてはめることを意味する。通常、lmmを指定して解析する。response = “quantitative”は、対象としている形質が量的形質であることを意味する。質的な二値形質の場合は、response = “binary”とする。また、test = “lrt”は、効果の検定法を表し、尤度比検定（Likelihood Ratio Test: LRT)を指定している。この他に、“wald”, “score”を選ぶことができる。これら2つの手法は、LRTよりも計算速度は速いが、計算に時間が問題にならないときには“lrt”を指定するとよい。eigenK = eigen(grm)、および、p = 4は、いずれも集団構造による偽陽性を抑えるための効果についてのオプションであり、前者は、ゲノム関係行列（の固有値分解の結果）を、後者は、遺伝的背景を説明するための主成分を上位何主成分までモデルに含めるかを指定している。ここでは、上位4主成分を遺伝的背景を表す変数としてモデルに組み込んでいる。
では、先頭の10 SNPsについての結果をみてみる。
head(gwa, 10)一番重要なのは、最後の列の\(p\)であり、これは、各SNPの効果の有意性検定の\(p\)値をあらわす。この値が小さいほど、効果が有意であることを意味する。
\(p\)値が小さい順に並べて、先頭10 SNPsについて結果をみてみる。
gwa[order(gwa$p)[1:10], ]第3染色体に高度に有意なSNPがあることが分かる。
GWASの結果を図示する。
次に、マンハッタンプロット（Manhattan plot）と呼ばれる図を描く。縦軸は、\(-log_{10}(p)\)を表す。\(p\)値が小さいと、\(-log_{10}(p)\)は大きくなる。
manhattan(gwa, pch = 20)
なお、非常に多くのSNPについて検定を行っているので、多重検定により生じる偽陽性をコントロールする必要がある。
多重検定に対して\(p\)値を補正して、偽陽性をコントロールする方法の一つとして、BenjaminiとHochberg (1995)が提案した偽発見率（false discovery rate: FDR）がある。これは、検出された陽性の結果のうち、偽陽性が含まれる割合の期待値を計算して、それに基づき偽陽性をコントロールする方法である。ここでは、偽陽性率（FDR）が5%以下になるように、\(p\)値の閾値を調整する。
fdr <- p.adjust(gwa$p, method = "BH")
gwa[fdr < 0.05, ]主に第3染色体と第5染色体に、有意なSNPが有ることが分かる。
この結果を、図示して確認する。有意なSNPを緑色で表す。
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)
分集団構造（遺伝的背景）による影響
GWASに用いる集団に分集団構造がある場合には、分集団構造が存在することにより偽陽性が生じやすくなる。
ここでは、主成分得点をGWASのモデルに組み込まない場合、あるいは、ゲノム関係行列をGWASのモデルに組み込まない場合について、陽性となるSNPがどのように増加するかを確認する。なお、ゲノム上の多くのSNPは、対象とする形質の表現型に関わっていない場合がほとんどである。そのような場合、形質の表現型に関わらないSNPの\(p\)値は一様分布に従うと期待される。そこで、\(p\)値が一様分布に従うという仮定のもとでQ-Q（Quantile-Quantile）プロットを描き、一様分布から逸脱するSNPを確認してみる。
まずは、上述したモデルでQ-Qプロットを描く。
qqplot.pvalues(gwa$p, pch = 20)
一部、逸脱するSNPがあるが、これらは表現型を支配する遺伝子の中、あるいは、近傍にあり、一様分布に従わないSNPと考えられる。
一方、主成分を一つも組み込まないモデルでGWASを行う。これは、p = 0とすればよい。
gwa.wopc <- association.test(bm.wp, Y = y,
method = "lmm",
response = "quantitative",
test = "lrt",
eigenK = eigen(grm),
p = 0)Q-Qプロットで、\(p\)値を図示する。
qqplot.pvalues(gwa.wopc$p, pch = 20)
結果は、主成分を含むモデルとほとんど変わらない。つまり、この形質では、主成分の効果をモデルに入れても入れなくても結果は大きく変わらないと考えられる。
実際にFDRを計算して、マンハッタンプロットを描いてみる。
fdr <- p.adjust(gwa.wopc$p, method = "BH")
col <- rep("black", nrow(gwa.wopc))
col[gwa.wopc$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa.wopc, pch = 20, col = col)
結果はほとんど変わらない。
次に、ゲノム関係行列で説明される効果を含まないモデルでGWASを実施する。method = “lm”、test = “wald”とする。
gwa.wogrm <- association.test(bm.wp, Y = y,
method = "lm",
response = "quantitative",
test = "wald")Q-Qプロットで、\(p\)値を図示する。
qqplot.pvalues(gwa.wogrm$p, pch = 20)
ほとんど全てのSNPで、\(p\)が一様分布から逸脱していることがわかる。このように、遺伝的背景を無視して、単純にSNPの遺伝子型と、形質の表現型の関連を解析すると、偽陽性が非常に多くなることが分かる。
マンハッタンプロットを描いてみる。
fdr <- p.adjust(gwa.wogrm$p, method = "BH")
col <- rep("black", nrow(gwa.wogrm))
col[gwa.wogrm$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa.wogrm, pch = 20, col = col)
染色体全体に有意なSNPが検出されてしまう。