生物測定担当分：ゲノムワイドアソシエーション研究（GWAS）

ゲノムワイドアソシエーション研究（GWAS）

利用するRのパッケージ

今回の学生実験ではgastonとよばれるgwas用のパッケージを使用する。なお、LDmap.RにはLDヒートマップ作製に必要な関数が定義されている

require(here)
require(gaston)
require(LDheatmap)          
source(here("source", "GWAScomp.R"))
source(here("source", "LDmap.R"))
source(here("source", "GWASpeak.R"))

利用するゲノムデータの読み込み

McCouchら(2016; Nature Communication 7:10532)が用いたHigh Density Rice Array (HDRA)を用いて下のタイピングされたSNP遺伝子型データのvcfファイルを読み込む。読み込みにはgastonのread.vcf関数を用いる。

bm <- read.vcf(
  here("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

## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are set

イネ遺伝資源悪セッション1568系統についての38769SNPsの遺伝子型データが含まれている。bmにはSNP遺伝子型が数値としてコードされたものが含まれている

as.matrix(bm[1:10,1:6])

##                       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.
## IRGC121316@c88dcbba.0            0            2
## 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            2
## IRGC121314@6c2404e4.0            0            0

SNP遺伝子型が、“0”と“2”でコードされていることが分かる。また、4列目の上から８番目のように、“1”でコードされているSNP遺伝子型もある。“0”, “1”, “2”は、それぞれ、A1/A1、A1/A2、A2/A2の遺伝子型を表し、“0”と“2”はホモ接合体、“1”はヘテロ接合体である。なお、栽培イネは自殖性のため、ヘテロ接合の頻度は低い。

なお、A1やA2の対立遺伝子が、それぞれ、どの塩基に対応するのかは、SNP情報で確認できる。

head(bm@snps)

##   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
##   N2.f NAs.f callrate        maf          hz
## 1   NA    NA        1 0.14126276 0.004464286
## 2   NA    NA        1 0.05293367 0.005102041
## 3   NA    NA        1 0.13265306 0.015306122
## 4   NA    NA        1 0.26179847 0.005739796
## 5   NA    NA        1 0.14030612 0.002551020
## 6   NA    NA        1 0.34917092 0.009566327

bm@snpsには、染色体番号（chr）、マーカー名（id）、物理位置（pos）、A1とA2の対立遺伝子（A1, A2)、3つのSNP遺伝子型（A1/A1, A1/A2, A2/A2）の観察数（N0, N1, N2）、非欠測率（callrate）、マイナーアリル頻度（maf）などの情報が含まれている。

利用する表現型データの読み込み

表現型データとしてはZhaoら(2011; Nature Communication 2:467)がイネ遺伝資源のGWASに用いた表現型データを使用する。

pheno <- read.csv(
  here("data", "RiceDiversityPheno4GWASGS.csv"), row.name = 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形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。

形質とマーカー遺伝子型データの抽出

36ある形質のうち、まずは種子長(Seed.length)のデータを用いて解析をする。一つの変数を取り出すと行名であった系統名が消えてしまうため、names関数を用いて再定義する必要がある。その上でomit関数を使ってNA(欠損値)を取り除く。

y <- pheno[ , "Seed.length"]
names(y) <- rownames(pheno)
y <- na.omit(y)

また、今回使用するマーカー遺伝子型データに含まれる系統数は表現型データの系統数を上回る。表現型データがない個体については解析ができないため、表現型データを持つ系統のマーカー遺伝子型データをここで抜き出す。

bm.wp <- bm[bm@ped$id %in% names(y),]
bm.wp

## A bed.matrix with 361 individuals and 38769 markers.
## snps stats are set
##   There is one monomorphic SNP
## ped stats are set

1568個体あった個体が361個体分に減っているのがわかる。

また、抽出しただけでは、表現型の並びとマーカー遺伝子型の並びが一致しないため以下のコードで\(y\)を並び替える。

y <- y[bm.wp@ped$id]

分散分析による解析

ここでは原因遺伝子を特定する方法のひつとして分散分析による手法を説明する。 ある表現型とマーカー情報が得られているとし、それぞれのマーカーは3種類(A1/A1,A1/A2,A2/A2,すなわち、0, 1, 2)の遺伝子型を持っているとする。このマーカー近傍の遺伝子が表現型に影響をしている場合、遺伝子型ごとにその表現型の平均値は異なると考えられる。これを分散分析で検定することで遺伝子の効果を調べることが可能である。

まずはじめに1番目のマーカーについて、３つの遺伝子型間の表現型値の違いを箱ひげ図を描いて確認する。遺伝子型を数値としてではなく、因子（factor）として質的に扱っていることに注意する。

geno.group <- factor(as.matrix(bm.wp[, 1])) #genotypes of the 1st marker
boxplot(y ~ geno.group)

つぎに、３つのグループ間の差が有意か検定を行って確認する。

model <- lm(y ~ geno.group)
anova(model)

## Analysis of Variance Table
## 
## Response: y
##             Df Sum Sq Mean Sq F value Pr(>F)
## geno.group   2   2.69  1.3426  1.4937 0.2259
## Residuals  358 321.81  0.8989

\(p\)値は5%よりも大きく、このマーカーの遺伝子型間に表現型\(y\)に差があるわけではないことがわかる。これをすべてのマーカーに対して行うことで、注目するマーカーの遺伝子型間で有意差があるマーカーを見つけ出すことができる。

さて、ここでマーカーの番号（marker.id）を指定すれば、その遺伝子型を要因として分散分析を行い、その\(p\)値が出力される関数geno.anovaを定義する。

geno.anova <- function(y, bm.wp, marker.id) {
  geno.group <- factor(as.matrix(bm.wp[, marker.id])) #genotype for anova
  model <- lm(y ~ geno.group)
  res <- anova(model)
  return(res$`Pr(>F)`[1])
}
geno.anova(y, bm.wp, 5)

## [1] 0.2777978

これを使って1〜400番目のマーカーの\(p\)値を調べてみる。ここではsapply関数を用いる。1つ目の引数の数字が2つ目の引数の関数に適用されるというものである。

p.vector <- sapply(1:400, geno.anova, y = y, bm.wp = bm.wp)

結果はp値で返される。p値は小さいほど有意であることを表すが、小さな値にみられる違いを確認するのは難しいため-log10(p)として変換してみる。さらにその変換した値をplot関数で視覚的に表してみる。

mlogp <- -log10(p.vector)
plot(mlogp)

こうすると、200番目付近のマーカーに有意なピークがあることが視覚的に見て取れる。このマーカーはwhich.max関数を使って調べることができる。

max(mlogp)

## [1] 6.678641

which.max(mlogp)

## [1] 221

221番目のマーカーが-log10(p)の値が6.68となっていることがわかる。

表現型を221番のマーカーの遺伝子型ごとに分類してプロットすると、確かに遺伝子型間での表現型の違いが認められる。

geno.group <- factor(as.matrix(bm.wp[, which.max(mlogp)]))
boxplot(y ~ geno.group)

さて、次のグラフは先の箱ひげ図をプロットで表したものである。今までは遺伝子型0,1,2を質的に捉えていたが、これを連続的に捉え、表現型に対して遺伝子型を回帰した結果も載せている。

geno.x <-as.matrix(bm.wp[, which.max(mlogp)])
plot(geno.x, y)
model <- lm(y ~ geno.x)
abline(model)

さらに、回帰モデルについて分散分析でその有意性を検定する。

anova(model)

## Analysis of Variance Table
## 
## Response: y
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## geno.x      1  26.254 26.2542  31.603 3.804e-08 ***
## Residuals 359 298.239  0.8307                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

連続データとして捉えた場合も分散分析の結果は有意である。「有意差がある=傾きが0でない」ととらえることができる。

単回帰分析に基づくGWAS (GWAS navieモデル)

これまでマーカー遺伝子型のスコア（0,1,2）を質的に捉え検定を行ってきたが、先に述べたようにスコアを0,1,2を連続的な値として扱い、形質の表現型値を回帰して、傾きの傾きの有意性検定を行うことでGWASを行うことができる。

最も単純なモデルは以下のように表される。ここで\(y_{i}\)はi番目の個体の表現型値、\(u\)は切片,\(\beta_{j}\)はj番目のマーカーに対する傾き,\(x_{ij}\)はマーカー遺伝子型のスコア(0,1,2),\(e_{i}\)は誤差を表している。

\[ y_{i}=u+\beta_{j}x_{ij}+e_{i} \] gastonのassotiation.test関数を用いると、すべてのマーカーに対して回帰を行った結果を得ることができる。

なお、ここでbm.wpは解析に用いるSNPデータ、Y = yは、解析に用いる表現型データを指定している。method = “lm”は（母数効果だけを含む）線形モデルを、method = “lmm” は、（母数効果と変量効果を含む）線形混合モデルをあてはめることを意味する。通常は“lmm”を指定して解析するが、ここでは、上のシンプルなモデルを当てはめるため“lm”とする。response = “quantitative”は、対象としている形質が量的形質であることを意味する。また、test = は、効果の検定法を表し、尤度比検定（Likelihood Ratio Test: LRT)“lrt”の他に、“wald”, “score”を選ぶことができる。ここでは、“wald”とする。

gwa<- association.test(bm.wp, Y = y[bm.wp@ped$id], 
                       method = "lm", response = "quantitative", 
                       test = "wald")
manhattan(gwa, pch = 20)

qqplot.pvalues(gwa, pch=20)

1つ目はマンハッタンプロットとよばれ、回帰係数の有意性を表す\(p\)値を\(-log_{10}(p)\)として変換して図示したものだる。横軸はマーカーの位置を表しており、ピークの付近に原因遺伝子が座乗する可能性がある。2つ目はqqplotとよばれ、全ての検定において帰無仮説が正しいときに期待される\(-log_{10}(p)\)の分布から、どれだけ逸脱しているかを読み取る。両者が一致する場合は、点が概ね直線上に散布されるが、帰無仮説からの逸脱があると直線から外れたところに散布される。

上の結果では点が直線から離れたところに散布されており、期待される-log10(p)の値よりも大きな値が得られていることがわかる。一部のSNPでは帰無仮説が正しくなく、本当に-log10(p)が大きくなるかもしれないが、多くのSNPで-log10(p)の値が過大になるのは、解析集団の背景にある分集団構造をモデルにうまく取り込めていないことに由来する。

グループワーク①: 分集団構造による影響の形質による違い (20分程度)

課題1: 形質を変えてqqplotを描き、分集団構造の影響（-log10(p)の値が過大になる程度）を形質間で比較せよ。分集団構造による影響の程度が大きい、または、小さい形質を調べて、どのような形質で特に影響が大きいか考察せよ。

形質の名前を番号とともに表示する。

paste(1:ncol(pheno), colnames(pheno))

##  [1] "1 Flowering.time.at.Arkansas"        "2 Flowering.time.at.Faridpur"       
##  [3] "3 Flowering.time.at.Aberdeen"        "4 FT.ratio.of.Arkansas.Aberdeen"    
##  [5] "5 FT.ratio.of.Faridpur.Aberdeen"     "6 Culm.habit"                       
##  [7] "7 Leaf.pubescence"                   "8 Flag.leaf.length"                 
##  [9] "9 Flag.leaf.width"                   "10 Awn.presence"                    
## [11] "11 Panicle.number.per.plant"         "12 Plant.height"                    
## [13] "13 Panicle.length"                   "14 Primary.panicle.branch.number"   
## [15] "15 Seed.number.per.panicle"          "16 Florets.per.panicle"             
## [17] "17 Panicle.fertility"                "18 Seed.length"                     
## [19] "19 Seed.width"                       "20 Seed.volume"                     
## [21] "21 Seed.surface.area"                "22 Brown.rice.seed.length"          
## [23] "23 Brown.rice.seed.width"            "24 Brown.rice.surface.area"         
## [25] "25 Brown.rice.volume"                "26 Seed.length.width.ratio"         
## [27] "27 Brown.rice.length.width.ratio"    "28 Seed.color"                      
## [29] "29 Pericarp.color"                   "30 Straighthead.suseptability"      
## [31] "31 Blast.resistance"                 "32 Amylose.content"                 
## [33] "33 Alkali.spreading.value"           "34 Protein.content"                 
## [35] "35 Year07Flowering.time.at.Arkansas" "36 Year06Flowering.time.at.Arkansas"

次の関数は、解析する形質の番号をtrait.idで指定するとqqplotを表示する。

qqplot.trait <- function(bm, pheno, trait.id) {
  y <- pheno[ , trait.id]
  names(y) <- rownames(pheno)
  y <- na.omit(y)
  bm.wp <- bm[bm@ped$id %in% names(y),]
  y <- y[bm.wp@ped$id]

  gwa<- association.test(bm.wp, Y = y,
                       method = "lm", response = "quantitative", 
                       test = "wald")
  qqplot.pvalues(gwa, pch = 20, main = colnames(pheno)[trait.id])
}

上に示した形質名から対象とする形質を選び、その番号をtrait.idに設定する。

trait.id <- 1 # 形質の番号
qqplot.trait(bm, pheno, trait.id)

分集団構造を抑えるための統計モデル

　このように、単回帰に基づくGWASでは、分集団構造により-log10(p)が過大になり、場合によっては偽陽性が生じる。分集団構造による影響を抑えるために、これまでに様々なモデルが提案された。ここでは、マクロな分集団構造(sub-population structure)を考慮したKモデルと、ミクロな分集団構造（家系構造）を考慮したQモデル、それらを組み合わせたQ+Kモデルを紹介する。その前に、\(p\)値、あるいは、-log10(p)の値の閾値の決め方について説明する。

有意水準

Manhattanプロットでは-log10(p)のピークが得られるが、それが有意であるかどうかを判断するためには何らかの基準が必要である。一般的な検定では5％水準などが使われるが、多くのマーカーを同時に検定するGWASでは検定の多重性が問題になる。例えば10,000個のマーカーを検定するとする。5％水準で検定した場合、すべてのマーカーについて帰無仮説が正しい場合でも500個のマーカーで\(p\)値が5%を下回ると期待される。多重検定のもとで適切に\(p\)値の基準を設定するために、さまざまな方法が提案されている。ここでは、GWASでよく用いられるBenjamini-Hochberg（BH）法を紹介する。

まずは各マーカーを通常の検定のように5%水準で検定し有意なマーカーを黄緑で表す。

gwa<- association.test(bm.wp, Y = y, 
                       method = "lm", response = "quantitative", 
                       test = "wald")
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[gwa$p < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)

このようにほとんどのマーカーが有意になる。これは分集団構造による影響で-log10(p)の値が過大になったためである。この結果に対してBH法を適用してみる。BH法では、偽発見率（検出された有意な結果のうちの偽陽性の割合）が指定した水準以下になるように基準を設定する手法である。p.adjust関数を使うことで、各マーカーに対する偽発見率（false discovery rate: FDR）を計算できる。

gwa<- association.test(bm.wp, Y = y, 
                       method = "lm", response = "quantitative", 
                       test = "wald")
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)

