本実験では、Rのパッケージagridatに含まれているオオムギの多環境試験データを例として解析を行う。 なお、表現型データについては以下からダウンロードすることもできる。 http://wheat.pw.usda.gov/ggpages/SxM/phenotypes.html
このデータセットには、オオムギ品種Steptoe×Morexの交配から得られた150の倍加半数体（doubled haploid、DH）系統のデータが含まれている。 遺伝子型データには7つの染色体全体に分布する218のマーカーの遺伝子型データが含まれる。 形質データは、最大16環境で計測された表現型データが含まれている。
まずは、必要なパッケージを読み込む。
require(qtl) # R package for QTL mapping
require(agridat) # QTL data
require(tidyverse) # for rearranging data structure 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”以降の列が形質データである。
ここでは、草丈”height”を解析対象とする。
target <- "height"
td <- pheno %>% select(gen, env, all_of(target))
head(td)## gen env height
## 1 Steptoe MN92 84.5
## 2 Steptoe MTi92 NA
## 3 Steptoe MTd92 75.5
## 4 Steptoe ID91 111.0
## 5 Steptoe OR91 90.0
## 6 Steptoe WA91 112.0なお、このデータは、復数の環境のデータが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… 84.5 NA 75.5 111 90 112 98 82 77.5 95 NA 78
## 2 Morex 97 NA 90.5 116 110 133 100 97 95 110 NA 94.5
## 3 SM1 114. 120 98.5 117. 100 146 120 93 85 100 80 106.
## 4 SM2 94.5 105 76.5 113 102. 123 105 84 82.5 85 62.5 87
## 5 SM3 88 105 82.5 117. 80 128 105 86 85 95 59.5 87
## 6 SM4 108. 113 86.5 119. 102. 126. 113 84 95 105 62.5 97
## # … with 4 more variables: SKo92 <dbl>, SKg92 <dbl>, SKk92 <dbl>, ID92 <dbl>
## # ℹ Use `colnames()` to see all variable namesなお、この形質では全16環境でデータが収集されているが、形質によっては環境の抜けがある。そこで、データが無い列（計測が行われなかった環境）を除いておく。この形質では、何も変化しない。
selector <- apply(!is.na(td.w), 2, sum) > 0 # 非欠測のデータが一つでもあればTRUEになる
td.ws <- td.w[, selector]
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… 84.5 NA 75.5 111 90 112 98 82 77.5 95 NA 78
## 2 Morex 97 NA 90.5 116 110 133 100 97 95 110 NA 94.5
## 3 SM1 114. 120 98.5 117. 100 146 120 93 85 100 80 106.
## 4 SM2 94.5 105 76.5 113 102. 123 105 84 82.5 85 62.5 87
## 5 SM3 88 105 82.5 117. 80 128 105 86 85 95 59.5 87
## 6 SM4 108. 113 86.5 119. 102. 126. 113 84 95 105 62.5 97
## # … with 4 more variables: SKo92 <dbl>, SKg92 <dbl>, SKk92 <dbl>, ID92 <dbl>
## # ℹ Use `colnames()` to see all variable names抽出した形質データをマーカー遺伝子型データと融合する。
元のデータには形質データが含まれていない（系統名のみがデータとして含まれている）。
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 98.7 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列目のデータ（1番目の環境のデータ）のデータを表示する。
plotPheno(cross, pheno.col = 2)
個別に図をながめるのもよいが、まとめてPDFに出力しておくと、あとで見返せて便利である。
pdf("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## quartz_off_screen
## 2PDF関数を用いたときには、必ず最後にファイルを閉じることを忘れないようにする。 QTL解析のために、染色体上に等間隔に取られた位置 （pseudo markerなどとよばれる）における遺伝子型の確率を計算しておく。 ここでは1 cM間隔でpseudo markerを配置する。
interval <- 1 # 1 cMに設定。この値を変えると計算する密度を変えられる
cross <- calc.genoprob(cross, step = interval) # 遺伝子型の確率を計算まずは、単純インターバルマッピング（simple interval mapping、SIM）を行う。 ここでは、16環境で計測されたデータのうち、6番目（env.id = 6）を解析する こととする。
env.id <- 6 # この値を変えると他の環境について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
## c1.loc98 1 98.0 4.075
## MWG858 2 39.3 9.007
## c3.loc55 3 55.0 14.684
## ABG397 4 129.7 0.885
## His3B 5 100.5 0.768
## MWG798A 6 156.6 0.339
## c7.loc190 7 190.0 0.896次に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.表現型データが無い個体が1個体あるために除いた、と警告がでるが、解析は問題なく行われる。
こうして求めた閾値を表示する。
summary(operm.em, alpha = 0.05) ## LOD thresholds (1000 permutations)
## lod
## 5% 2.69得られた閾値をもとにLOD値のピークのうち有意なものを求めると、このデータの場合、およそ4つ（第1染色体は厳密には複数に分かれるかもしれない）の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と同様にして並べ替え検定を行って 閾値を決定する。ここでは1000回の並べ替えを行う。少し時間がかかる。
opermcim.em <- cim(cross, pheno.col = env.id + 1, method = "em", n.marcovar = n.covar, window = window.size, n.perm = 1000)5%の有意水準で閾値を決める。
summary(opermcim.em, alpha = 0.05)## LOD thresholds (1000 permutations)
## [,1]
## 5% 4.27その閾値のもとで有意なマーカーをリストアップする。
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## c1.loc107 1 107.0 9.61
## MWG858 2 39.3 18.97
## c3.loc55 3 55.0 26.99結果は、無作為並び替えであるため、若干変わるかもしれない。
最後にQTLが検出された第1, 2、第3染色体をクローズアップしてみる。
chr.id <- c(1, 2,3)
plot(outcim.em, chr = chr.id)
abline(h = summary(opermcim.em, alpha = 0.05))
次に検出されたQTLの効果の推定を行ってみよう。 makeqtl関数で検出されたQTLの位置の遺伝子型を抜き出し、 それを独立変数（説明変数）として、表現型データ（従属変数、被説明変数）を回帰する。
まずは、上で検出された3つの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 + Q3
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 3 12430.235 4143.41153 36.36847 67.50379 0 0
## Error 145 5983.893 41.26823
## Total 148 18414.128
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
## 1@107.0 1 2120 9.812 11.51 51.37 0 3.59e-11 ***
## 2@39.3 1 3788 15.867 20.57 91.78 0 < 2e-16 ***
## 3@55.0 1 5874 22.129 31.90 142.34 0 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 126.0304 0.5303 237.663
## 1@107.0 -4.1257 0.5756 -7.167
## 2@39.3 5.0599 0.5282 9.580
## 3@55.0 6.4390 0.5397 11.9311つ目のQTLの説明する分散は全分散の11.51%、2つ目のQTLは20.57%、3つ目のQTLは31.90%である。 また、3つのQTLを合わせた変動は、67.50%である。 さらに、1つ目のQTLについて、遺伝子型AAのBBに対する効果は、-4.13、2つ目のQTLは5.06、3つ目のQTLは6.44である。
付録：ファイルの読み込みかた
ここでは、マーカー遺伝子型データと表現型データをコンマ区切りテキスト（csv）ファイルから読み込む方法を説明する。なお、ここで用いるデータは、agridatと同一のソースからのデータであるが、agridatと一部（マーカー数など）が異なる。
dataディレクトリに含まれている、geno_map.csvとyieldを読み込んで、解析データセットcrossを作成する。
なお、このデータは、倍化半数体のためcrosstype = “dh”（doubled haploids）としてデータを読み込む。
cross <- read.cross(format = "csvs", genfile = "geno_map.csv", phefile = "height.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 99.3 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 = "geno_bc.csv", phefile = "height.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 99.3 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 99.3 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)