研究メモブログの目次

【疑問】一遺伝子座あたり1SNPをランダム抽出するには?(解決済み)

ipyrad かなりいい!stacks よりSNP たくさん拾える!と喜んでいたのも束の間、タイトル通りの疑問が発生してしまいました。
ipyrad にもランダムサンプリングオプションはあるのですが、ランダムサンプリングで出せるのは.phy 形式とstructure 形式だけなんですよね。
わたしはSNP コールの後、Tassel5 を使って個体のフィルタリングを行っているのですが、structure ファイルではフィルタリングができない。
さらに.phy のランダムサンプリングファイルは出力できず…(原因不明)。
そういうことなら、別のソフトウェアで全SNP ファイルをいじるのが早そう!となったわけです。
ひとまず、調べるしかありません。


【解1】
全SNP からランダム抽出するのはvcftoolsで出来そうです。
speciationgenomics.github.io
ただし、これはstacks やipyrad のものとは違って各遺伝子座から一つ抽出する、ということはできなさそう…?


【解2】
先輩からstacks のpopulations はvcf をインプットに使えるから、ipyrad で出したvcf をpopulations で処理してみたら?とのアドバイスをいただき、試してみました。
結果「できなかった」
同じvcf でもstacks のvcf のフォーマットは若干ipyrad で出したものと違うみたい…。
自力でテキストエディタでいじってみましたが、どうにもうまくいきませんでした。無念。

【解3】完全にランダムサンプリングではないですが、1遺伝子座あたり1SNP をフィルタリングしてくれるプログラムは見つかりました。
github.com

使い方(参考:Filtering your SNPs | 2019-06-11_GBS_EE

①vcftools とPLINK とR 、およびR パッケージのdplyr/readr/tibble/stringr をインストールしておく

github からGBS_SNP_filter フォルダをダウンロード$ git clone https://github.com/laninsky/GBS_SNP_filter.git

➂フィルタリングしたいvcf ファイルをダウンロードしたフォルダにコピー
cp ../working/denovo/stacks/populations.snps.vcf GBS_SNP_filter

④popmap.txt という名前で集団情報ファイルを作ります
作り方はStacks のpopulations の集団情報と同じで良さげ
作り方はこれを見てね:【RAD-seq 解析】stacks ― denovo-map.pl 編 ― - 研究メモ
このファイルも同じフォルダに入れておきます

⑤GBS_SNP_filter フォルダに入ります

⑥GBS_SNP_filter.txt というファイルをテキストエディタで作成します

GBS_SNP_filter の中身
1.vcf ファイルの名前
2.SNP欠損割合の閾値 (e.g. 0.85 なら85% かそれ以上の個体が持つSNP のみになる)
3.サンプルの欠損割合の閾値 (e.g. 0.9 なら最大90%が欠損値でもそのサンプルは残る)
4.HWEから外れているかどうかを判断するために使用されるp値(例:0.05、大きいほど多くの遺伝子を削除)
5.互いにLDにあるかどうかを判断するために使用されるr^2カットオフ (例: 0.5, 小さい値ほど多くの遺伝子を削除)
6.ある遺伝子座において何個体群がHWE/in LDから外れたらその遺伝子座を削除するか(例:3、小さい値ほど多くの遺伝子座を削除)
7.vcfファイル内のカラムヘッダに遺伝子座ID
8.遺伝子座IDの正規表現パターンを指定 (何だかわかりません)※指定の必要がない場合も、空白行を残しておく必要がある

⑦実行
bash GBS_SNP_filter.sh


で、できるとチュートリアルには書いてあるのですが、わたしはできませんでした(なんで!)


【解4】
ようやく見つけました、ちゃんと動くランダムサンプリング用スクリプトgithub.com

このスクリプトを目的vcf ファイルがあるフォルダに保存して、あとは
python vcf_single_snp.py Hoge.vcf > out.vcf
を実行するだけです。

は~~~つかれた!