本実験では、Rのパッケージagridatに含まれているオオムギの多環境試験データを例として解析を行う。
なお、表現型データについては以下からダウンロードすることもできる。 http://wheat.pw.usda.gov/ggpages/SxM/phenotypes.html また、遺伝子型データは、以前は以下からダウンロードできたが、現在はできないようである。 http://www.genenetwork.org/genotypes/SXM.geno
このデータセットには、オオムギ品種Steptoe×Morexの交配から得られた150の倍加半数体（doubled haploid、DH）系統のデータが含まれている。 遺伝子型データには7つの染色体全体に分布する218のマーカーの遺伝子型データが含まれる。 形質データは、最大16環境で計測された表現型データが含まれている。
まずは、必要なパッケージを読み込む。
require(qtl) # R package for QTL mapping## Loading required package: qtlrequire(agridat) # QTL data## Loading required package: agridatrequire(tidyverse) # for rearranging data structure ## Loading required package: tidyverse## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()tidyverseは、データの整形に用いると便利なパッケージである。
次に、agridatからデータを読み込む。なお、このプログラムの最後に、csvファイルからデータを読み込む方法も紹介する。
まずは、マーカー遺伝子型を読み込む。
geno <- steptoe.morex.geno
geno## This is an object of class "cross".
## It is too complex to print, so we provide just this summary.
## Doubled haploids
##
## No. individuals: 150
##
## No. phenotypes: 1
## Percent phenotyped: 100
##
## No. chromosomes: 7
## Autosomes: 1 2 3 4 5 6 7
##
## Total markers: 223
## No. markers: 37 37 31 33 29 22 34
## Percent genotyped: 96
## Genotypes (%): AA:50.9 BB:49.1次に、表現型データを読み込む。
pheno <- steptoe.morex.pheno
head(pheno)## gen env amylase diapow hddate lodging malt height protein yield
## 1 Steptoe MN92 22.7 46 149.5 NA 73.6 84.5 10.5 5.5315
## 2 Steptoe MTi92 30.1 72 178.0 10 76.5 NA 11.2 8.6403
## 3 Steptoe MTd92 26.7 78 165.0 15 74.5 75.5 13.4 5.8990
## 4 Steptoe ID91 26.2 74 179.0 NA 74.1 111.0 12.1 8.6290
## 5 Steptoe OR91 19.6 62 191.0 NA 71.5 90.0 11.7 5.3440
## 6 Steptoe WA91 23.6 54 181.0 NA 73.8 112.0 10.0 6.2700この表現型データには、全ての形質のデータが含まれている。ここから、解析対象とする形質のデータを抜き取って解析をする。
まずは、表現型データのサマリを確認する。
summary(pheno)## gen env amylase diapow hddate
## Morex : 16 ID91 : 152 Min. :14.90 Min. : 35.00 Min. :143.0
## SM1 : 16 ID92 : 152 1st Qu.:25.62 1st Qu.: 70.00 1st Qu.:174.5
## SM10 : 16 MA92 : 152 Median :28.50 Median : 84.00 Median :183.5
## SM103 : 16 MN92 : 152 Mean :29.04 Mean : 87.25 Mean :181.2
## SM104 : 16 MTd91 : 152 3rd Qu.:32.00 3rd Qu.:101.00 3rd Qu.:191.0
## SM105 : 16 MTd92 : 152 Max. :49.30 Max. :229.00 Max. :217.0
## (Other):2336 (Other):1520 NA's :1066 NA's :1066
## lodging malt height protein
## Min. : 0.00 Min. :69.00 Min. : 34.00 Min. : 8.00
## 1st Qu.: 15.00 1st Qu.:73.30 1st Qu.: 82.50 1st Qu.:11.90
## Median : 35.00 Median :74.50 Median : 95.00 Median :13.00
## Mean : 37.13 Mean :74.72 Mean : 94.95 Mean :12.92
## 3rd Qu.: 55.00 3rd Qu.:75.90 3rd Qu.:109.00 3rd Qu.:14.00
## Max. :100.00 Max. :83.00 Max. :151.00 Max. :17.50
## NA's :1521 NA's :1067 NA's :5 NA's :1067
## yield
## Min. : 1.390
## 1st Qu.: 4.092
## Median : 5.272
## Mean : 5.293
## 3rd Qu.: 6.338
## Max. :11.526
## “amylose”以降の列が形質データである。
ここでは、“yield”を解析対象とする。
target <- "yield"
td <- pheno %>% select(gen, env, all_of(target))
head(td)## gen env yield
## 1 Steptoe MN92 5.5315
## 2 Steptoe MTi92 8.6403
## 3 Steptoe MTd92 5.8990
## 4 Steptoe ID91 8.6290
## 5 Steptoe OR91 5.3440
## 6 Steptoe WA91 6.2700なお、このデータは、復数の環境のデータが1列で表されている。これを、環境ごとに1列になるように整形する。
td.w <- td %>% pivot_wider(id_cols = gen,
names_from = env,
values_from = all_of(target))
head(td.w)## # A tibble: 6 × 17
## gen MN92 MTi92 MTd92 ID91 OR91 WA91 MTi91 MTd91 NY92 ON92 WA92 MA92
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Stept… 5.53 8.64 5.90 8.63 5.34 6.27 4.10 7.07 6.05 3.70 4.10 6.43
## 2 Morex 4.13 4.18 3.71 7.95 8.33 5.97 1.74 4.67 5.50 3.14 3.26 7.51
## 3 SM1 4.61 6.58 5.23 6.10 6.29 4.39 5.07 2.50 6.83 2.56 3.88 7.23
## 4 SM2 5.16 5.38 5.34 7.20 6.01 4.41 5.49 2.81 5.78 2.07 4.55 6.33
## 5 SM3 4.70 6.07 5.87 6.01 4.51 5.85 5.77 2.48 6.10 3.18 3.02 6.53
## 6 SM4 5.30 6.32 4.80 8.09 6.72 6.49 5.36 3.57 6.02 2.93 2.42 7.47
## # … with 4 more variables: SKo92 <dbl>, SKg92 <dbl>, SKk92 <dbl>, ID92 <dbl>なお、この形質では全16環境でデータが収集されているが、他の形質では抜けがある。そこで、データが無い列を除いておく。この形質では、何も変化しない。
selector <- apply(!is.na(td.w), 2, sum) > 0 # 非欠測のデータが一つでもあればTRUEになる
td.ws <- td.w[, selector]抽出した形質データをマーカー遺伝子型データと融合する。
元のデータには形質データが含まれていない（系統名のみがデータとして含まれている）。
head(geno$pheno)## gen
## 1 SM1
## 2 SM2
## 3 SM3
## 4 SM4
## 5 SM5
## 6 SM6そこで、この系統名をもとに、2つのデータを融合する。
cross <- geno
cross$pheno <- geno$pheno %>% left_join(td.ws, by = "gen")
cross## This is an object of class "cross".
## It is too complex to print, so we provide just this summary.
## Doubled haploids
##
## No. individuals: 150
##
## No. phenotypes: 17
## Percent phenotyped: 100 99.3 99.3 99.3 99.3 99.3 99.3 99.3 99.3 99.3 99.3
## 99.3 99.3 99.3 99.3 99.3 99.3
##
## No. chromosomes: 7
## Autosomes: 1 2 3 4 5 6 7
##
## Total markers: 223
## No. markers: 37 37 31 33 29 22 34
## Percent genotyped: 96
## Genotypes (%): AA:50.9 BB:49.1表現型の数（No. phenotypes）が17になった。実際は、16で、1つは系統名である。
では、種々のplot関数を使って （連鎖地図、欠測値、表現型データのヒストグラム）を表示してみる。
連鎖地図を表示する。
plotMap(cross)
欠測しているマーカー遺伝子型データを表示する。
plotMissing(cross)
黒い点が欠測している場所である。
表現型データのヒストグラムを描く。
plotPheno(cross)
なんだか奇妙なヒストグラムが描かれるが、 これは、表現型データの最初にある系統名のヒストグラムである。図示する意味はない。 2列目のデータ（環境１の収量）のデータを表示する。
plotPheno(cross, pheno.col = 2)
個別に図をながめるのもよいが、まとめてPDFに出力しておくと、あとで見返せて便利である。
pdf("out/cross_summary.pdf") # set the name of pdf file
plotMap(cross)
plotMissing(cross)
for(i in 2:nphe(cross)) {
plotPheno(cross, pheno.col = i)
}
dev.off() # should close the file at the end## png
## 2PDF関数を用いたときには、必ず最後にファイルを閉じることを忘れないようにする。 QTL解析のために、染色体上に等間隔に取られた位置 （pseudo markerなどとよばれる）における遺伝子型の確率を計算しておく。 ここでは2 cM間隔でpseudo markerを配置する。
interval <- 2 # 2 cMに設定。この値を変えると計算する密度を変えられる
cross <- calc.genoprob(cross, step = interval) # 遺伝子型の確率を計算まずは、単純インターバルマッピング（simple interval mapping、SIM）を行う。 ここでは、16環境で計測された収量データのうち、16番目（env.id = 16）を解析する こととする。
env.id <- 16 # この値を変えると他の環境についてQTL解析ができるRのqtlパッケージでは、いくつかの方法でインターバルマッピングを行う ことができる。ここでは、EMアルゴリズムに基づく方法（“em”）、Haley and Knottの 回帰に基づく方法（“hk”）の2つを試す。 qtlパッケージに含まれるscanoneという関数を用いてSIMを行う。 まずは、EMアルゴリズムに基づく方法でQTL解析。
out.em <- scanone(cross, pheno.col = env.id + 1, method = "em")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.plot(out.em)
大きなQTLが、第2、第3染色体に検出される。
では、次に、Haley and
Knott回帰でQTL解析を行ってみる。
out.hk <- scanone(cross, pheno.col = env.id + 1, method = "hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.plot(out.hk)
EMアルゴリズムとほぼ同じ結果である。 2つの手法による結果を同時に描いて、比較をしてみよう。
par(mfrow = c(3,1))
plot(out.em, col = "blue", lty = "solid", main = "EM-algorithm")
plot(out.hk, col = "red", lty = "dotted", main = "Haley and Knott regression")
plot(out.em, out.hk, col = c("blue", "red"),
lty = c("solid", "dotted"), main = "Comparison")
par(mfrow = c(1,1)) 上述したように、2つの方法の結果はほとんど同じである。 これは、このデータではマーカーが高密度に配置されているために、 マーカー間に存在するQTLの計算方法が違う（マーカー上での計算は同じ） 2手法の違いが出なかったためだと考えられる。そこで、これ以降の計算は、 2手法のうち標準的に用いられる“em”を用いた解析を行う。 EMを使った解析で検出されたQTLをリストアップする。
summary(out.em)## chr pos lod
## dRcs1 1 32.4 0.347
## MWG858 2 39.3 9.625
## c3.loc58 3 58.0 7.217
## c4.loc80 4 80.0 1.938
## ABC261 5 139.1 2.755
## cMWG652A 6 26.6 1.192
## ABC482 7 171.6 2.092次にLOD値の大きさに基づきQTLを検出するためのしきい値を決定する。 まず、しきい値を決める前にQTLのリストを表示する。するとLOD値のピークから 判断されたQTLのリストが表示される。次に、しきい値を決定するために並べ替え 検定（permutation test）を行う。並べ替え検定では、表現型データを無作為に 並び替えながらQTL解析を繰り返し行う。無作為に並び替えることにより、遺伝子 型データと表現型データの対応関係が撹乱され、データが本来もっていたQTLに関 する情報が消失してしまう。逆に見れば、並び替えを行ったデータでなお検出され るようなLOD値のピークは偽物のQTL（偽陽性）に起因するものと考えられる。 このようなアイデアに基づき、並び替える毎にQTL解析を行い、染色体全域の中でLOD が最大になったときの値を記録する。これを複数回（例えば1,000回）繰り返すこと により、QTLが無い場合のLOD値の分布（経験的帰無分布）を求める。 しきい値は、例えば、こうして得られる経験的帰無分布の上位5%点等に設定する。 例えば、上位5％点に設定すると、QTLが無いのに誤って偽のQTLを検出してしまう 確率が「ゲノム全体で5%」になるように設定される。
operm.em <- scanone(cross, pheno.col = env.id + 1, method = "em", n.perm = 1000, verbose = F)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.こうして求めた閾値を表示する。
summary(operm.em, alpha = 0.05) ## LOD thresholds (1000 permutations)
## lod
## 5% 2.6得られた閾値をもとにLOD値のピークのうち有意なものを求めると、このデータの場合、2つのQTL領域が検出される。
plot(out.em)
abline(h = summary(operm.em, alpha = 0.05))
次に、コンポジットインターバルマッピング（composite interval mapping、CIM）を 行ってみよう。CIMでは、マーカー共変量の数（ここでは5に設定）と解析の窓の大きさ（その窓の中には共変量マーカーを配置しない。ここでは10（cM）に設定）を指定して解析を行う必要がある。
n.covar <- 5
window.size <- 10
outcim.em <- cim(cross, pheno.col = env.id + 1, method = "em", n.marcovar = n.covar, window = window.size)SIMとCIMの結果を図示して比較すると、前者に比べ後者はQTLの位置の解像度が高いことが分かる。これはCIMでは共変量を用いることで、注目している場所以外に位置しているQTLの影響をモデルに取り込むことができ、それらが引き起こす変動を押さえ込むことで、注目している場所にあるQTLのシグナルをより明瞭に捉えられるようになる。
par(mfrow = c(3,1))
plot(out.em, col = "blue", main = "Simple Interval Mapping")
plot(outcim.em, col = "red", main = "Composite Interval Mapping")
add.cim.covar(outcim.em, col = "green")
plot(out.em, outcim.em, col = c("blue", "red"), main = "Comparison")
add.cim.covar(outcim.em, col = "green")
par(mfrow = c(1,1))ここでも、SIMと同様にして並べ替え検定を行って 閾値を決定する。ここでは時間の節約のために200回の繰り返しで計算している。 実際にはもう少し多く、例えば、1000回くらいは繰り返したほうがよい。
opermcim.em <- cim(cross, pheno.col = env.id + 1, method = "em", n.marcovar = n.covar, window = window.size, n.perm = 200)5%の有意水準で閾値を決める。
summary(opermcim.em, alpha = 0.05)## LOD thresholds (200 permutations)
## [,1]
## 5% 4.4その閾値のもとで有意なマーカーをリストアップする。
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## c2.loc40 2 40 14.2
## c3.loc58 3 58 10.4結果は、（繰り返しが少ないため）毎回変わるかもしれない。
最後にQTLが検出された第2、第3染色体だけをクローズアップしてみる。
chr.id <- c(2,3)
plot(outcim.em, chr = chr.id)
abline(h = summary(opermcim.em, alpha = 0.05))
次に検出されたQTLの効果の推定を行ってみよう。 makeqtl関数で検出されたQTLの位置の遺伝子型を抜き出し、 それを独立変数（説明変数）として、表現型データ（従属変数、被説明変数）を回帰する。
まずは、上で検出された2つのQTLについて情報をtempという入れ物に保存する。
temp <- summary(outcim.em, perms = opermcim.em, alpha = 0.05)
qtl <- makeqtl(cross, chr = temp$chr, pos = temp$pos, what = "prob")
res <- fitqtl(cross, qtl = qtl, get.ests = T, method = "hk", pheno.col = env.id + 1)## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 1 individuals with missing phenotypes.QTL解析の結果を表示する
plot(qtl) 
回帰分析における\(F\)検定の結果を表示する。
summary(res) ##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 149
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 2 47.02602 23.513011 17.44145 41.67073 0 0
## Error 146 65.82542 0.450859
## Total 148 112.85144
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
## 2@40.0 1 25.07 10.441 22.22 55.61 0 7.25e-12 ***
## 3@58.0 1 18.85 8.149 16.71 41.82 0 1.41e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 4.95809 0.05539 89.510
## 2@40.0 -0.41588 0.05577 -7.457
## 3@58.0 -0.36318 0.05616 -6.4671つ目のQTLの説明する分散は全分散の22.45%、2つ目のQTLは16.17%である。 また、2つのQTLを合わせた変動は、41.9%である。 さらに、1つ目のQTLについて、遺伝子型AAのBBに対する効果は、-0.414、 2つ目のQTLについては、-0.358である。
付録：ファイルの読み込みかた
ここでは、マーカー遺伝子型データと表現型データをコンマ区切りテキスト（csv）ファイルから読み込む方法を説明する。なお、ここで用いるデータは、agridatと同一のソースからのデータであるが、agridatと一部（マーカー数など）が異なる。
dataディレクトリに含まれている、geno_map.csvとyieldを読み込んで、解析データセットcrossを作成する。
なお、このデータは、倍化半数体のためcrosstype = “dh”（doubled haploids）としてデータを読み込む。
cross <- read.cross(format = "csvs", genfile = "data/geno_map.csv", phefile = "data/pheno/yield.csv", crosstype = "dh")## --Read the following data:
## 150 individuals
## 495 markers
## 17 phenotypes## Warning in summary.cross(cross): Some markers at the same position on chr 2; use
## jittermap().## Warning in summary.cross(cross): Invalid genotypes.
## Observed genotypes: 1 3## Warning in summary.cross(cross): Strange genotype pattern on chr 7.## --Cross type: dhcross## This is an object of class "cross".
## It is too complex to print, so we provide just this summary.## Warning in summary.cross(x): Some markers at the same position on chr 2; use
## jittermap().## Warning in summary.cross(x): Invalid genotypes.
## Observed genotypes: 1 3## Warning in summary.cross(x): Strange genotype pattern on chr 7.## Doubled haploids
##
## No. individuals: 150
##
## No. phenotypes: 17
## Percent phenotyped: 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 100 100 100
##
## No. chromosomes: 7
## Autosomes: 1 2 3 4 5 6 7
##
## Total markers: 495
## No. markers: 60 78 81 60 93 56 67
## Percent genotyped: 95.8
## Genotypes (%): AA:100.0 BB:0.0## Warning in summary.cross(x): Some markers at the same position on chr 2; use
## jittermap().## Warning in summary.cross(x): Invalid genotypes.
## Observed genotypes: 1 3## Warning in summary.cross(x): Strange genotype pattern on chr 7.おかしな遺伝子型のパターンを示すマーカーがあるようだ。Genotypeのパターンを見ると、BBとして読み込まれているデータがない。実は、“B”として入力しているものが正しく”BB”として認識されていない（バグだと思われる）。
“B”は”H”として入力すると、正しくBBとして認識されるようである。
そこで、“B”を”H”に置き換えたファイルを作成して、以下のようにして読み込む。
cross <- read.cross(format = "csvs", genfile = "data/geno_bc.csv", phefile = "data/pheno/yield.csv", crosstype = "dh")## --Read the following data:
## 150 individuals
## 495 markers
## 17 phenotypes## Warning in summary.cross(cross): Some markers at the same position on chr 2; use
## jittermap().## --Cross type: dhcross## This is an object of class "cross".
## It is too complex to print, so we provide just this summary.## Warning in summary.cross(x): Some markers at the same position on chr 2; use
## jittermap().## Doubled haploids
##
## No. individuals: 150
##
## No. phenotypes: 17
## Percent phenotyped: 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 100 100 100
##
## No. chromosomes: 7
## Autosomes: 1 2 3 4 5 6 7
##
## Total markers: 495
## No. markers: 60 78 81 60 93 56 67
## Percent genotyped: 95.8
## Genotypes (%): AA:51.2 BB:48.8## Warning in summary.cross(x): Some markers at the same position on chr 2; use
## jittermap().まだ、警告がでるが、これは、第2染色体の同じ位置に2つのマーカーが座乗していることによる。
このような場合は、jittermap関数によりこれらマーカー間に小さな隙間をつくる。
cross <- jittermap(cross)
cross## This is an object of class "cross".
## It is too complex to print, so we provide just this summary.
## Doubled haploids
##
## No. individuals: 150
##
## No. phenotypes: 17
## Percent phenotyped: 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 100 100 100
##
## No. chromosomes: 7
## Autosomes: 1 2 3 4 5 6 7
##
## Total markers: 495
## No. markers: 60 78 81 60 93 56 67
## Percent genotyped: 95.8
## Genotypes (%): AA:51.2 BB:48.8警告はなくなる。
読み込まれたデータについて、連鎖地図を表示する。
plotMap(cross)