GWASチュートリアル
本チュートリアルでは、Rを用いて 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データの読み込み、GWA解析にgastonを用いる。Rのパッケージインストーラで簡単にダウンロード・インストールできる。
require(gaston) # vcfフィアルの読み込みと、gwasの実行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イネ遺伝資源アクセッション1,568系統についての、38,769 SNPsの遺伝子型データが含まれている。
bmには、SNP遺伝子型が数値としてコードされたものが含まれている。
# 行列として最初の10サンプルの、最初の5 SNPの情報を表示する
as.matrix(bm[1:10, 1:5]) ## 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.
## IRGC121316@c88dcbba.0 0
## IRGC121578@950ba1f2.0 0
## IRGC121491@0b7d031b.0 0
## IRGC121440@d02dfe4e.0 0
## IRGC121896@89688646.0 0
## IRGC121323@57c71a3b.0 0
## IRGC121251@392016d6.0 0
## IRGC121902@f1be6ab1.0 0
## IRGC124483@9a9187e6.0 0
## IRGC121314@6c2404e4.0 0SNP遺伝子型が、“0”と“2”でコードされていることが分かる。また、4列目の上から８番目のように、“1”でコードされているSNP遺伝子型もある。“0”, “1”, “2”は、それぞれ、A1/A1、A1/A2、A2/A2の遺伝子型を表し、“0”と“2”はホモ接合体、“1”はヘテロ接合体である。なお、イネは自殖のため、ヘテロ接合の頻度は低い。
なお、A1やA2の対立遺伝子が、それぞれ、どの塩基に対応するのかは、SNP情報で確認できる。
# 最初の10 SNPの情報を示す
head(bm@snps, 10)## chr id dist pos A1 A2 quality filter N0 N1 N2 NAs N0.f N1.f
## 1 1 SNP-1.8562. 0 9563 A T 0 PASS 1343 7 218 0 NA NA
## 2 1 SNP-1.18594. 0 19595 G A 0 PASS 1481 8 79 0 NA NA
## 3 1 SNP-1.24921. 0 25922 C T 0 PASS 1348 24 196 0 NA NA
## 4 1 SNP-1.25253. 0 26254 A T 0 PASS 1153 9 406 0 NA NA
## 5 1 SNP-1.29213. 0 30214 T A 0 PASS 1346 4 218 0 NA NA
## 6 1 SNP-1.30477. 0 31478 C T 0 PASS 1013 15 540 0 NA NA
## 7 1 SNP-1.31732. 0 32733 T G 0 PASS 1133 10 425 0 NA NA
## 8 1 SNP-1.32666. 0 33667 A G 0 PASS 129 1 1438 0 NA NA
## 9 1 SNP-1.40457. 0 41458 C T 0 PASS 1341 9 218 0 NA NA
## 10 1 SNP-1.75851. 0 76852 G A 0 PASS 1346 4 218 0 NA NA
## N2.f NAs.f callrate maf hz
## 1 NA NA 1 0.14126276 0.0044642857
## 2 NA NA 1 0.05293367 0.0051020408
## 3 NA NA 1 0.13265306 0.0153061224
## 4 NA NA 1 0.26179847 0.0057397959
## 5 NA NA 1 0.14030612 0.0025510204
## 6 NA NA 1 0.34917092 0.0095663265
## 7 NA NA 1 0.27423469 0.0063775510
## 8 NA NA 1 0.08258929 0.0006377551
## 9 NA NA 1 0.14190051 0.0057397959
## 10 NA NA 1 0.14030612 0.0025510204bm@snpsには、染色体番号（chr）、マーカー名（id）、物理位置（pos）、A1とA2の対立遺伝子（A1,A2)、マイナーアリル頻度（maf）などの情報が含まれている。
系統/品種の属性情報と表現型データの読み込み
次に、おなじ遺伝資源アクセッションについて得られた表現型データを読み込む。
# 表現型データを読み込む
pheno <- read.csv("data/RiceDiversityPheno4GWASGS.csv", row.names = 1)
str(pheno) # 読み込んだデータの構造を示す## '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形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
ここでは、表現型データとして玄米の縦横比（Brown.rice.length.with.ration）を用いる。玄米の縦横比のデータを抜き出した後、各要素（各系統に対応）に名前を付ける。その後、na.omit関数を用いて欠測値を除く。
y <- pheno[, "Brown.rice.length.width.ratio"] # "###"を解析対象とする形質とする
# 後にSNP遺伝子型と対応付けるために、yの要素のそれぞれに名前をつけておく
names(y) <- rownames(pheno)
y <- na.omit(y) # 欠測しているところを除く
length(y) # yのサイズを示す## [1] 361解析に必要なマーカー遺伝子型データの抽出
実は、全ての系統について表現型データが得られているわけではない（例えば、玄米の縦横比では361系統のデータが得られている）。そこで、bmから、表現型データをもつ系統のSNPデータを抜き出す。
# 表現型データがある系統のデータだけを抜き出す
bm.wp <- bm[bm@ped$id %in% names(y),]
bm.wp # データのsummaryを示す## 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と同じ順番で並び替える。
# 表現型をSNP遺伝子型と同じ順番に並び替える
y <- y[bm.wp@ped$id] GWASの準備
GWASを行うための準備を行う。
まず、マイナーアリル頻度（MAF）が5%以上のSNPだけを選ぶ。MAFが小さいSNPは、ノイズのような表現型の摂動にも当てはまってしまうために除いておく。
# mafが0.05以上のものだけを抜き出す
bm.wp <- select.snps(bm.wp, maf > 0.05)
bm.wp # データの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となるように標準化する。
standardize(bm.wp) <- "mu_sigma"
grm <- GRM(bm.wp)GWASの実行
準備したデータを用いて、GWASを実行する。
gwa <- association.test(bm.wp, Y = y, method = "lmm", response = "quantitative", test = "wald", eigenK = eigen(grm), p = 4)ここで、簡単に関数の引数について説明しておく。
bm.wpは解析に用いるSNPデータ、Y = yは、解析に用いる表現型データを指定している。method = “lmm”は、線形混合モデルをあてはめることを意味する。通常、lmmを指定して解析する。response = “quantitative”は、対象としている形質が量的形質であることを意味する。また、test = は、効果の検定法を表し、尤度比検定（Likelihood Ratio Test: LRT)lrtの他に、“wald”, “score”を選ぶことができる。ここでは、“wald”とする。eigenK = eigen(grm)、および、p = 4は、いずれも集団構造による偽陽性を抑えるための効果についてのオプションであり、前者は、ゲノム関係行列（の固有値分解の結果）を、後者は、遺伝的背景を説明するための主成分を上位何主成分までモデルに含めるかを指定している。ここでは、上位4主成分を遺伝的背景を表す変数としてモデルに組み込んでいる。
では、先頭の10 SNPsについての結果をみてみる。
head(gwa, 10)## chr pos id A1 A2 freqA2 h2 beta sd
## 1 1 9563 SNP-1.8562. A T 0.13988920 0.9061210 0.008391307 0.04058884
## 2 1 25922 SNP-1.24921. C T 0.13434903 0.9062347 0.015191355 0.04055168
## 3 1 26254 SNP-1.25253. A T 0.29501385 0.9060718 0.003571187 0.04045482
## 4 1 30214 SNP-1.29213. T A 0.13711911 0.9061029 0.005529634 0.04111036
## 5 1 31478 SNP-1.30477. C T 0.22437673 0.9057650 0.005767405 0.02649280
## 6 1 32733 SNP-1.31732. T G 0.30332410 0.9059233 -0.003916660 0.04063249
## 7 1 33667 SNP-1.32666. A G 0.86980609 0.9036403 -0.033153583 0.02735686
## 8 1 41458 SNP-1.40457. C T 0.14127424 0.9060878 0.012587092 0.04027855
## 9 1 76852 SNP-1.75851. G A 0.13850416 0.9061264 0.009118406 0.04066683
## 10 1 80943 SNP-1.79942. C T 0.08033241 0.9058468 0.033856963 0.05116239
## p
## 1 0.8362135
## 2 0.7079452
## 3 0.9296574
## 4 0.8930016
## 5 0.8276652
## 6 0.9232090
## 7 0.2255535
## 8 0.7546597
## 9 0.8225844
## 10 0.5081283一番重要なのは、最後の列の\(p\)であり、これは、各SNPの効果の有意性検定の\(p\)値をあらわす。この値が小さいほど、効果が有意であることを意味する。
\(p\)値が小さい順に並べて、先頭10 SNPsについて結果をみてみる。
gwa[order(gwa$p)[1:10], ]## chr pos id A1 A2 freqA2 h2 beta
## 10562 3 16733441 SNP-3.16732086. G T 0.3725762 0.8788667 0.2533487
## 15910 5 5371772 SNP-5.5371749. G A 0.4778393 0.8606853 0.2057749
## 15911 5 5372955 SNP-5.5372932. A G 0.4903047 0.8689087 0.1784448
## 15904 5 5327913 SNP-5.5327890. A G 0.3850416 0.8804746 0.1865606
## 10561 3 16719368 SNP-3.16718013. G A 0.2562327 0.8924705 -0.2011179
## 15896 5 5249951 SNP-5.5249928. G A 0.4030471 0.8865718 0.1646579
## 15925 5 5516194 SNP-5.5516131. G A 0.6426593 0.8995778 0.1513896
## 15897 5 5265914 SNP-5.5265891. C T 0.4044321 0.8871808 0.1628942
## 10557 3 16688126 SNP-3.16686770. G A 0.3130194 0.8962656 -0.1510675
## 15920 5 5458379 SNP-5.5458356. A T 0.4418283 0.8930123 0.1674133
## sd p
## 10562 0.02665234 1.987375e-21
## 15910 0.02500894 1.902808e-16
## 15911 0.02443748 2.833153e-13
## 15904 0.02960887 2.960458e-10
## 10561 0.03560530 1.618182e-08
## 15896 0.02974932 3.115051e-08
## 15925 0.02745448 3.503629e-08
## 15897 0.02957011 3.614025e-08
## 10557 0.02782328 5.650034e-08
## 15920 0.03103676 6.889287e-08第3, 5染色体に特に高度に有意なSNPがあることが分かる。
GWASの結果を図示する。
次に、マンハッタンプロット（Manhattan plot）と呼ばれる図を描く。縦軸は、\(-log_{10}(p)\)を表す。\(p\)値が小さいと、\(-log_{10}(p)\)は大きくなる。
manhattan(gwa, pch = 20)
なお、非常に多くのSNPについて検定を行っているので、多重検定により生じる偽陽性をコントロールする必要がある。
多重検定に対して\(p\)値を補正して、偽陽性をコントロールする方法の一つとして、BenjaminiとHochber(1995)が提案した偽発見率（false discovery rate: FDR）がある。これは、検出された陽性の結果のうち、偽の陽性（偽陽性、本当は陰性だが陽性として検出されているもの）が含まれる割合の期待値を計算して、それに基づき偽陽性をコントロールする方法である。ここでは、偽発見率（FDR）が5%以下になるように、\(p\)値の閾値を調整する。
fdr <- p.adjust(gwa$p, method = "BH")
gwa[fdr < 0.05, ]## chr pos id A1 A2 freqA2 h2 beta
## 5016 2 2679995 SNP-2.2679992. C T 0.05263158 0.9443238 0.3694488
## 8883 3 1596846 SNP-3.1595841. A G 0.35318560 0.9260435 0.2121902
## 10557 3 16688126 SNP-3.16686770. G A 0.31301939 0.8962656 -0.1510675
## 10561 3 16719368 SNP-3.16718013. G A 0.25623269 0.8924705 -0.2011179
## 10562 3 16733441 SNP-3.16732086. G T 0.37257618 0.8788667 0.2533487
## 10572 3 16780671 SNP-3.16779544. A T 0.30332410 0.8955280 -0.1741554
## 15895 5 5247991 SNP-5.5247968. C A 0.40027701 0.8907179 0.1444242
## 15896 5 5249951 SNP-5.5249928. G A 0.40304709 0.8865718 0.1646579
## 15897 5 5265914 SNP-5.5265891. C T 0.40443213 0.8871808 0.1628942
## 15899 5 5298292 SNP-5.5298269. G A 0.11634349 0.8991583 0.2535089
## 15901 5 5303070 SNP-5.5303047. C T 0.26315789 0.8941448 0.1408904
## 15902 5 5307410 SNP-5.5307387. G C 0.26315789 0.8937388 0.1374524
## 15903 5 5321857 SNP-5.5321834. G A 0.12188366 0.8997981 0.2423396
## 15904 5 5327913 SNP-5.5327890. A G 0.38504155 0.8804746 0.1865606
## 15905 5 5338205 SNP-5.5338182. T C 0.27562327 0.8907620 0.1801285
## 15907 5 5362484 SNP-5.5362461. G A 0.21052632 0.8950362 -0.1836103
## 15910 5 5371772 SNP-5.5371749. G A 0.47783934 0.8606853 0.2057749
## 15911 5 5372955 SNP-5.5372932. A G 0.49030471 0.8689087 0.1784448
## 15915 5 5388363 SNP-5.5388340. C T 0.11080332 0.8939576 -0.1587393
## 15917 5 5431060 SNP-5.5431037. G A 0.30609418 0.8989630 0.1474012
## 15919 5 5442823 SNP-5.5442800. G A 0.24515235 0.8919045 0.1839881
## 15920 5 5458379 SNP-5.5458356. A T 0.44182825 0.8930123 0.1674133
## 15925 5 5516194 SNP-5.5516131. G A 0.64265928 0.8995778 0.1513896
## 15928 5 5547500 SNP-5.5547437. G T 0.49445983 0.9017788 0.1436518
## 17886 5 28533147 SNP-5.28470501. A G 0.20914127 0.8973239 0.2001066
## 17887 5 28548001 SNP-5.28485355. C T 0.07063712 0.8967777 0.2420380
## 22972 7 22057564 SNP-7.22056570. A G 0.56786704 0.8942887 -0.1859165
## 22977 7 22113010 SNP-7.22112016. A T 0.10941828 0.8905355 0.2233007
## 23002 7 22252173 SNP-7.22251179. A T 0.08725762 0.8952676 0.2046841
## 23081 7 23048478 SNP-7.23047484. A G 0.11495845 0.8922213 0.1688938
## 23254 7 24845026 SNP-7.24844031. C T 0.09556787 0.8950154 0.2021734
## 23271 7 25214651 SNP-7.25213656. C T 0.08587258 0.8890976 0.2796072
## 26130 8 24258192 SNP-8.24255477. G A 0.12326870 0.9012496 0.1657622
## 26187 8 24948137 SNP-8.24945422. G T 0.16066482 0.9328161 0.2777845
## sd p
## 5016 0.07368416 5.332328e-07
## 8883 0.04699912 6.338538e-06
## 10557 0.02782328 5.650034e-08
## 10561 0.03560530 1.618182e-08
## 10562 0.02665234 1.987375e-21
## 10572 0.04071803 1.893372e-05
## 15895 0.03089362 2.941161e-06
## 15896 0.02974932 3.115051e-08
## 15897 0.02957011 3.614025e-08
## 15899 0.05918868 1.843370e-05
## 15901 0.03385445 3.159540e-05
## 15902 0.03364821 4.407975e-05
## 15903 0.05211250 3.314260e-06
## 15904 0.02960887 2.960458e-10
## 15905 0.03499722 2.647748e-07
## 15907 0.03404130 6.900232e-08
## 15910 0.02500894 1.902808e-16
## 15911 0.02443748 2.833153e-13
## 15915 0.03791274 2.827086e-05
## 15917 0.03560261 3.470347e-05
## 15919 0.03869302 1.983764e-06
## 15920 0.03103676 6.889287e-08
## 15925 0.02745448 3.503629e-08
## 15928 0.03460767 3.312190e-05
## 17886 0.04121932 1.205809e-06
## 17887 0.04768319 3.855429e-07
## 22972 0.04133729 6.874128e-06
## 22977 0.04364940 3.124481e-07
## 23002 0.04806738 2.060116e-05
## 23081 0.04117478 4.098155e-05
## 23254 0.04603103 1.122551e-05
## 23271 0.05288558 1.243234e-07
## 26130 0.03902775 2.163708e-05
## 26187 0.06454360 1.678730e-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を実施する。method = “lm”、test = “wald”、p=0とする。
gwa0 <- association.test(bm.wp, Y = y, method = "lm", response = "quantitative", test = "wald", p = 0)マンハッタンプロットを描いてみる。
fdr <- p.adjust(gwa0$p, method = "BH")
#gwa0[fdr < 0.05, ]
col <- rep("black", nrow(gwa0))
col[gwa0$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa0, pch = 20, col = col)
\(-log_{10}(p)\)の値が全体に非常に大きくなり、関連遺伝子の検出を目的とするGWA解析が全く機能しなくなる。GWASでは、遺伝的背景を適切にモデルに組み込み、偽陽性を抑え込むことが極めて重要である。