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

Size: px

Start display at page:

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

ありかつまつかた
5 years ago
Views:

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

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) # 欠損データをプロットする

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.

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

p 値が極端に小さいマーカーは分離比が理論値よりも大きく歪んでおり誤判定している可能性がありますそれらのマーカーは再度ジェノタイピングするまたは代替マーカーを使用するなど対処する必要があるでしょうただし例データのように同一染色体上の隣接するマーカーごと歪む場合 ( 染色体 4 の RM307, RM261, RM185) や分離比が歪む原因が遺伝的なものである場合はそのま

4 p 値が極端に小さいマーカーは分離比が理論値よりも大きく歪んでおり誤判定している可能性がありますそれらのマーカーは再度ジェノタイピングするまたは代替マーカーを使用するなど対処する必要があるでしょうただし例データのように同一染色体上の隣接するマーカーごと歪む場合 ( 染色体 4 の RM307, RM261, RM185) や分離比が歪む原因が遺伝的なものである場合はそのまま使用すべきでしょう例えばボンフェローニ補正前の p 値が以下の歪みの激しいマーカーを削除する場合には以下のようにします 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 の結果を表示する plot(out.simhk, out.cimhk, col=c("blue", "red"),show.marker.names=true)

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 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 ピークを表示

av1 は Add モデル -One モデルの LOD スコア上三角に相互作用モデル (Full-Add) 下三角は Full モデルの LOD をプロットするには以下のコマンドを使います

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 モデル ) に plot(out2.hk, lower="fv1", upper="av1", zscale=true)

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)

しますつまり操作をなぞる場合は赤字部分をコピー&ペーストします入力コマンドの各行の# 以降右はコメントでなくても操作に影響しません library(qtl) # qtl ライブラリーを起動する testmap <- read.cross(

しますつまり操作をなぞる場合は赤字部分をコピー&ペーストします入力コマンドの各行の# 以降右はコメントでなくても操作に影響しません library(qtl) # qtl ライブラリーを起動する testmap <- read.cross( R/qtl による連鎖地図の作成清水顕史 R の利用 R の qtl ライブラリーは連鎖地図も作成できます 1 そのための遺伝子型および表現型データは次のフォーマットで準備しますデータは表現型 (QTL 解析用 ) 遺伝子型の順に上から並べます 1 列目は表現型名またはマーカー名