なお、分集団構造による影響により、有意水準を補正しても偽陽性を抑えられない。分集団構造による影響を抑えるには、有意水準の補正ではなく、統計モデルそのものを改良する必要がある。

ミクロな分集団構造（家系構造）を考慮したモデル (Kモデル)

　品種群のようなマクロな分集団構造に加え、家系構造のようなミクロな構造も偽陽性の原因となる。ゲノムワイドマーカーで見られる類似度をもとに、こうしたミクロな家系構造による影響も考慮に入れるのがKモデルであり、先のNaiveモデルに新たに\(\alpha_{i}\)という項がここに加わる。

\[ y_{i}=u+\beta_{j}x_{ij}+\alpha_{i}+e_{i} \] この新たに加わった\(\alpha_{i}\)はi番目の個体における遺伝効果を表している。この効果は遺伝的に近い個体同士は近い値をもつとして、\(\alpha_{i}\)に対して次のような仮定をおく。

\[ \alpha \sim N(0, K\sigma^2) \] ここで\(\alpha={\alpha_{1},\alpha_{2},...\alpha_{n}}\)であり、\(K\)はゲノム関係行列と呼ばれる。ゲノム関係行列におけるi行j列の要素はi番目とj番目の個体の遺伝的な距離を表しており、これによりミクロな分集団構造である家系構造をモデルに取り込むことができる。

