読者です 読者をやめる 読者になる 読者になる

マクロ生物学徒の備忘録

某国立大で進化生態学を学ぶポンコツ研究者の備忘録

plinkでペアワイズ連鎖不平衡解析|非モデル生物の場合

はじめに

RAD-seqやMIG-seqなどで取得できるSNP情報は、もろもろの解析に使用する前に連鎖不平衡を調べて除去してやる必要がある。

ここではplinkを使用してSNP間の相関係数をペアワイズで計算し、閾値(こちらで設定できる)よりも相関係数が大きな組み合わせを連鎖不平衡のあるSNPとして検出する。

また、今回は非モデル生物でRAD-seqやMIG-seqをおこなってde novoでアライメントした場合を想定し、各locusの染色体番号や位置情報を使わずに連鎖解析をおこなう。

plinkのインストール方法(mac)はこちら

 

インプットファイル

.pedファイルと.mapファイルが必要

.pedファイルは.strなどからPGDspiderで変換可能

ただし、同じ座位に3つ以上アリルがあるとはじかれるので、事前に取り除いておく必要がある

.mapファイルはエクセルで自分で作れる

左から順に

染色体番号|座位の番号|遺伝的距離|ポジション

を半角スペース区切りで記述していくが、de novoの場合、座位の番号以外は0と記入すれば良い

0 loc1 0 0
0 loc2 0 0
0 loc3 0 0
0 loc4 0 0

といった感じ

.pedファイルに記述されている座位と.mapファイルに記述されている座位は一致する必要があるかも

 

実行

.pedファイルと.mapファイルを同じファイル名で(拡張子は異なる)同じディレクトリに格納する

例えば、

mkdir ~/Desktop/plink

で、できたフォルダに

○○.ped
○○.map

を移すといった感じ

その後、

plink --file ○○ --r2 --ld-window 1000000 --ld-window-r2 0

※〇〇:拡張子をつけてはいけない

と打ち込めば、同ディレクトリに

plink.ld

が作成される

コマンドのオプションを解説すると、

-- file ○○ インプットファイルの指定

-- r2 相関係数の計算

-- ld-window 1000000 何個隣のSNPまで相関係数を計算するか …※

--ld-window-r2 0 相関係数が指定した数値以上のものを出力する(全て出力するときは0)

※入力したデータが実際の配列順で並んでいる場合、あまり遠くの座位と相関係数を出しても意味が無い。そのため、デフォルトでは9個隣の座位までしか相関係数を計算しないようになっている。de novoの場合、実際の配列順に並んでいるわけではないので、全てのSNPを総当りで調べる必要がある。よって、ここには手持ちのSNP数以上の数字を入れる必要がある。SNP数以上であれば数字はなんでもよい。

 

アウトプットファイル

plink.ldをmiでもエクセルでもなんでもいいので開いてみると、

CHR_A   BP_A   SNP_A   CHR_B   BP_B   SNP_B           R2 
     0              0        loc1           0               0         loc2       0.0152985
     0              0        loc1           0               0         loc3       0.0515134
     0              0        loc1           0               0         loc4                1

こんな感じになっているはず

見たままだが、loc1とloc2の相関係数が0.0152985であるという感じに読む

ちなみに、たとえば

--ld-window-r2 0

の0の数字を0.2にすれば、相関係数が0.2以上のものだけをリストアップしてくれる