1
R/qtl による連鎖地図の作成
清水顕史R の利用
R の qtl ライブラリーは連鎖地図も作成できます1。そのための遺伝子型および表現型デ ータは次のフォーマットで準備します。データは表現型（QTL 解析用）、遺伝子型の順に上 から並べます。 1 列目は表現型名またはマーカー名を記入します。表現型データは最低 1 つ記入しておく必 要があるのでデータがない場合も何かデタラメな値を入力します。2 列目はマーカーの座乗 する染色体番号または連鎖群番号（任意の数字）を記入します（不明の場合も何か記入し ておく）。3 列目は必要な場合に cM 単位で連鎖地図位置を記入します（連鎖地図情報が記 入されていない場合、マーカーは連鎖地図順に並べられていると見なしてデータが読み込 まれます。マーカーは物理地図位置の順に入力します）。表現型データの2、3 列目は空白 にしておきます。遺伝子型データは、母親ホモ接合体をA、父親ホモ接合体を B、ヘテロ接 合体をH とし、優性マーカーの場合は、A と C(B or H)か、B と D(A or H)と表記します。 （注意1）表現型データや遺伝子型データは、測定の際にコンタミやエラーを起こさないよ うに注意を払うのは当然ですが、記録したものをパソコンに入力する際にもタイプミスが 生じます。入力したデータは印刷し、必ず見直しましょう（複数人で確認し合うとよい）。 （注意2）解析（マッピング）は一つの個体に由来する分離集団ごとに計算します。たとえ 同じ交配組合せのF2集団でも、由来するF1個体が異なる場合は別の集団として解析します （減数分裂は環境要因の影響を受けることがあり、分離比が歪むかもしれない）。 2-1)データ入力 例題用データ(test.csv)を、例えば C ドライブの Rdata フォルダに置きます。R 起動後に[フ ァイル]メニューの[ディレクトリの変更…]からデータのあるフォルダ（ここでは C:¥Rdata） を指定します。以下では、R への入力コマンドは赤字で記し、出力結果は青字で示すことに 1 研究室のパソコンに導入済みです。R を自力でインストールする場合は、http://www.okada.jp.org/Rwiki/などを参考 にしましょう。その後、qtl ライブライ-もインストールします。R の使い方は研究室に参考図書があります。2
します。つまり、操作をなぞる場合は赤字部分をコピー&ペーストします。入力コマンドの 各行の#以降右はコメントで、なくても操作に影響しません。
library(qtl) # qtl ライブラリーを起動する
testmap <- read.cross("csvr", file="test.csv", estimate.map=FALSE) #データ読み込み summary(testmap) #交配タイプ、個体数、マーカー数、表現型数、を確認する データが正常に読み込めているかどうか、交配タイプや個体数、マーカー数、表現型数 を確認しましょう。染色体数や染色体毎のマーカー数もみておきましょう。 2-2)欠損データの確認 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")
3
欠損データの多い2マーカーや個体は、とりあえず解析から外します。必要があれば、ジ
ェノタイピングし直すか近傍マーカーを作り直します。データの削除は以下の通りです3。
例えばデータ数が 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)
2-3)分離比の歪んだマーカーを検出する 各マーカー遺伝子座の分離比が期待値から大きくずれる場合4、遺伝子型の判別方法に問題5 がある可能性が高い。期待比からのズレは geno.table 関数を使って調べることができ、ボン フェローニ補正6した 5%水準以下のマーカーは以下のように検出できる。 gt <- geno.table(testmap) gt[gt$P.value < 0.05/totmar(testmap),] この例の場合、染色体 3 のこの領域の隣接マーカー（RM7 や RM251）の p 値も小さいた め、歪みの原因はマーカー遺伝子型の問題ではなく遺伝的な要因であるとみなせます。こ の種のマーカーの歪みは連鎖地図の作成や形質マッピングの妨げにはならないでしょう。 隣接マーカーから突出した歪みを示すマーカーは削除するのが無難です。もう一度ジェ ノタイピングしなおすか、代替マーカーを探すのがよいですが、さしあたっては解析から 外します。例えばボンフェローニ補正前の p 値が 1×10-10以下の歪みの激しいマーカーを削 除する場合には以下のようにします。
todrop <- rownames(gt[gt$P.value < 1e-10,]) testmap <- drop.markers(testmap, todrop) 2 明確な基準はありませんが、データの20％以上が欠損する場合は怪しいと考えましょう。 3 Excel で、入力ファイルを編集してもよいでしょう。 4 F2分離集団では、共優性マーカーの遺伝子型の分離がA：H：B=1:2:1 に、優性マーカーは A:C または B:D=1:3 に分離することが期待でされる。 5 マーカーが単一遺伝子座由来でない場合や、ホモ・ヘテロの判別が難しく誤った遺伝子型を含んでいる 場合などが考えられる。 6 ボンフェローニ補正は、χ2適合度検定によるp値に、調査した全マーカー数を掛ける。
4 2-4) 連鎖地図の検定
formLinkageGroups()関数で、連鎖群を推定できます。指定したmax.rf（組換価の上限）と min.lod（LODスコアの下限）を満たす任意の2つのマーカーは同じ連鎖群に割り当てます7。
lg <- formLinkageGroups(testmap, max.rf=0.45, min.lod=5) table(lg[,2])
連鎖群別に割り当てられたマーカー数が、降順に示されます。例データの場合、連鎖群 数はちょうどイネの12の染色体数に等しくなりました。その連鎖群順に8、マーカー対の組
換価とLODスコアをプロットすると連鎖群ごとのマーカー順序の精度を図示できます。 testmap <- formLinkageGroups(testmap, max.rf=0.45, min.lod=5, reorgMarkers=TRUE)
plot.rf(testmap, alternate.chrid=TRUE) 図の左上半分は組換価のプロットで、右下半分はLOD スコアのプロットです。マーカー 間で、連鎖群作成の閾値よりも組換価が小さい、もしくはLOD スコアが大きい場合を赤色 で示し（連鎖する可能性が非常に高い）、青は逆に組換価が大きいまたはLOD スコアが小 さいことを示します（独立である可能性が非常に高い）。上例では、対角線上に赤色がプロ 7 置換系統を用いたマッピングの場合は、より厳しい条件でも連鎖地図は作成できるでしょう。 8 組換価とLOD スコアのプロットは、染色体番号順ではないことに注意
5 ットされ、連鎖群の異なるマーカー間は全て青色となり、連鎖群の割り当ておよびマーカ ーの並びに矛盾がないきれいな連鎖地図であることを示しています。ちなみに、例データ で、マーカー順序をわざとシャッフルした間違った連鎖地図では以下のようなプロットに なります。
連鎖地図の検証
これまでの操作で欠損データや分離比の歪みがみられたら、入力データからおかしなマー カーを削除します。必要に応じてジェノタイピング・データのエラーチェック、再ジェノ タイピング、マーカーの差し替えを行います。連鎖群内のマーカー順序の検証
(修正)データの再読み込みtestmap <- read.cross("csvr", file="test.csv", estimate.map=FALSE) #データ再読込
例データでは染色体毎に、連鎖群が定まりマーカー順序も物理地図通りになりました（p6）。 そこで以下では染色体ごとにマーカー順序を検定します。そのさい、染色体番号（予め入力用 データの2 列目に記入しておく、test.csv を参照）を”chr”オプションで指定します。 連鎖群内のマーカー数が 30 以下程度の場合のマーカー順序の探索 orderMarkers()関数では、任意のマーカー対を取り上げ、無作為に選んだマーカーを最良の 場所に１つずつ追加します。
6 testmap <- orderMarkers(testmap, chr=12) pull.map(testmap, chr=12) ripple()関数を使うと、指定したwindowsサイズの数だけマーカー順序の並べ替えを検討し てくれます。ripple()関数の"likelihood"メソッドを用いると、尤度が一番大きいマーカー順序 を最善の並びとして検索してくれます。その際error.probを指定することで、遺伝子型デー タに含まれるエラー確率（推定値）を考慮できます。windowsサイズは大きい程（最大値は 連鎖群内のマーカー数）より緻密な検索ができます。ただし、windowsサイズを7以上9にす ると、（用いるPCの性能によっては）計算に数時間以上要する場合があります。組換価と LODスコアのプロットで、連鎖の順序に不具合が見られない場合は、windowsサイズを4～7 程度で並びを探索するとよいでしょう。
rip12lik <- ripple(testmap, chr=12, window=7, method="likelihood", error.prob=0.005) summary(rip12lik) 例の場合、2 と 1 を入れ替えた並びが 2 番目に LOD が高く、染色体長の推定値が少し長 くなることが分かります。マーカー順序はLOD が高いほどよく、染色体長はなるべく短い ものを選びます。 選ばれるマーカー順序は error.prob で指定した遺伝子型データのエラー率に影響を受けま す。エラー率を変えた場合に、最適マーカー順序が変わるか否かを検討する場合、 compareorder 関数を用いて、新たに指定したマーカー順序(new)と初期入力順序(orig)を比較 します。染色体 12 で、マーカー2 と 1 を入れ替えた順序との比較は、エラー率 0.01 で以下 のように計算します。 compareorder(testmap, chr=12, c(2,1,3:10), error.prob=0.01) compareorder(testmap, chr=12, c(2,1,3:10), error.prob=0.001) 9 オーバーナイトで計算させる場合があるとしても、せいぜいwindows サイズは 10 程度までにする
7 compareorder(testmap, chr=12, c(2,1,3:10), error.prob=0) 結果からマーカー順序は初期順序のLOD が最大で、染色体の長さも短いことがわかりま す。もし初期順序よりもLOD スコアの高いマーカーがみつかった場合は、switch.order()関 数でマーカー順序を入れ替えます。この関数は est.map()と同じ引数を使用して、マーカー 順序を交換したマーカー順序で染色体地図を推定します。
testmap <- switch.order(testmap, chr=X, c(2,1,3:10), error.prob=0.005)
例データの場合、全ての染色体で最良になるマーカー順序を初期値としています。
連鎖地図の検証
ギャップの検証 summary.map()関数を用いると、各染色体の最大ギャップがわかります。 summary.map(testmap) プロットplot.map()関数で、引数にshow.marker.names=TRUEを用いると、マーカー名付きの 連鎖地図でギャップ（連鎖地図のマーカー密度が低い箇所）を表示できます。 plot.map(testmap, show.marker.names=TRUE)8 例データの一番大きなギャップは34.0cM でした。94 個体程度以上のデータ数がある場 合、50cM 未満のマーカー・インターバルであればマーカー間の連鎖は有意といえるでしょ う。ただし、形質マッピングで、40cM 以上のマーカー間に QTL が検出された場合は、マ ーカー間（ギャップの中央付近がよい）に座乗する新たなマーカーを追加しマーカー密度を 高めた方がよいでしょう。 1 マーカー除去 ジェノタイピング・エラーを起こしやすいマーカー（例えばヘテロが区別しにくい）を検 出する方法に、droponemarker()関数を用いた 1 マーカー除去があります。一度に一つのマー カーだけ除去した場合の LOD スコアと染色体長の減少程度を観察することで、怪しいマー カーを検出することができます。
dropone <- droponemarker(testmap, error.prob=0.005) par(mfrow=c(2,1))
plot(dropone, lod=1, ylim=c(-100,0))
9 結果は横軸に染色体1~12 が連鎖地図順に並び、地図上のその位置のマーカーを取り除い たときのLOD スコアおよび染色体長の減少程度がプロットされています。取り除くと LOD スコアが上昇するマーカーや、取り除くと染色体長が大きく減少するマーカーが怪しいマー カーです。各染色体で一番染色体長を減少させるマーカーは次のように検出できます。 summary(dropone, lod.column=2) 染色体の両端に位置するマーカーは取り除くと染色体長が減少するのは当たり前なので、 これは割り引いて考える必要があります。例の場合、染色体6 と 11 で検出されたマーカー のみ染色体の内側に位置するマーカーになりました。試しにこれらのマーカーを取り除いて みる場合は以下のようにします。
badmar <- c(rownames(summary(dropone, lod.column=2))[6],rownames(summary(dropone, lod.column=2))[11]) testmap <- drop.markers(testmap, badmar)
10 newmap <- est.map(testmap, error.prob=0.005)
testmap <- replace.map(testmap, newmap)
summary.map(testmap) #連鎖地図の再推定
取り除くと遺伝距離が10cM 以上減少する染色体の内側のマーカーについては、遺伝子型 の再チェックか、必要に応じて代わりになるマーカーの遺伝子型を調べます。
問題個体の探索
countXO()関数を用いると、各個体の染色体の乗換え数（推定値）を表示できます。
plot(countXO(testmap), ylab="Number of crossovers")
集団の中で突出して染色体乗換え数が多い個体は怪しいので、削除します。例えば染色体 乗換え数が80 以上のマーカーを除いて連鎖地図をつくるには、
testmap <- subset(testmap, ind=(countXO(testmap) < 80))
とします。削除後のマーカー順序は連鎖群毎のripple 関数を用いた検証（p7）を再度行い ます。 ジェノタイピング・エラー率を見積もる 実際に調査したマーカーの遺伝子型データにはどうしてもジェノタイピングのエラーが含 まれます。以下のコードはエラー率0.001~0.0250 の範囲で連鎖地図の LOD を計算し、プ ロットしてくれます。LOD が最大になるエラー率を採用します。 loglik <- err <- c(0.001, 0.0025, 0.005, 0.0075, 0.01, 0.0125, 0.015, 0.0175, 0.02) for(i in seq(along=err)) {
cat(i, "of", length(err), "\n")
tempmap <- est.map(testmap, error.prob=err[i]) loglik[i] <- sum(sapply(tempmap, attr, "loglik")) }
11 lod <- (loglik - max(loglik))/log(10)
plot(err, lod, xlab="Genotyping error rate", xlim=c(0,0.02), ylab=expression(paste(log[10], " likelihood"))) 例データの場合はエラー率0.015 で LOD スコアが最大（付近）になりました。データに よっては、エラー率の検索範囲を広げる必要があるでしょう。 ジェノタイピング・エラーの可能性がある個体の探索 calc.errorlod()関数を用いると、ジェノタイピング・エラーの可能性がある個体を検出する ことができます。top.errorlod()関数で、エラーを持つ可能性が高いマーカーおよび個体番号 をリスト表示できます。LODスコアが4~5以上の場合、エラーの可能性が高いと考えるのが 一般的です。例えばLOD5以上のリスト表示は以下のようにします（ここでは、マーカー順 序をデタラメに並べたerror.csvを使用します）。
testmap <- read.cross("csvr", file="error.csv", estimate.map=FALSE) #データ読み込み testmap <- calc.errorlod(testmap, error.prob=0.015)
print(toperr <- top.errorlod(testmap, cutoff=4))
plot.geno()関数を用いると、検出された怪しい遺伝子型を表示できます。染色体1についてみ てみると、以下のようになります。結果の図では、左縦軸の個体番号(5,9,11)中に見出され た怪しい遺伝子型が赤色四角で示されています。横軸は染色体の連鎖地図位置(cM)です。
12
ジェノタイピング・エラーの検出では、隣接マーカーが二重組換え型のものを検出していま す（Lincoln and Lander 1992 Genomics 14:604–610)。検出されたエラー遺伝子型については、 再調査を行うべきです。error.csv の場合、全ジェノタイピング数（94 個体×138 マーカー－ 51 欠損データ）のうち LOD≧5 の閾値で検出されるジェノタイピング・エラーは 4 つのみ （0.03%）であり、連鎖地図と形質マッピングには影響を及ぼしません10。 検出されたジェノタイピング・エラーを削除する場合は次のコードが使えます。 testmap.clean <- testmap for(i in 1:nrow(toperr)) { chr <- toperr$chr[i] id <- toperr$id[i] mar <- toperr$marker[i] testmap.clean$geno[[chr]]$data[testmap$pheno$id==id, mar] <- NA } #分離比の歪みを再チェック 最後に、もう一度分離比の歪みをチェックします。
testmap <- read.cross("csvr", file="test.csv", estimate.map=FALSE) #データ再読み込み gt <- geno.table(testmap, scanone.output=TRUE)
par(mfrow=c(2,1))
plot(gt, ylab=expression(paste(-log[10], "P-value"))) plot(gt, lod=3:5, ylab="Genotype frequency")
abline(h=c(0.25, 0.5), lty=2, col="gray")
13 単一マーカーによる分離比の歪みはみられません（上段）。遺伝子型の比（下段の黒、青、 赤線は、それぞれ AA, AB, BB に対応）は、染色体 7 に１つ、染色体 8 に１つ、染色体 9 に 2 つある優性マーカーで歪んで見えますが、分離比自体は歪みません。 最終的な地図を描きます。 plot.map(testmap, show.marker.names=TRUE)
14 最終的に得られた連鎖地図のマーカー順序や染色体群および地図位置は、入力データに反 映させましょう（その際、マーカー名と遺伝子型データが入れ違いを起こさないように注意 すること）。 形質マッピングに必要な連鎖地図情報とは、染色体全域に渡り20cM 程度の遺伝距離で連 鎖するDNA マーカーの遺伝子型情報です。当然ですが、染色体断片置換系統と反復親との 交雑後代の場合は、置換領域のみの遺伝子型情報が必要になります。