このモデルで検定するためにはゲノム関係行列を計算する必要があるので、まずこれを計算する。

grm <- GRM(bm.wp)

次にこのゲノム関係行列を使って検定を行うが、先ほど動かしたコードに対して、method = “lmm”とし、eigenK=eigen(grm)を新しく加えればよい。lmmは線形混合モデルを表しており、eigen関数を用いてGWASでの計算に内部的に利用されているゲノム関係行列の固有値分解の結果を与えている。

gwa <- association.test(bm.wp, Y = y, 
                        method = "lmm", response = "quantitative", 
                       test = "wald", eigenK=eigen(grm))
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)

qqplot.pvalues(gwa,pch=20)

navieモデルに比べ、分集団構造による影響を抑え込むことに成功し、3番染色体上にきれいなピークが現れた。

分集団による影響を効果として取り込んだモデル (Qモデル)

　次に、マクロな分集団構造を意識したQモデルを紹介する。分集団というのはsub-populationと英語で書くことからもわかるように種や集団内に存在する二次的な集団のことである。たとえばイネではIndica(インド型)やJaponica(日本型)などがある。今回解析の対象である種子の長さが両集団で異なることは有名であるが、そのほかにも分集団間で大きく異なる形質があり、この存在がGWASの検出力を低下させる原因になっている。従ってこの分集団を適切にモデルに考慮することが必要である。

