Licenses
CC BY 4.0
はじめに
対象の生物種における、個人/系統/品種/遺伝資源からゲノム多型情報と特定の表現型を取得すれば、その表現型の違い(自然変異；Natural variation) を説明する多型を回帰分析により推定することができる．この手法のことを、ゲノムワイド関連解析 (GWAS; Genome Wide Association Study) という．GWASにより、形質の原因遺伝子/QTLsやDNAマーカーの同定が期待できる．
GWASの計算モデルについて
GWASでは、ゲノム多型データと表現型データの回帰解析を行う． 回帰モデルとしては、主に以下の２つが適用される．
一般化線形モデル (GLM, generalized linear model)
混合線形モデル (MLM, generalized linear mixed model)
MLM(Q+K)の場合、対象集団のゲノム多型データにおける集団構造(K)と血縁関係(Q)の情報をモデルに組み込む． つまり、対象集団において、分集団や血縁関係が認められる場合は、MLM(Q+K)のモデルでGWASを実行するのがベターである． 集団構造(K)と血縁関係(Q)の計算方法についても下記で解説する．
Step01_パッケージのインストール・ロードと解析準備
GWASを実行するための関数がまとめられたgastonパッケージを使用する． gastonのインストールが未実施の場合、install.packages関数により、まずはインストールを行う．
library関数を用いて、gastonを呼び出す．
また、本解析で出力される結果ファイルを格納するためのディレクトリ(フォルダ)を作成しておく．
Step02_ゲノムデータの読み込み
この演習では、「世界のイネ」コアコレクションのデータを用いる． 「世界のイネ」コアコレクション (通称WRC) とは、世界の代表的なイネ系統を揃えたセットであり、イネ69品種から構成される． この系統集団については、Tanaka et al., 2020. Plant Cell Physiol.において、全ゲノムリシーケンスが実施され、ゲノムと表現型データが公開されている．
gaston::read.vcf関数により、ゲノム多型データであるvcfファイル (WRC.Chr7.vcf.gz) を読み込む．ここでは、変数bmにvcfファイルの情報を格納している． 演習上、データサイズを減らすため、染色体7番 (Chr7) のデータのみにしている．
#read vcf
bm<-gaston::read.vcf("vcf/WRC.Chr7.vcf.gz")
bm@snps$id<-paste(bm@snps$chr,"_",bm@snps$pos,sep="")読み込んだvcfファイル (bm)の情報について確認する
まず、変数bmに存在する系統名を確認してみる． 「世界のイネ」コアコレクションの系統は、WRC01-99の系統番号が付いている．
## [1] "WRC01" "WRC02" "WRC03" "WRC04" "WRC05" "WRC06" "WRC07" "WRC09"
## [9] "WRC10" "WRC100" "WRC11" "WRC12" "WRC13" "WRC14" "WRC15" "WRC16"
## [17] "WRC17" "WRC18" "WRC19" "WRC20" "WRC21" "WRC22" "WRC23" "WRC24"
## [25] "WRC25" "WRC26" "WRC27" "WRC28" "WRC29" "WRC30" "WRC31" "WRC32"
## [33] "WRC33" "WRC34" "WRC35" "WRC36" "WRC37" "WRC38" "WRC39" "WRC40"
## [41] "WRC41" "WRC42" "WRC43" "WRC44" "WRC45" "WRC46" "WRC47" "WRC48"
## [49] "WRC49" "WRC50" "WRC51" "WRC52" "WRC53" "WRC55" "WRC57" "WRC58"
## [57] "WRC59" "WRC60" "WRC61" "WRC62" "WRC63" "WRC64" "WRC65" "WRC66"
## [65] "WRC67" "WRC68" "WRC97" "WRC98" "WRC99"変数bmでの多型情報(bm@snps)について確認してみる． 染色体 (chr) のどの位置 (pos) にどういう変異が存在しているか記載されている．
| chr | id | dist | pos | A1 | A2 | quality | filter | N0 | N1 | N2 | NAs | N0.f | N1.f | N2.f | NAs.f | callrate | maf | hz |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7 | 7_10802 | 0 | 10802 | G | A | 33599.3 | . | 38 | 1 | 30 | 0 | NA | NA | NA | NA | 1 | 0.4420290 | 0.0144928 |
| 7 | 7_11514 | 0 | 11514 | T | C | 84766.4 | . | 15 | 0 | 54 | 0 | NA | NA | NA | NA | 1 | 0.2173913 | 0.0000000 |
| 7 | 7_11559 | 0 | 11559 | G | GA | 15259.4 | . | 53 | 0 | 16 | 0 | NA | NA | NA | NA | 1 | 0.2318841 | 0.0000000 |
| 7 | 7_11619 | 0 | 11619 | G | A | 34148.0 | . | 47 | 2 | 20 | 0 | NA | NA | NA | NA | 1 | 0.3043478 | 0.0289855 |
| 7 | 7_11876 | 0 | 11876 | G | T | 11893.8 | . | 63 | 0 | 6 | 0 | NA | NA | NA | NA | 1 | 0.0869565 | 0.0000000 |
| 7 | 7_12045 | 0 | 12045 | G | A | 14380.0 | . | 54 | 0 | 15 | 0 | NA | NA | NA | NA | 1 | 0.2173913 | 0.0000000 |
| 7 | 7_12234 | 0 | 12234 | A | G | 98912.9 | . | 4 | 0 | 65 | 0 | NA | NA | NA | NA | 1 | 0.0579710 | 0.0000000 |
| 7 | 7_12687 | 0 | 12687 | G | A | 20892.2 | . | 58 | 1 | 10 | 0 | NA | NA | NA | NA | 1 | 0.1521739 | 0.0144928 |
| 7 | 7_12795 | 0 | 12795 | G | A | 21320.0 | . | 58 | 1 | 10 | 0 | NA | NA | NA | NA | 1 | 0.1521739 | 0.0144928 |
| 7 | 7_12932 | 0 | 12932 | C | T | 19598.7 | . | 58 | 1 | 10 | 0 | NA | NA | NA | NA | 1 | 0.1521739 | 0.0144928 |
次に、これら多型情報をテーブルで数値化した情報を確認してみる． データを情報処理する際は、バイナリーな数字で表現したほうが扱いやすいので、各多型を以下のような、多型スコア行列で表現するのが一般的である．
「0」or「2」のバイナリーな変数で各系統の多型が表現されている． 「0」はREFホモ (A/A: 参照配列と同じ塩基配列がホモ) であり、「2」はALTホモ (B/B: 参照配列と異なる塩基配列がホモ) を意味する． 今回のデータではほとんど存在しないが、多型がヘテロの場合(A/B)は、「1」で表現する．
ちなみに、イネは自殖性の植物なので、ほとんどの変異がホモで固定されている．
| 7_10802 | 7_11514 | 7_11559 | 7_11619 | 7_11876 | 7_12045 | 7_12234 | 7_12687 | 7_12795 | 7_12932 | |
|---|---|---|---|---|---|---|---|---|---|---|
| WRC01 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| WRC02 | 0 | 2 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 |
| WRC03 | 0 | 2 | 2 | 0 | 0 | 0 | 2 | 0 | 0 | 0 |
| WRC04 | 0 | 2 | 2 | 0 | 0 | 2 | 2 | 0 | 0 | 0 |
| WRC05 | 2 | 2 | 0 | 2 | 0 | 0 | 2 | 0 | 0 | 0 |
Step03_ゲノム多型データの集団構造 (Q) の確認
GWASでは、解析集団の遺伝的集団構造が偽陽性(実際には原因遺伝子ではないのに、GWASで検出されてしまう)を生み出す原因となる． この問題を考慮するために、上述したMLMが適用される． そこで、解析集団の集団構造について確認しておく必要がある．
まず、多型データについて、集団中のマイナーアレル頻度(MAF; Minor Allele Frequency)が5%を上回るものだけを抽出する．MAFが小さい多型は、レアバリアントかノイズか統計的に判定することが難しく、MAF > 0.05を満たすような集団内でよくみられる多型のみを対象とするためである．
gaston::select.snps関数により、MAF > 0.05を満たす多型を抽出し、変数bm5に格納する．
prcomp関数により、ゲノム多型の行列データに対して、主成分分析(PCA; Principal component analysis)を行う．
第一主成分(PC1)と第二主成分(PC2)についてプロットしてみる．
ざっくりと分けると、プロットが3つに固まっているのがわかる(クラスター)． つまり、この解析集団(WRC)には、遺伝的な分集団が存在することが考えられる． Tanaka et al., 2020. Plant Cell Physiol.でも報告されているとおり、世界中のイネの遺伝背景を反映するWRCは3つの分集団を形成している． イネ (Oryza sativa L.) は、その種内でSub-populationとして、japonica型、indica型、aus型の3つの分集団に大別される． したがって、WRCでは、これら3つの分集団が反映されているため、GWASの計算モデルにおいて、分集団 (Q) を考慮する必要がある．
PCA結果における各主成分(PC)の説明率について確認する．
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 208.520 150.4870 80.84470 71.58684 67.75733 66.43833
## Proportion of Variance 0.302 0.1573 0.04539 0.03559 0.03189 0.03066
## Cumulative Proportion 0.302 0.4593 0.50466 0.54025 0.57214 0.60280
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 60.05267 56.74225 51.38184 49.26134 47.00669 45.59467
## Proportion of Variance 0.02505 0.02236 0.01834 0.01685 0.01535 0.01444
## Cumulative Proportion 0.62784 0.65020 0.66854 0.68539 0.70074 0.71518
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 44.55439 42.55259 42.15352 41.53076 41.26226 39.81280
## Proportion of Variance 0.01379 0.01258 0.01234 0.01198 0.01182 0.01101
## Cumulative Proportion 0.72896 0.74154 0.75388 0.76586 0.77769 0.78869
## PC19 PC20 PC21 PC22 PC23 PC24
## Standard deviation 38.02992 37.14676 35.82275 34.52836 33.81134 33.2920
## Proportion of Variance 0.01004 0.00958 0.00891 0.00828 0.00794 0.0077
## Cumulative Proportion 0.79874 0.80832 0.81724 0.82552 0.83346 0.8411
## PC25 PC26 PC27 PC28 PC29 PC30
## Standard deviation 32.71154 32.65372 31.54804 30.90553 29.42656 29.18277
## Proportion of Variance 0.00743 0.00741 0.00691 0.00663 0.00601 0.00591
## Cumulative Proportion 0.84858 0.85599 0.86290 0.86954 0.87555 0.88147
## PC31 PC32 PC33 PC34 PC35 PC36
## Standard deviation 28.84213 28.49068 27.96846 27.78149 27.27637 26.64315
## Proportion of Variance 0.00578 0.00564 0.00543 0.00536 0.00517 0.00493
## Cumulative Proportion 0.88724 0.89288 0.89831 0.90367 0.90884 0.91377
## PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 26.2998 25.53650 25.26069 24.71791 24.40961 23.91712
## Proportion of Variance 0.0048 0.00453 0.00443 0.00424 0.00414 0.00397
## Cumulative Proportion 0.9186 0.92310 0.92754 0.93178 0.93592 0.93989
## PC43 PC44 PC45 PC46 PC47 PC48
## Standard deviation 23.78521 23.13783 22.4468 22.14411 22.06729 21.72095
## Proportion of Variance 0.00393 0.00372 0.0035 0.00341 0.00338 0.00328
## Cumulative Proportion 0.94382 0.94754 0.9510 0.95444 0.95782 0.96110
## PC49 PC50 PC51 PC52 PC53 PC54
## Standard deviation 21.1254 20.83648 20.46686 19.94377 19.58438 19.20497
## Proportion of Variance 0.0031 0.00302 0.00291 0.00276 0.00266 0.00256
## Cumulative Proportion 0.9642 0.96722 0.97013 0.97289 0.97555 0.97811
## PC55 PC56 PC57 PC58 PC59 PC60
## Standard deviation 18.75045 18.1969 17.87525 17.7979 17.3781 16.89119
## Proportion of Variance 0.00244 0.0023 0.00222 0.0022 0.0021 0.00198
## Cumulative Proportion 0.98056 0.9829 0.98507 0.9873 0.9894 0.99135
## PC61 PC62 PC63 PC64 PC65 PC66
## Standard deviation 15.57301 15.36727 14.64631 14.36091 12.39385 9.51609
## Proportion of Variance 0.00168 0.00164 0.00149 0.00143 0.00107 0.00063
## Cumulative Proportion 0.99304 0.99468 0.99617 0.99760 0.99867 0.99930
## PC67 PC68 PC69
## Standard deviation 9.40607 3.59905 6.685e-12
## Proportion of Variance 0.00061 0.00009 0.000e+00
## Cumulative Proportion 0.99991 1.00000 1.000e+00累積寄与率についてプロットしてみる．
PC8-10くらいからカーブが寄与率がプラトーに近づきはじめているようにみえる．
続いて、この累積寄与率について変化量(Δ)を計算して、プロットしてみる．
deltacum<-c()
for (i in 2:length(cumeff)){
deltacum<-c(deltacum,cumeff[i]-cumeff[i-1])
}
plot(deltacum,type="b")
この寄与率の変化量が小さくなり始めているところで、ある程度主成分得点の寄与率は収束している．つまり、ここでは、PC8-10までの情報で、ゲノム多型情報の集団構造をおおかた説明していると考えられる． 以下にも示すが、今回の解析では、PC8(累積寄与率65%)までの変数をGWAS計算のMLMに組み込むことにする．
Step04_表現型データの読み込み
read.csv関数により、表現型データ (WRC_pheno.csv) を読み込む．
解析対象としたい表現型を抽出する． ここでは、米の種子色形質である「Grain_color」の列を抽出する． 「Grain_color」列は、変数phenoの4列目である．
すべての系統について、表現型データが得れていない場合もあるので、変数bm5から、表現型データが存在する系統の多型データを変数bm.wpに取り出す． また、解析に用いる最終的な表現型値のベクトルをyに格納しておく．
以下の通り、今回のデータでは、69系統における205,798の多型データ (SNPs/Indel) を解析に用いる．
## A bed.matrix with 69 individuals and 205798 markers.
## snps stats are set
## ped stats are setStep05_ゲノム関係 (血縁関係) 行列 (K) の計算
上述したとおり、MLM(Q+K)での解析では、集団構造(Q)と血縁関係(K)の情報が必要である． ゲノム関係行列を計算する前に、多型スコア行列であるbm.wpを標準化しておく必要がある．
gaston::standardize関数により標準化を行う．
ゲノム関係行列は、gaston::GRMにより計算する． この計算結果は、grmに格納しておく．
これで、GWASの計算に必要なデータが出揃った．次のステップでGWASを実行する．
Step06_GWASの実行
gaston::association.test関数により、GWASを実行する． オプションは以下の通りである．
- x : 多型スコア行列 (
bm.wp) - Y : 表現型値のベクトル (
y) - method : “lm”(GLM, generalized linear model: 線形モデル) or “lmm”(MLM, generalized linear mixed model: 線形混合モデル)
- response : “quantitative” or “binary”
- eigenK : ゲノム関係行列 (K)(
grm) - p : 主成分得点の数 (Q)
gwa_lmm <- gaston::association.test(x = bm.wp,
Y = y,
method = "lmm",
response = "quantitative",
test = "lrt",
eigenK=eigen(grm),
p=8)変数gwa_lmmにGWASの解析結果が格納されている．
GWASの結果をgaston::manhattan関数を用いて、マンハッタンプロットにより可視化する． ちなみに蛇足ですが、マンハッタンプロットという名称は、プロットの様子がマンハッタンのビル群に似ているから名付けられたそうです．
全多型マーカーをプロットするとデータサイズが重くなるので、ここでは、p<0.1を満たすものだけを抽出してプロットしている．
#marker filtering
plotlmm<-gwa_lmm[gwa_lmm$p<0.1,]
#manhattan plot
gaston::manhattan(plotlmm,
pch = 20,
chrom.col = c("deepskyblue4", "orangered3"),
main = "Grain_color_Chr7"
)
この図は、縦軸に-log10(p-value)、横軸が染色体上の位置を示しており、各プロットが多型である． 縦軸は、p-valueにlogをとって、マイナスを付けているので、値が高いほど統計的に確からしいことを意味する． 統計的に確からしい、つまりp-value (帰無仮説が成立する確率)が低い多型(SNPs/Indel)が表現型の差異に寄与している多型だと考えられる．
この解析では、染色体7番において、一つ顕著なピークが観察されている． したがって、このゲノム領域に原因遺伝子が座乗している可能性が考えられる．
作図したマンハッタンプロットをpdf関数で保存しておく．
#作図をpdfで保存
pdf("GWAS/GWAS_result_WRC_Chr7/Manhattanplot_WRC_GrainColor.pdf", width = 6, height = 5)
#manhattan plot
gaston::manhattan(plotlmm,
pch = 20,
chrom.col = c("deepskyblue4", "orangered3"),
main = "Grain_color_Chr7"
)
dev.off()続いて、実際のGWAS結果が格納されたデータテーブルをみてみる． すべての結果を閲覧するとデータサイズが大きいので、ここでは関連上位100多型のみを抽出して、GWAS結果を確認してみる．
#関連上位100多型に関する結果を抽出
gwa_lmm <- gwa_lmm[order(gwa_lmm$p),]
gwa_lmm_top100 <- gwa_lmm[1:100,]
#結果の確認
head(gwa_lmm_top100)| chr | pos | id | A1 | A2 | freqA2 | h2 | LRT | p | |
|---|---|---|---|---|---|---|---|---|---|
| 41026 | 7 | 6084536 | 7_6084536 | T | C | 0.3405797 | 0 | 79.15878 | 0 |
| 41028 | 7 | 6084550 | 7_6084550 | A | G | 0.3405797 | 0 | 79.15878 | 0 |
| 40992 | 7 | 6068071 | 7_6068071 | A | AACGCGAAAAGTCGG | 0.3478261 | 0 | 75.46361 | 0 |
| 40971 | 7 | 6064303 | 7_6064303 | A | G | 0.3478261 | 0 | 68.92314 | 0 |
| 40972 | 7 | 6064370 | 7_6064370 | G | A | 0.3478261 | 0 | 68.92314 | 0 |
| 40976 | 7 | 6065094 | 7_6065094 | A | G | 0.3478261 | 0 | 68.92314 | 0 |
結果の表でも示されている通り、染色体7番(chr7)の6,064,000-6,068,000bpあたりに存在する多型が原因だと考えられる． 特に、染色体7番の6,068,071bp (7_6068071)に、14bpの挿入配列(insertion)が検出されている． 3の倍数でない挿入配列なので、コドン読み枠がずれるフレームシフト変異により、タンパク機能に影響を及ぼしている可能性が考えられる．
この変異の効果や、どのような遺伝子における変異なのかを同定することが次のステップとなる．これについては、各種データベースや遺伝子情報などを参考に同定していく．
本解析で用いたイネコアコレクションについては、TASUKE+ (WRC+JRC ver.) というゲノム多型情報と遺伝子情報を閲覧できるブラウザが開発されており、これを利用することで、ユーザフレンドリーに解析することができる． これについては、別資料で解説する．
この結果がまとまったテーブルについて、write.csv関数で出力し、保存しておく．
Step07_QQ-plotの作成
GWASによる統計解析が妥当かどうか (偽陽性が考えられるか) 判定するために、QQ-plot(quantile-quantile plot)という手法が用いられる．
QQ-plotでは、x軸とy軸に、以下をプロットする．
- 帰無仮説から期待される-log10(p-value) (expected)
- GWASの多型で実際に計測された-log10(p-value) (observed)
これらに関連性が認められない場合、45°のプロットになるが、真の関連性がある場合、ほとん どの結果は45°の線上に収まるが、一部は外れて上方にプロットされる．
gaston::qqplot.pvalues関数でQQ-plotを作成する．
結果が示す通り、一部のプロットが上方に外れており、十分偽陽性が抑制されたモデルでGWASが実行されていることが示されている．
作図したQQ-plotをpdf関数で保存しておく．