R で QTL 解析 以下で R への入力コマンドはゴシック赤字で表記しています # より右はコメントなの で入力の必要はありません 操作を再現する際 タイプミスに注意しましょう データの読み込み qtl ライブラリーを起動し ファイル IN-RIL.csv を読み込みます library(qtl)

R で QTL 解析 以下で R への入力コマンドはゴシック赤字で表記しています # より右はコメントなの で入力の必要はありません 操作を再現する際 タイプミスに注意しましょう データの読み込み qtl ライブラリーを起動し ファイル IN-RIL.csv を読み込みます library(qtl) # ライブラリー起動 testmap <- read.cross("csvr", file="in-ril.csv", estimate.map=false)# データ読み込み testmap <- convert2riself(testmap) # RI self としてデータを処理する summary(testmap) # 交配タイプ 個体数 マーカー数 表現型数 を確認する 読み込んだデータを plot.rf 関数で確認します plot.rf(testmap, alternate.chrid=true) # マーカー間連鎖の確認 出力される図は 連鎖地図上に並べたマーカー間の連鎖の程度に応じて色分けされていて 赤色が連鎖 青色が独立であることを意味します 連鎖地図上の位置に問題がない場合 隣接するマーカーが連鎖するので図の対角線上のみが赤くなります イネの 12 種類の染色体は 上端および右端の 1~12 で示した区画に分かれて配置されています

読み込んだデータを plot 関数で確認します plot(testmap) # 欠損データ 連鎖地図 形質別ヒストグラムを表示 一般的な QTL 解析は 形質が正規分布に従う統計モデルを仮定しています 分布が大き くズレる場合は 外れ値の有無や 形質値データの精度 再現性に問題がないか確認します 欠損データの確認遺伝子型が判別できなかった欠損データを確認しておきます plot.missing(testmap) # 欠損データをプロットする

par(mfrow=c(1,2), las=1) plot(ntyped(testmap), ylab="no. typed markers", main="no. genotypes by individual") plot(ntyped(testmap, "mar"), ylab="no. typed individuals", main="no. genotypes by marker") 欠損データの多い 1 マーカーや個体は とりあえず解析から外します 必要に応じて ジェノタイピングし直すか近傍マーカーを作り直します データの削除は以下の通りです 2 例えばデータ数が 50 以上の個体のみを採用する場合 testmap <- subset(testmap, ind=(ntyped(testmap)>50)) 例えば個体数が 70 未満のマーカーを解析から除外する場合 nt.bymar <- ntyped(testmap, "mar") todrop <- names(nt.bymar[nt.bymar < 70]) testmap <- drop.markers(testmap, todrop) 遺伝子型情報のチェック分離比の歪んだマーカーを検出する 各マーカー遺伝子座の分離比が期待値から大きくずれる場合 3 遺伝子型の判別方法に問題 4 がある可能性が高い 期待比からのズレは geno.table 関数を使って調べることができ ボン フェローニ補正 5 した 5% 水準以下のマーカーは以下のように検出できる gt <- geno.table(testmap) gt[gt$p.value < 0.05/totmar(testmap),] 1 明確な基準はありませんが データの 20% 以上が欠損する場合は怪しいと考えましょう 2 入力ファイルを編集してもよいでしょう 3 F2 分離集団では 共優性マーカーの遺伝子型の分離が A:H:B=1:2:1 に 優性マーカーは A:C または B:D=1:3 に分離することが期待でされる 4 マーカーが単一遺伝子座由来でない場合や ホモ ヘテロの判別が難しく誤った遺伝子型を含んでいる場合などが考えられる 5 ボンフェローニ補正は χ 2 適合度検定による p 値に 調査した全マーカー数を掛ける

p 値が極端に小さいマーカーは 分離比が理論値よりも大きく歪 んでおり 誤判定している可能性があります それらのマーカーは 再度ジェノタイピングする または代替マーカーを使用するなど対 処する必要があるでしょう ただし 例データのように 同一染色 体上の隣接するマーカーごと歪む場合 ( 染色体 4 の RM307, RM261, RM185) や 分離比が歪む原因が遺伝的なものである場合はそのま ま使用すべきでしょう 例えばボンフェローニ補正前の p 値が 1 10-10 以下の歪みの激し いマーカーを削除する場合には以下のようにします todrop <- rownames(gt[gt$p.value < 1e-10,]) testmap <- drop.markers(testmap, todrop) 連鎖地図の検定 プロット plot.map() 関数で 引数に show.marker.names=true を用いると マーカー名付きの 連鎖地図でギャップ ( 連鎖地図のマーカー密度が低い箇所 ) を表示できます plot.map(testmap, show.marker.names=true) single-qtl スキャン 連鎖地図の隣接マーカー間 ( マーカー インターバル ) に着目し QTL を検定します LOD スコアが高いほど そのインターバル間に QTL が存在する確率が高いと考えます