Qモデルには先ほどのnavieモデルに対して新たに\(\sum_{k=1}^{K} v_{k}q_{ik}\)という項が加わる。ここで、\(v_{k}\)は回帰係数、\(q_{ik}\)はk番目の分集団の影響の割合を表している。分集団構造の次元（正確には主成分の数）\(K\)はある程度自由に変えられる。 \[ y_{i}=u+\beta_{j}x_{ij}+\sum_{k=1}^{K} v_{k}q_{ik}+ e_{i} \] このモデルでGWASを行うには、method=“lm”とするとともに、p=5という分集団の数を指定するオプションを設定する。5は主成分数であり1や3など、自由に変えることができる。

gwa<- association.test(bm.wp, Y = y, 
                       method = "lm", response = "quantitative", 
                       test = "wald", eigenK = eigen(grm), 
                       p = 5)
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)

qqplot.pvalues(gwa,pch=20)

このデータでは、Naiveモデルよりは直線からの逸脱は小さいが、Kモデルほどの改善が見られない。

両方を含めたモデル (Q+Kモデル)

QモデルとKモデル両方を組み合わせることも可能で、これはQ+Kモデルとして知られている。これを実行するにはKモデルに対してlmをlmmと置き換えるだけでよい。

gwa<- association.test(bm.wp, Y = y, 
                       method = "lmm", response = "quantitative", 
                       test = "wald", eigenK = eigen(grm), 
                       p = 5)
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)

