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

(1)

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 で示した区画に分かれて配置されています。

(2)

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

(3)

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 _F₂_{分離集団では、共優性マーカーの遺伝子型の分離が}_{A：H：B=1:2:1 に、優性マーカーは A:C または} B:D=1:3 に分離することが期待でされる。 4 _{マーカーが単一遺伝子座由来でない場合や、ホモ・ヘテロの判別が難しく誤った遺伝子型を含んでいる} 場合などが考えられる。 5 _{ボンフェローニ補正は、χ}2_{適合度検定による}_p値に、調査した全マーカー数を掛ける。

(4)

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 が存在する確率が高いと考えます。

(5)

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 の結果を表示する

(6)

有意になったマーカー間の相互作用効果をプロットしてみます。インターバルマッピン

グで有意になった2 地点(染色体 3 の 147cM と染色体 4 の 22cM)の遺伝子型別に表現型値

をプロットし、線分が平行でないとき交互作用があるとみなします。 effectplot(testmap, mname1="3@147", mname2="4@22") #交互作用

Two-QTL スキャン

scantwo関数を用いると、一度に2領域のスキャンを行い、インターバル間の相互作用を探索することができます。以下のモデルに対してそれぞれ対数尤度(LOD)を計算し、その差を LODスコアとして検定に用います。 Fullモデル：

y









₁

q

₁





₂

q

₂





₃



q

₁



q

₂







Addモデル：

y









₁

q

₁





₂

q

₂





Full-Addモデル： y



₃



q₁q₂







Oneモデル :

y









₁

q

₁





Nullモデル：

y









testmap <- calc.genoprob(testmap, step=1) #1cM 毎に QTL を検索する out2.hk <- scantwo(testmap, method="hk") # two QTL のスキャン

(7)

(注）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 モデル) に

(8)

#並べ替え検定で有意水準を推定する

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)