必要なパッケージの読み込み

require(qtl)
## Loading required package: qtl
require(agridat)
## Loading required package: agridat
require(tidyverse)
## 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()

最小コードでQTL解析を実行

このプログラムでは、最小限のコードで上述した内容を実行する。

agridatからデータを読み込む。対象形質は”amylase”とするが、そこを別の形質に変えるとその形質を解析できる。

target <- "amylase"

tmp <- steptoe.morex.pheno %>% select(gen, env, all_of(target)) %>% 
  pivot_wider(id_cols = gen, names_from = env, values_from = all_of(target))
selector <- apply(!is.na(tmp), 2, sum) > 0 # 非欠測のデータが一つでもあればTRUEになる
tmp <- tmp[, selector]

cross <- steptoe.morex.geno
cross$pheno <- steptoe.morex.geno$pheno %>% left_join(tmp, 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:     10 
##     Percent phenotyped: 100 99.3 99.3 98 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

2cM間隔で確率を計算する

interval <- 2  # 2 cM、これを変えると間隔を変えられる
cross <- calc.genoprob(cross, step = interval) # 確率計算

計測された環境を調べる

colnames(cross$pheno) # 表現型データの列名をチェック
##  [1] "gen"   "MN92"  "MTi92" "MTd92" "ID91"  "OR91"  "WA91"  "MTi91" "WA92" 
## [10] "ID92"

全部で9環境で計測されている。

CIMを行う

env.id <- 2  # MTi92は2番目。ここを変えるといろいろな環境に対して解析できる。
n.covar <- 5 # 共変量マーカーの数
window.size <- 10 # 10 cM以内には、共変量マーカーをおかない
outcim.em <- cim(cross, pheno.col = env.id + 1, method = "em", 
                 n.marcovar = n.covar, window = window.size) 
# QTL解析の実行

並び替え検定。ここでは、少し時間はかかるが200回並び替える。

opermcim.em <- cim(cross, pheno.col = env.id + 1, method = "em", 
                   n.marcovar = n.covar, window = window.size, 
                   n.perm = 200)
# n.permで繰り返し数を指定している。

5%の有意水準でしきい値を決める。

summary(opermcim.em, alpha = 0.05) # しきい値を示す。
## LOD thresholds (200 permutations)
##    [,1]
## 5% 4.08

そのしきい値のもとで有意なマーカーをリストアップする。

summary(outcim.em, perms = opermcim.em, alpha = 0.05) # QTlのリストアップ
##          chr  pos  lod
## c2.loc76   2 76.0 4.09
## ABC302     7 78.3 9.41

プロットをする

plot(outcim.em) # QTL解析結果のLOD得点のプロット
abline(h = summary(opermcim.em, alpha = 0.05)) # しきい値のところに線をひく