qqplot.pvalues(gwa,pch=20)

Kモデルの結果と同様、3番染色体に顕著に有意なピークがみられる。

さて、今回の形質(種子の長さ)ではKモデルにより分集団による影響が抑えられることがわかったが、これは対象とする形質により変わる可能性がある。形質とモデルを変えて、その違いを確認してみよう。

グループワーク② 形質とモデルと分集団構造による影響

解析の対象とする形質と解析に用いるモデル、また、QモデルとQ+Kモデルで設定する主成分の数を変化させて、分集団構造による影響の抑制程度の違いを確認して、自由に考察せよ。

trait.id <- 18 # 解析する形質の番号
model <- "N" # 用いるモデル N:Navie, K:K, Q:Q, QK:Q+K
n.pca <- 1 # Q, Q+Kを指定した場合は考慮する主成分数を指定。
           # それ以外では無視される

GWAS_comp(bm, pheno, trait.id, model, n.pca, gwa.res = FALSE)

原因遺伝子の調査

最後に原因遺伝子を調べる方法を紹介する。先ほどと同様、種子長を対象とし、Q+Kモデルの結果を用いる。先の関数と同じものをつかうが、gwa.res = TRUEとすると結果の出力もできる。

trait.id <- 18 # 種子の長さ
model <- "QK" # モデルの指定 N:Navie, K:K, Q:Q, QK:Q+K
n.pca <- 5 # Q, Q+Kを指定した場合は考慮する主成分数を指定

gwa <- GWAS_comp(bm, pheno, trait.id, model, n.pca, gwa.res = TRUE)

まずはピークに対応するマーカーを調べる。

fdr <- p.adjust(gwa$p, method = "BH")
gwa[fdr < 0.05, ]

