連鎖不平衡、GWAS
講義の内容
本講義では、Rを用いてGWASを実際に行い、分集団構造による偽陽性や、偽陽性に対応すためのモデル、および、多数のSNPに関する検定法について説明を行います。
本講義では、Iwataら（2015; PLoS ONE 10:e0120610）)でがイネの種子形のゲノミック予測に用いた楕円フーリエ記述子のデータと、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データの読み込み、GWASにはgastonパッケージを用います。いずれも、Rのパッケージインストーラで簡単にダウンロード・インストールが可能です。階層構造をもつフォルダにあるデータファイルの指定に便利なhereもロードしておきます。
require(gaston) # for reading a vcf file and performing GWAS
require(here)
require(tidyverse)GWASの実行
データの読み込み
ゲノムワイドマーカーの遺伝子型データは、多くの場合、Variant Call Format (vcf)とよばれる形式で保存されています。今回もSNP遺伝子型データはvcf形式で保存されています。そこでま、ず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込みます。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用います。
# read a vcf file
bm <- read.vcf(here("data2", "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: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列目の上から8番目のように、“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)## 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)、3つのSNP遺伝子型（A1/A1, A1/A2, A2/A2）の観察数（N0, N1, N2）、非欠測率（callrate）、マイナーアリル頻度（maf）などの情報が含まれています。
系統/品種の属性情報と表現型データの読み込み
次に、楕円フーリエ記述子のデータと、遺伝資源のパスポートデータを読み込んでみましょう。
# read phenotypic data
ef.org <- read.csv(here("data2", "FourierDescriptor_datasetC.csv"), row.names = 1)
line.org <- read.csv(here("data2", "RiceDiversityLine4GWASGS.csv"), row.names = 1, stringsAsFactors = T)
head(ef.org)[,1:8] # show the data summary## a1 b1 c1 d1 a2 b2 c2
## 1 1 -2.70e-17 4.80e-17 0.6111694 0.001039338 -0.005764824 -0.009906241
## 3 1 -4.45e-18 9.03e-18 0.5004835 0.002268561 -0.010033800 -0.012319853
## 4 1 -1.95e-17 3.62e-17 0.4762106 0.004475037 -0.010257737 -0.016015754
## 5 1 -2.68e-17 8.04e-17 0.3628859 0.000684247 -0.005984455 -0.019866794
## 6 1 7.11e-18 -5.04e-18 0.6107227 0.003886584 0.002052005 0.001065063
## 7 1 -1.76e-17 3.64e-17 0.4894189 -0.001395340 -0.001280532 -0.013163296
## d2
## 1 0.010476213
## 3 0.006070962
## 4 0.007118144
## 5 -0.006931079
## 6 0.016307643
## 7 0.003520876head(line.org)## NSFTV.ID GSOR.ID IRGC.ID Accession.Name
## NSFTV1@86f75d2b.0 1 301001 To be assigned Agostano
## NSFTV3@5ef1be74.0 3 301003 117636 Ai-Chiao-Hong
## NSFTV4@81d03b86.0 4 301004 117601 NSF-TV 4
## NSFTV5@5533f406.0 5 301005 117641 NSF-TV 5
## NSFTV6@0d125c0e.0 6 301006 117603 ARC 7229
## NSFTV7@e37be9e5.0 7 301007 To be assigned Arias
## Country.of.origin Latitude Longitude Sub.population
## NSFTV1@86f75d2b.0 Italy 41.871940 12.56738 TEJ
## NSFTV3@5ef1be74.0 China 27.902527 116.87256 IND
## NSFTV4@81d03b86.0 India 22.903081 87.12158 AUS
## NSFTV5@5533f406.0 India 30.472664 75.34424 AROMATIC
## NSFTV6@0d125c0e.0 India 22.903081 87.12158 AUS
## NSFTV7@e37be9e5.0 Indonesia -0.789275 113.92133 TRJefには、77次元の楕円フーリエ記述子のデータが含まれています。 lineには、これら遺伝資源の起源国や、その緯度・経度、および、分集団への属性データが含まれています。
まず、efとlineを融合してみましょう。efの列名とlineのNSFTV.IDが対応しています。
line <- line.org %>% filter(NSFTV.ID %in% rownames(ef.org)) %>% arrange(NSFTV.ID)
head(line) ## NSFTV.ID GSOR.ID IRGC.ID Accession.Name
## NSFTV1@86f75d2b.0 1 301001 To be assigned Agostano
## NSFTV3@5ef1be74.0 3 301003 117636 Ai-Chiao-Hong
## NSFTV4@81d03b86.0 4 301004 117601 NSF-TV 4
## NSFTV5@5533f406.0 5 301005 117641 NSF-TV 5
## NSFTV6@0d125c0e.0 6 301006 117603 ARC 7229
## NSFTV7@e37be9e5.0 7 301007 To be assigned Arias
## Country.of.origin Latitude Longitude Sub.population
## NSFTV1@86f75d2b.0 Italy 41.871940 12.56738 TEJ
## NSFTV3@5ef1be74.0 China 27.902527 116.87256 IND
## NSFTV4@81d03b86.0 India 22.903081 87.12158 AUS
## NSFTV5@5533f406.0 India 30.472664 75.34424 AROMATIC
## NSFTV6@0d125c0e.0 India 22.903081 87.12158 AUS
## NSFTV7@e37be9e5.0 Indonesia -0.789275 113.92133 TRJef <- ef.org[as.character(line$NSFTV.ID),]
rownames(ef) <- rownames(line) # 行名をlineと同じにそろえてしまいます。
head(ef)[,1:8]## a1 b1 c1 d1 a2 b2
## NSFTV1@86f75d2b.0 1 -2.70e-17 4.80e-17 0.6111694 0.001039338 -0.005764824
## NSFTV3@5ef1be74.0 1 -4.45e-18 9.03e-18 0.5004835 0.002268561 -0.010033800
## NSFTV4@81d03b86.0 1 -1.95e-17 3.62e-17 0.4762106 0.004475037 -0.010257737
## NSFTV5@5533f406.0 1 -2.68e-17 8.04e-17 0.3628859 0.000684247 -0.005984455
## NSFTV6@0d125c0e.0 1 7.11e-18 -5.04e-18 0.6107227 0.003886584 0.002052005
## NSFTV7@e37be9e5.0 1 -1.76e-17 3.64e-17 0.4894189 -0.001395340 -0.001280532
## c2 d2
## NSFTV1@86f75d2b.0 -0.009906241 0.010476213
## NSFTV3@5ef1be74.0 -0.012319853 0.006070962
## NSFTV4@81d03b86.0 -0.016015754 0.007118144
## NSFTV5@5533f406.0 -0.019866794 -0.006931079
## NSFTV6@0d125c0e.0 0.001065063 0.016307643
## NSFTV7@e37be9e5.0 -0.013163296 0.003520876次に、楕円フーリエ記述子の主成分分析を行います。最初の3列（\(a_1, b_1, c_1\)）は定数なので解析から除きます。
ef.pca <- prcomp(ef[, -(1:3)])ここからは、形質例として楕円フーリエ記述子の第1主成分を用います。第1主成分得点を主成分分析の結果からベクトルとして抜き出し、その後、ベクトルの各要素（各系統の形質の表現型値に対応）に名前を付けます。その後、na.omit関数を用いて欠測値をもつ要素を除きます。
y <- ef.pca$x[,1] # select seed length with ratio as the target trait
names(y) <- rownames(ef.pca$x) # naming the elements of y解析に必要なマーカー遺伝子型データの抽出
マーカー遺伝子型をもつ全系統について表現型データが得られているわけではありません。そこで、全データであるbmから、表現型データをもつ系統のSデータを抜き出します。
bm.wp <- bm[bm@ped$id %in% names(y),] # select samples with ids contained in rownames(pheno)
bm.wp # show the summary of the selected data## A bed.matrix with 358 individuals and 38769 markers.
## snps stats are set
## There is one monomorphic SNP
## ped stats are setrm(bm) # 必要がなくなったbmをメモリ節約のために削除しておきます。この結果、358品種/系統のマーカー遺伝子型が抜き出されます。
最後に、表現型データ\(y\)とパスポートデータlineを、 SNPデータbm.wpと同じ順番で並び替えます。
y <- y[bm.wp@ped$id] # reorder the elements of y as same as the bm@ped$id
line.wp <- line[bm.wp@ped$id, ]このように、データの対応関係を常に意識して解析の準備をすることが重要です。
分集団構造解析
次に、361系統の遺伝的背景を調べるために、マーカー遺伝子型データに基づく主成分分析を行いましょう。マーカー遺伝子型データは、上述したように“0”, “1”, “2”でスコア化され、行列のようなかたちでbm内に保存されています。
このデータをもとに主成分分析を行うには、prcompのような主成分分析のための関数を用いてもいいですが、以前の講義で説明したように計算に時間がかかります。そこで、以前説明したように内積行列の固有値分解を通して、主成分分析を行います。
bm.wpから行列を取り出し、内積行列を計算してみましょう。なお、解析をおこなう前に多型の無いSNPをデータから除いておくひつようがあります。
X <- as.matrix(bm.wp)
X <- X[, apply(X, 2, sd) > 0] # 標準編纂が0以上、つまり多型がある列だけを取り出す。
XXt <- tcrossprod(scale(X))次に、内積行列の固有値分解を行い、主成分得点を計算します。
eigen.XXt <- eigen(XXt)
pcs <- eigen.XXt$vectors[,1:4] %*% diag(sqrt(eigen.XXt$values[1:4]))
colnames(pcs) <- paste0("PC", 1:ncol(pcs))
rm(X)
rm(XXt)
rm(eigen.XXt)なお、マーカースコア行列や内積行列、固有値分析の結果は、サイズが大きいので、主成分得点を計算した後に削除しておきます。
主成分分析の結果を、所属する分集団で色付けして、図示してみます。
op <- par(mfrow = c(1,2))
plot(pcs[,1:2], col = line.wp$Sub.population) # a scatter plot graph between pc1 and pc2
plot(pcs[,3:4], col = line.wp$Sub.population) # 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, fg = line.wp$Sub.population)
symbols(pcs[,3:4], circle = y - min(y), inches = 0.1, fg = line.wp$Sub.population)
par(op)遺伝的背景との関係がありそうです。
連鎖不平衡の確認
では、次に、SNP間の連鎖不平衡（linkage disequilibrium: LD）を計算してみましょう。 GWASにおける原因遺伝子検出には、原因遺伝子とSNP間にLDが生じていることで、SNPと形質間に間接的なアソシエーションが生じるという現象が利用されていたことを思い出しましょう。
連鎖不平衡は、gastonのLDという関数を用いて計算することができます。なお、LDにはいくつかの指標がありますが、授業で説明した\(r^2\)を計算してみます。全SNP間でLDを計算すると膨大な数になりますので、ここでは、第1染色体上のSNPでLDを計算します。
x <- select.snps(bm.wp, chr == 1)
ld.x <- LD(x, c(1, ncol(x)), measure = "r2")では、物理的距離とLD（\(r2\)）間の関係を見てみましょう。全ての組合せについて視覚化することは難しいので、例として、100番目のSNPとそれ以外のSNP間での関係を図示してみます。
snp.id <- 100 # 100, 200, 300 ...
d.x <- abs(x@snps$pos[snp.id] - x@snps$pos) # 物理距離
plot(d.x, ld.x[snp.id,])
近傍のSNPと高いLDを示していますが、近傍のSNPでもLDが低い場合や、遠くにあるSNPとのLDが高い場合があることが分かります。snp.idを100以外に変えて見てみると、様々なパターンがあることが分かります。
op <- par(mfrow = c(3, 3))
for(i in 1:9) {
snp.id <- 100 * i # 100, 200, 300 ...
d.x <- abs(x@snps$pos[snp.id] - x@snps$pos)
plot(d.x, ld.x[snp.id,], main = snp.id)
}
par(op)遺伝資源などで見られるLDのパターンは、このように連鎖にだけ依存していません（連鎖不平衡、LD）が現象につけられた名前であり、原因につけられた名前でないことに注意しましょう）。
rm(list = c("x", "ld.x", "d.x")) # メモリ節約のためGWASの準備
ここからは、GWASを行います。まずは、そのための準備をします。
まず、マイナーアリル頻度（minor allele frequency: MAF, 2対立遺伝子（アリル）のうち頻度が小さいほうのアリルの頻度）が5%以上のSNPsだけを選びます。これは、MAFの小さいSNPは、GWASにおいてその多型が表現型変異に当てはまりやすく、偽陽性を生じやすいためです。
bm.wp <- select.snps(bm.wp, maf > 0.05) # filter with MAF
bm.wp # show the data summary## A bed.matrix with 358 individuals and 36463 markers.
## snps stats are set
## ped stats are setこれにより、36,463のSNPマーカーがGWAS解析用に残ります。
次に、遺伝的背景を説明するために、ゲノム関係行列を計算します。この行列を含めて解析することで、偽陽性を抑え込むことができます。ゲノム関係行列を計算する前に、SNPスコア行列であるbm.wpを標準化（standardize）しておくしておく必要があります。ここでは、平均0、標準偏差1となるように標準化します（なお、自分で標準化を行わない場合は、平均\(2p\)、標準偏差はHardy-Weinberg平衡に従うときの期待値\(\sqrt{2p(1-p)}\)で標準化されます。ここで、\(p\)はA2対立遺伝子の頻度です。）。
standardize(bm.wp) <- "mu_sigma"
grm <- GRM(bm.wp)GWASの実行
準備しておいたデータを用いて、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”を選ぶことができます。eigenK = eigen(grm)、および、p = 4は、いずれも集団構造による偽陽性を抑えるための効果についてのオプションであり、前者は、ゲノム関係行列（の固有値分解の結果）を、後者は、遺伝的背景を説明するための主成分を上位何主成分までモデルに含めるかを指定しています。ここでは、上位4主成分を遺伝的背景を表す変数としてモデルに組み込んでいます。
では、先頭の10 SNPsについての結果をみてみましょう。
head(gwa, 10)## chr pos id A1 A2 freqA2 h2 LRT p
## 1 1 9563 SNP-1.8562. A T 0.14106145 0.8717404 0.06893811 0.7928892
## 2 1 25922 SNP-1.24921. C T 0.13547486 0.8718531 0.06280951 0.8021093
## 3 1 26254 SNP-1.25253. A T 0.28631285 0.8721909 0.01531777 0.9015014
## 4 1 30214 SNP-1.29213. T A 0.13826816 0.8716741 0.08842823 0.7661849
## 5 1 31478 SNP-1.30477. C T 0.24022346 0.8722125 0.16376373 0.6857147
## 6 1 32733 SNP-1.31732. T G 0.29748603 0.8718513 0.05770915 0.8101541
## 7 1 33667 SNP-1.32666. A G 0.87709497 0.8701317 0.13743568 0.7108437
## 8 1 41458 SNP-1.40457. C T 0.14106145 0.8718214 0.04657251 0.8291385
## 9 1 76852 SNP-1.75851. G A 0.13966480 0.8717461 0.07670073 0.7818193
## 10 1 80943 SNP-1.79942. C T 0.09217877 0.8702694 0.88628493 0.3464860最も重要なのは、最後の列の\(p\)であり、これは、各SNPの効果の有意性検定の\(p\)値をあらわします。この値が小さいほど、効果が有意であることを意味します。
そこで、\(p\)値が小さい順に並べて、先頭10 SNPsについて改めて結果をみてみましょう。
gwa[order(gwa$p)[1:10], ]## chr pos id A1 A2 freqA2 h2 LRT
## 16070 5 5371772 SNP-5.5371749. G A 0.4734637 0.7827426 64.31886
## 10665 3 16733441 SNP-3.16732086. G T 0.3645251 0.8345868 55.20341
## 16071 5 5372955 SNP-5.5372932. A G 0.4832402 0.7927668 52.55020
## 16065 5 5327913 SNP-5.5327890. A G 0.3715084 0.8214964 36.54417
## 16085 5 5516194 SNP-5.5516131. G A 0.6396648 0.8549593 34.84871
## 16068 5 5362484 SNP-5.5362461. G A 0.2122905 0.8386957 29.64287
## 16080 5 5458379 SNP-5.5458356. A T 0.4245810 0.8497937 29.54718
## 16066 5 5338205 SNP-5.5338182. T C 0.2667598 0.8465115 25.90144
## 23197 7 22113010 SNP-7.22112016. A T 0.1103352 0.8410718 25.51440
## 16058 5 5265914 SNP-5.5265891. C T 0.3882682 0.8431956 23.36040
## p
## 16070 1.058278e-15
## 10665 1.086798e-13
## 16071 4.193902e-13
## 16065 1.492459e-09
## 16085 3.563460e-09
## 16068 5.194340e-08
## 16080 5.457201e-08
## 16066 3.593006e-07
## 23197 4.390937e-07
## 16058 1.343152e-06第3、5、7染色体に高度に有意な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, ]## chr pos id A1 A2 freqA2 h2 LRT
## 10660 3 16688126 SNP-3.16686770. G A 0.31145251 0.8575286 19.28535
## 10664 3 16719368 SNP-3.16718013. G A 0.25837989 0.8521074 21.43197
## 10665 3 16733441 SNP-3.16732086. G T 0.36452514 0.8345868 55.20341
## 16056 5 5247991 SNP-5.5247968. C A 0.37849162 0.8454511 19.01029
## 16057 5 5249951 SNP-5.5249928. G A 0.38826816 0.8421626 23.31025
## 16058 5 5265914 SNP-5.5265891. C T 0.38826816 0.8431956 23.36040
## 16060 5 5298292 SNP-5.5298269. G A 0.10893855 0.8553497 18.59121
## 16064 5 5321857 SNP-5.5321834. G A 0.11452514 0.8572776 19.99156
## 16065 5 5327913 SNP-5.5327890. A G 0.37150838 0.8214964 36.54417
## 16066 5 5338205 SNP-5.5338182. T C 0.26675978 0.8465115 25.90144
## 16068 5 5362484 SNP-5.5362461. G A 0.21229050 0.8386957 29.64287
## 16070 5 5371772 SNP-5.5371749. G A 0.47346369 0.7827426 64.31886
## 16071 5 5372955 SNP-5.5372932. A G 0.48324022 0.7927668 52.55020
## 16075 5 5388363 SNP-5.5388340. C T 0.10893855 0.8402135 21.83398
## 16079 5 5442823 SNP-5.5442800. G A 0.25000000 0.8377329 20.30884
## 16080 5 5458379 SNP-5.5458356. A T 0.42458101 0.8497937 29.54718
## 16085 5 5516194 SNP-5.5516131. G A 0.63966480 0.8549593 34.84871
## 16088 5 5547500 SNP-5.5547437. G T 0.46787709 0.8645725 19.49374
## 19181 6 9030490 SNP-6.9029490. A G 0.31983240 0.8613269 17.81457
## 23197 7 22113010 SNP-7.22112016. A T 0.11033520 0.8410718 25.51440
## 23486 7 24845026 SNP-7.24844031. C T 0.09776536 0.8565890 20.63318
## 23503 7 25214651 SNP-7.25213656. C T 0.08379888 0.8495891 19.00131
## p
## 10660 1.125667e-05
## 10664 3.666064e-06
## 10665 1.086798e-13
## 16056 1.300153e-05
## 16057 1.378637e-06
## 16058 1.343152e-06
## 16060 1.619656e-05
## 16064 7.778458e-06
## 16065 1.492459e-09
## 16066 3.593006e-07
## 16068 5.194340e-08
## 16070 1.058278e-15
## 16071 4.193902e-13
## 16075 2.972877e-06
## 16079 6.589528e-06
## 16080 5.457201e-08
## 16085 3.563460e-09
## 16088 1.009298e-05
## 19181 2.435112e-05
## 23197 4.390937e-07
## 23486 5.562373e-06
## 23503 1.306288e-05第3, 5, 6, 7染色体上に、有意な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\)値はその分布に従ったかたちで生じるため）。そこで、\(p\)値が一様分布に従うという仮定のもとでQ-Q（Quantile-Quantile）プロットを描き、一様分布から逸脱するSNPがどの程度あるのか確認してみる。
まずは、上述したモデルでQ-Qプロットを描いてみます。
qqplot.pvalues(gwa$p, pch = 20)
一部のSNPで、一様分布から\(p\)値が逸脱するようにみえますが、これらは表現型を支配する遺伝子の中、あるいは、近傍にあり、一様分布に従わないSNPと考えられます。
一方、主成分を一つも組み込まないモデルでGWASを行ってみましょう。これは、p = 0とすれば実施できます。
gwa.wopc <- association.test(bm.wp, Y = y, method = "lmm", response = "quantitative", test = "lrt", eigenK = eigen(grm), p = 0)まずはマンハッタンプロットを描いてみましょう。
fdr <- p.adjust(gwa.wopc$p, method = "BH")
col <- rep("black", nrow(gwa))
col[gwa.wopc$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa.wopc, pch = 20, col = col)
さきほどと結果は大きくは変わりません。
次に、Q-Qプロットで、\(p\)値を図示します。
qqplot.pvalues(gwa.wopc$p, pch = 20)
結果は、主成分を含むモデルとあまり大きくは変わっていません。つまり、この形質では、主成分の効果をモデルに入れても入れなくても結果が変わらないと考えられます。
次に、ゲノム関係行列で説明される効果を含まないモデルでGWASを実施してみます。ここでは、method = “lm”、test = “wald”とします。
gwa.wogrm <- association.test(bm.wp, Y = y, method = "lm", response = "quantitative", test = "wald")まずはマンハッタンプロットを描いてみましょう。
fdr <- p.adjust(gwa.wogrm$p, method = "BH")
col <- rep("black", nrow(gwa))
col[gwa.wogrm$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa.wogrm, pch = 20, col = col)
多くのマーカーでアソシエーションが有意です。これらは、ほとんどが偽陽性と考えられます。
では、Q-Qプロットで、\(p\)値を図示してみましょう。
qqplot.pvalues(gwa.wogrm$p, pch = 20)
すると、ほとんど全てのSNPで、\(p\)が一様分布から逸脱していることがわかります。このように、遺伝的背景を無視して、単純にSNPの遺伝子型と、形質の表現型の関連を解析すると、偽陽性が非常に多くなり、肝心の遺伝子との関連が見えなくなることが分かります。
効果の推定と図示
改めてGWASを行い、ゲノム関係行列で説明される効果を含まないモデルでGWASを実施してみます。ここでは、method = “lmm”、test = “wald”とします。
gwa <- association.test(bm.wp, Y = y, method = "lmm", response = "quantitative", test = "wald", eigenK = eigen(grm), p = 4)改めて、マンハッタンプロットを描いてみましょう。
fdr <- p.adjust(gwa$p, method = "BH")
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)
有意なSNPのうち、第3, 5, 7染色体上で最も\(-\log_{10}(p)\)が大きな3つのSNPについて、その効果を形状として再現して確認してみましょう。
id <- rep(NA, 3); beta <- rep(NA, 3)
tmp <- gwa %>% filter(chr == 3) %>% arrange(p) %>% slice_head(n = 1)
id[1] <- as.character(tmp$id)
beta[1] <- tmp$beta
tmp <- gwa %>% filter(chr == 5) %>% arrange(p) %>% slice_head(n = 1)
id[2] <- as.character(tmp$id)
beta[2] <- tmp$beta
tmp <- gwa %>% filter(chr == 7) %>% arrange(p) %>% slice_head(n = 1)
id[3] <- as.character(tmp$id)
beta[3] <- tmp$beta
g <- as.matrix(bm.wp)[, id] %*% diag(beta)
head(g)## [,1] [,2] [,3]
## NSFTV80@02d095ba.0 0.05252120 0.04180259 -0.01428254
## NSFTV87@64fa0112.0 0.05252120 -0.03758906 -0.01428254
## NSFTV96@32a6808e.0 0.05252120 -0.03758906 -0.01428254
## NSFTV100@1f10be3d.0 -0.03012755 0.04180259 -0.01428254
## NSFTV114@eac19fb8.0 -0.03012755 0.04180259 0.11516424
## NSFTV140@85551f9c.0 -0.03012755 0.04180259 -0.01428254フーリエ記述子から、輪郭を再構成するための関数を定義しておきます。これは、前々回の講義で用いた関数です。
ef2coord <- function(ef, theta = seq(0, 2*pi, 0.01)) {
x <- 0; y <- 0
for(i in 1:(length(ef)/4)) {
x <- x + ef[i * 4 - 3] * cos(i * theta) + ef[i * 4 - 2] * sin(i * theta)
y <- y + ef[i * 4 - 1] * cos(i * theta) + ef[i * 4 ] * sin(i * theta)
}
coord <- list(x = x, y = y)
return(coord)
}では、1SNPずつ図示していきましょう。まずは、第3染色体のSNPです。
z.u <- max(g[,1])
z.l <- min(g[,1])
x.u <- z.u * ef.pca$rotation[,1]
x.l <- z.l * ef.pca$rotation[,1]
ef.u <- c(1, 0, 0, ef.pca$center + x.u)
ef.l <- c(1, 0, 0, ef.pca$center + x.l)
coord.u <- ef2coord(ef.u)
coord.l <- ef2coord(ef.l)
plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.l, col = "orange")
次に、第5染色体のSNPです。
z.u <- max(g[,2])
z.l <- min(g[,2])
x.u <- z.u * ef.pca$rotation[,1]
x.l <- z.l * ef.pca$rotation[,1]
ef.u <- c(1, 0, 0, ef.pca$center + x.u)
ef.l <- c(1, 0, 0, ef.pca$center + x.l)
coord.u <- ef2coord(ef.u)
coord.l <- ef2coord(ef.l)
plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.l, col = "orange")
最後に第7染色体のSNPです。
z.u <- max(g[,3])
z.l <- min(g[,3])
x.u <- z.u * ef.pca$rotation[,1]
x.l <- z.l * ef.pca$rotation[,1]
ef.u <- c(1, 0, 0, ef.pca$center + x.u)
ef.l <- c(1, 0, 0, ef.pca$center + x.l)
coord.u <- ef2coord(ef.u)
coord.l <- ef2coord(ef.l)
plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.l, col = "orange")
第7染色体の変異が思いのほか大きいことがわかります。推定された遺伝子型の効果をヒストグラムを描いて確認してみます。
hist(g[,3])
対立遺伝子（アリル）頻度が偏っているようです。この7番染色体のSNPが真のアソシエーションを捉えているのだとすると、アリル頻度の偏りにより第3, 5番染色体のSNPよりも検出が難しかったとも考えられます。ただ、アリル頻度が低いために、その推定値も過度に当てはめられた値になっている可能性もあります。より詳細な検証が必要です。
形状の第2主成分のGWAS
なお、前々回の講義で、微細な形の変異を評価していることがわかった楕円フーリエ記述子の第2主成分についても、関連するSNPがあるかどうか、GWASを行って確認してみましょう。
どのような形状変異であったか、主成分分析の結果から改めてかたちの変動を逆変換で求めて、視覚化してみましょう。
z.u <- 0 + 2 * ef.pca$sdev[2] # 標準偏差なので平方根。平均は0。
z.l <- 0 - 2 * ef.pca$sdev[2]
x.u <- z.u * ef.pca$rotation[,2]
x.l <- z.l * ef.pca$rotation[,2]
ef.u <- c(1, 0, 0, ef.pca$center + x.u)
ef.m <- c(1, 0, 0, ef.pca$center)
ef.l <- c(1, 0, 0, ef.pca$center + x.l)
coord.u <- ef2coord(ef.u)
coord.m <- ef2coord(ef.m)
coord.l <- ef2coord(ef.l)
plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.m, col = "gray")
lines(coord.l, col = "orange")
やはり、かなり微細な変異です。このような変異に関連するSNPは検出されるでしょうか？解析を進めていきましょう。まずは、第2主成分の得点を\(y\)として取り出します。また、ゲノムデータbm.wpに合わせて並べ替えておきます。
y <- ef.pca$x[,2] # select seed length with ratio as the target trait
names(y) <- rownames(ef.pca$x)
y <- y[bm.wp@ped$id]早速、GWASを行ってみましょう。
gwa <- association.test(bm.wp, Y = y, method = "lmm",
response = "quantitative",
test = "lrt", eigenK = eigen(grm), p = 4)偽発見率（FDR）を計算し、FDRが5%以下のSNPに色付けしてマンハッタンプロットを描きます。
fdr <- p.adjust(gwa$p, method = "BH")
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)
第3, 6染色体上に2つほど有意なSNPがあるようですが、\(-\log10(p)\)の値は、第1主成分のピークに比べると顕著ではないことがわかります。ただし、この結果は、非常に微細に見える形状も、何らか遺伝的支配を受けていることを示唆しています。