testmap <- calc.genoprob(testmap, step=1) #1cM 毎に QTL 検索 #interval mapping(im) を行う (Haley-Knott 回帰に基づく method="hk") #1 番目の形質を解析する場合に pheno.col=1 とする 2 番目の形質なら pheno.col=2 out.simhk <- scanone(testmap, pheno.col=1, method="hk") QTL の有意水準 (pval) は 対象形質の分布の形 ( 遺伝様式を反映 ) や連鎖地図の大きさによって変動し 各形質の並替え検定によって推定するのが一般的です QTL 解析のような多重比較検定では マーカー毎に有意確率を決めると連鎖地図全体の有意確率が不当に高くなってしまう ( 本当は有意でない QTL を検出してしまう ) ため 地図全体の有意確率を適切に制御する必要があります 並べ替え検定では形質値の無作為な並べ替えを 1000 回以上行い それぞれ最も有意確率が低くなる値を集めた分布の α% 点の LOD スコアを水準 α の閾値とみなします # 計算時間の都合上 実習の中では 100 回のシミュレート (n.perm=100) とする # もしくは マルチタスク機能が使える PC の場合は (n.perm=1000,n.cluster=2~) とする operms <- scanone (testmap, n.perm=100, method="hk") #operms <- scanone (testmap, n.perm=1000, method="hk",n.cluster=2) summary(out.simhk, perms=operms, alpha=0.05, pvalues=true) #composite interval mapping(cim) を行う out.cimhk <- cim(testmap, n.marcovar=5, pheno.col=1, method="hk") # 並べ替え検定を行う,IM と同じ取扱い operm <- cim(testmap,n.perm=100,n.marcovar=5,method="hk") #operm <- cim(testmap,n.perm=1000,n.marcovar=5,method="hk",n.cluster=2) #5% 水準で有意な QTL ピークを書き出す summary(out.cimhk, perms=operm, alpha=0.05, pvalues=true) #IM と CIM の結果を表示する plot(out.simhk, out.cimhk, col=c("blue", "red"),show.marker.names=true)

有意になったマーカー間の相互作用効果をプロットしてみます インターバルマッピングで有意になった 2 地点 ( 染色体 3 の 147cM と染色体 4 の 22cM) の遺伝子型別に表現型値をプロットし 線分が平行でないとき交互作用があるとみなします effectplot(testmap, mname1="3@147", mname2="4@22") # 交互作用 Two-QTL スキャン scantwo 関数を用いると 一度に2 領域のスキャンを行い インターバル間の相互作用を探索することができます 以下のモデルに対してそれぞれ対数尤度 (LOD) を計算し その差を LODスコアとして検定に用います Fullモデル : y q q q 1 1 2 2 3 1 q2 1q1 2q2 3 q1 q2 q 1 1 Addモデル : y Full-Addモデル : y Oneモデル : y Nullモデル : y testmap <- calc.genoprob(testmap, step=1) #1cM 毎に QTL を検索する out2.hk <- scantwo(testmap, method="hk") # two QTL のスキャン summary(out2.hk) # 染色体対毎の最大 LOD ピークを表示

( 注 )lod.full は Full モデル lod.fv1 は Full モデル -One モデル lod.int は Full モデル - Add モデル lod.add は Add モデルの lod.av1 は Add モデル -One モデルの LOD スコア 上三角に相互作用モデル (Full-Add) 下三角は Full モデルの LOD をプロットするには 以下のコマンドを使います カラーマップが赤いほど LOD スコアが高く 連鎖地図上の位 置を表す横軸の点と縦軸の点 plot(out2.hk, zscale=true) # 下三角を Full-One( 条件付き相互作用モデル ) 上三角を Add-One( 条件付き add モデル ) に plot(out2.hk, lower="fv1", upper="av1", zscale=true)

# 並べ替え検定で有意水準を推定する operm2 <- scantwo(testmap, method="hk", n.perm=10) #operm2 <- scantwo(testmap, method="hk", n.perm=1000,n.cluster=2) summary(out2.hk, perms=operm2, alpha=0.05, pvalues=true)