##       chr      pos              id   A1   A2    freqA2        h2       beta
## 11250   3 15991434 SNP-3.15990080.    G    A 0.6011080 0.9900000  0.3086103
## 11256   3 16056312 SNP-3.16054958.    T    C 0.5983380 0.9900000  0.3454859
## 11279   3 16202379 SNP-3.16201026.    A    G 0.4584488 0.9900000  0.2741251
## 11282   3 16210385 SNP-3.16209032.    A    G 0.1980609 0.9900000  0.2969402
## 11293   3 16274770 SNP-3.16273416.    C    A 0.1925208 0.9900000  0.2978098
## 11304   3 16310295 SNP-3.16308941.    G    T 0.1759003 0.9900000  0.3370494
## 11308   3 16325371 SNP-3.16324017.    C    T 0.3033241 0.9900000  0.2442254
## 11354   3 16664927 SNP-3.16663571.    G    A 0.1121884 0.9900000  0.3838574
## 11357   3 16688126 SNP-3.16686770.    G    A 0.3130194 0.9891084 -0.2966389
## 11360   3 16711055 SNP-3.16709700.    G    T 0.1149584 0.9900000  0.3674665
## 11361   3 16719368 SNP-3.16718013.    G    A 0.2562327 0.9878470 -0.4900465
## 11362   3 16733441 SNP-3.16732086.    G    T 0.3725762 0.9796198  0.6112241
## 11387   3 16889939 SNP-3.16888812.    G    A 0.2229917 0.9897760  0.3738233
## 11423   3 17348645 SNP-3.17347454.    A    G 0.7368421 0.9887392  0.3473225
## 11433   3 17453008 SNP-3.17451789.    C    T 0.3296399 0.9892100  0.2803271
## 17117   5  5516194  SNP-5.5516131.    G    A 0.6426593 0.9900000  0.2949359
## NA     NA       NA            <NA> <NA> <NA>        NA        NA         NA
##               sd            p
## 11250 0.06883634 7.351585e-06
## 11256 0.07127162 1.250564e-06
## 11279 0.05519982 6.832974e-07
## 11282 0.06269481 2.176772e-06
## 11293 0.06328219 2.525496e-06
## 11304 0.06460729 1.819479e-07
## 11308 0.05705953 1.867372e-05
## 11354 0.07105893 6.591717e-08
## 11357 0.06220005 1.850415e-06
## 11360 0.07015447 1.623585e-07
## 11361 0.08295382 3.474180e-09
## 11362 0.05491850 8.997947e-29
## 11387 0.06972127 8.245066e-08
## 11423 0.06621642 1.560605e-07
## 11433 0.06358282 1.039214e-05
## 17117 0.05777390 3.307889e-07
## NA            NA           NA

posは物理位置を表しており、idはマーカーのidを表している。これらが原因遺伝子を調べるための手掛かりとなる。

なお、緑のピーク(有意なピーク)がない場合でも次のようにコードを書くことで、manhattanプロットの中で値の高いマーカーを調べることができる。

gwa[order(gwa$p)[1:10], ]

##       chr      pos              id A1 A2    freqA2        h2       beta
## 11362   3 16733441 SNP-3.16732086.  G  T 0.3725762 0.9796198  0.6112241
## 11361   3 16719368 SNP-3.16718013.  G  A 0.2562327 0.9878470 -0.4900465
## 11354   3 16664927 SNP-3.16663571.  G  A 0.1121884 0.9900000  0.3838574
## 11387   3 16889939 SNP-3.16888812.  G  A 0.2229917 0.9897760  0.3738233
## 11423   3 17348645 SNP-3.17347454.  A  G 0.7368421 0.9887392  0.3473225
## 11360   3 16711055 SNP-3.16709700.  G  T 0.1149584 0.9900000  0.3674665
## 11304   3 16310295 SNP-3.16308941.  G  T 0.1759003 0.9900000  0.3370494
## 17117   5  5516194  SNP-5.5516131.  G  A 0.6426593 0.9900000  0.2949359
## 11279   3 16202379 SNP-3.16201026.  A  G 0.4584488 0.9900000  0.2741251
## 11256   3 16056312 SNP-3.16054958.  T  C 0.5983380 0.9900000  0.3454859
##               sd            p
## 11362 0.05491850 8.997947e-29
## 11361 0.08295382 3.474180e-09
## 11354 0.07105893 6.591717e-08
## 11387 0.06972127 8.245066e-08
## 11423 0.06621642 1.560605e-07
## 11360 0.07015447 1.623585e-07
## 11304 0.06460729 1.819479e-07
## 17117 0.05777390 3.307889e-07
## 11279 0.05519982 6.832974e-07
## 11256 0.07127162 1.250564e-06

このコードは、p値が低い順にマーカーを表示するコードである。今回はp値が最も小さかったSNP-3.16732086.を対象に解析を行う

LDプロットの作成

　次に先ほど求めたidをもとにLDプロットと呼ばれるプロットを作成する。LDプロットのLDとは連鎖不平衡(Linkage Disequilibrium)の略であるが、このLDとは直観的にはマーカー間の非独立的な関係である。同一染色体上に近接しているマーカー間では、それらは連鎖していることで親から子にともに伝えられる可能性が高い。この程度を表したのがLDである。

GWASでピークが見られた場合、そのピークに近接した位置に原因遺伝子がある可能性が高いが、“近接”の範囲がどの程度なのかが明らかではない。そこで、ピークを中心として近接しているSNPとのLDのパターンを図示することで、どの範囲までを探索すればよいかを確認する。

LD.map関数はLDプロットを書くための関数であり、phenotypeには表現型を、id.makerにはピークが検出されたマーカーのidを指定する。プロットではidで指定されたマーカーの左右\(n\)個分がプロットされるが、この数を指定しているのがn.makerである。(phenoは本来必要ないが計算の都合上考慮している)。

LD.map(bm, y, gwa,
       id.marker = "SNP-3.16732086.",
       n.marker = 50)

一つ目は通常のLDプロット、二つ目はそこにGWASの結果を対応付けて表示したものである。ピークと強く相関するマーカーが多く隣接する場合はそれらを含めた領域を探索する必要があるが、今回は比較的マーカーの独立性が高く、ピークに対応する物理位置のごく周辺だけ(3番染色体16733441kb付近)を探索すればよさそうである。

rap-dbによる遺伝子探索

最後にRAP-dbと呼ばれるデータベースを使って3番染色体16733441kb付近に原因となる遺伝子の候補がないかを確認する。RAP-dbはRice Annotation Projectと呼ばれるプロジェクトにより明らかになった遺伝子のデータベースでありイネに関する膨大なゲノム情報がまとめられている。今回はその中でもAgriGeneリストと呼ばれるサイトを用いて物理位置から原因遺伝子候補を探してみる。

AgriGeneリスト https://rapdb.dna.affrc.go.jp/agri_genes/agri_gene_list.html このサイトを開くとデータベースが現れる。ここの3列目がpositionとなっており、物理位置に対応する。ここから3番染色体 16,733,441 kb 付近を探すと、Os03g0407400というIDの領域が見つかった。これは粒形に関する遺伝子をコードしている領域で、GWASにより原因遺伝子を検出できることが示された。

グループワーク③ 遺伝子探索

特に顕著なピークがみられた形質について、顕著なピークに対応する原因遺伝子候補を以下の流れに沿って調べよ。

ピークのマーカーIDを検索

次のコードでは、形質のIDとモデル、主成分数(Q/Q+Kモデルのみ)、FDRの閾値を指定すると、閾値を超えたピークのうち、各染色体で最も高いピークの情報を出力できる。

trait.id <- 18 # 解析対象とする形質の番号。
method <- "QK" # モデル N:Navie, K:K, Q:Q, Q+K:Q+K
n.pca <- 5 # Q, Q+Kを指定した場合は考慮するpcaの数を指定

gwa <- GWAS_comp(bm, pheno, trait.id, model, n.pca, gwa.res = TRUE)

GWAS_peak(gwa)

##       chr      pos              id A1 A2    freqA2        h2      beta
## 11362   3 16733441 SNP-3.16732086.  G  T 0.3725762 0.9796202 0.6112237
## 17117   5  5516194  SNP-5.5516131.  G  A 0.6426593 0.9900000 0.2949360
##               sd            p
## 11362 0.05491855 8.999661e-29
## 17117 0.05777391 3.307892e-07

Step2 LDプロットの表示

次の先ほどの結果からマーカーIDを選び次のコードを実行する。このコードによりLDプロットが出力される。

id.marker <- "SNP-3.16732086." # 検索するマーカー
n.marker <- 50 # 考慮する左右のマーカー数
LD.map(bm, y, gwa,
       id.marker = id.marker,
       n.marker = 50)

Step3 AgriGeneでの検索

2つ前の結果から物理位置と染色体番号を割り出しAgriGeneリストでその付近の遺伝子を調べる。見つかったらその遺伝子がピークと関連ありそうかをLDプロットをもとに考察する(その遺伝子がピークのマーカーと強く相関していそうかという観点が重要)。