pyRADの使い方
はじめに
pyRADって何？どうやってインストールするの？って方はこちら↓
pyRADのインストール｜OS X ElCapitan (10.11.6) - マクロ生物学徒の備忘録
ここではインストールが済んでいる前提で書きます。
pyRADの使い方は公式サイトのチュートリアルにわりとわかりやすく書いてある。
読めばだいたいわかるはず。
解釈に自信がないところもあるので、間違ってたらご指摘ください。
（間違っていても責任はとれませんのであしからず。）
大まかな流れは、インプットファイルを用意して、pyRADでパラメータファイルを作成し、パラメータファイルをテキストエディタで改変し、pyRADを実行する、という感じ。
この項では、
・インプットファイルについて
・パラメータファイルの作り方
・パラメータの各項目について
・pyRADの実行
・パラメータをいじってみた印象
について触れていく。
インプットファイルについて
fastq形式で、
・fastqファイルが個体ごとにわかれている状態
・fastqファイルが１つにまとまっており、個体識別用のバーコードが付いている状態
・圧縮されて〇〇.fastq.gzとなっている状態
のどれでもできる。
ちなみに、次世代シーケンサから返ってきたファイルをそのまま使ってもいいが、FASTX-toolkitなどを使って、処理しておいた方がSNPが多くとれてくる。そんな気がする。
FASTX-toolkitのインストール｜OS X ElCapitan (10.11.6) - マクロ生物学徒の備忘録
一般に、次世代シーケンサが出力するデータは、リードの最後の方のクオリティが低くなる傾向がある。
FASTX-toolkitを使うと、後ろ側からクオリティが低い部分を削る処理が可能なので、これだけでもやっておくのがオススメ。
パラメータファイルの作り方
pythonが立ち上がっているとして、
cd ~/Desktop/MIG/pyrad
python ~/Analysis/pyrad/pyrad/pyRAD.py -n ##pyRADをインストールしたディレクトリを指定する→params.txtが作られるはず。
このparams.txtを好きなテキストエディタで開くと、
==** parameter inputs for pyRAD version 3.0.66 **======================== affected step ==
./ ## 1. Working directory (all)
./*.fastq.gz ## 2. Loc. of non-demultiplexed files (if not line 18) (s1)
./*.barcodes ## 3. Loc. of barcode file (if not line 18) (s1)
vsearch ## 4. command (or path) to call vsearch (or usearch) (s3,s6)
muscle ## 5. command (or path) to call muscle (s3,s7)
TGCAG ## 6. Restriction overhang (e.g., C|TGCAG -> TGCAG) (s1,s2)
2 ## 7. N processors (parallel) (all)
6 ## 8. Mindepth: min coverage for a cluster (s4,s5)
4 ## 9. NQual: max # sites with qual < 20 (or see line 20)(s2)
.88 ## 10. Wclust: clustering threshold as a decimal (s3,s6)
rad ## 11. Datatype: rad,gbs,pairgbs,pairddrad,(others:see docs)(all)
4 ## 12. MinCov: min samples in a final locus (s7)
3 ## 13. MaxSH: max inds with shared hetero site (s7)
c88d6m4p3 ## 14. Prefix name for final output (no spaces) (s7)
==== optional params below this line =================================== affected step ==
## 15.opt.: select subset (prefix* only selector) (s2-s7)
## 16.opt.: add-on (outgroup) taxa (list or prefix*) (s6,s7)
## 17.opt.: exclude taxa (list or prefix*) (s7)
## 18.opt.: loc. of de-multiplexed data (s2)
## 19.opt.: maxM: N mismatches in barcodes (def= 1) (s1)
## 20.opt.: phred Qscore offset (def= 33) (s2)
## 21.opt.: filter: def=0=NQual 1=NQual+adapters. 2=strict (s2)
## 22.opt.: a priori E,H (def= 0.001,0.01, if not estimated) (s5)
## 23.opt.: maxN: max Ns in a cons seq (def=5) (s5)
## 24.opt.: maxH: max heterozyg. sites in cons seq (def=5) (s5)
## 25.opt.: ploidy: max alleles in cons seq (def=2;see docs) (s4,s5)
## 26.opt.: maxSNPs: (def=100). Paired (def=100,100) (s7)
## 27.opt.: maxIndels: within-clust,across-clust (def. 3,99) (s3,s7)
## 28.opt.: random number seed (def. 112233) (s3,s6,s7)
## 29.opt.: trim overhang left,right on final loci, def(0,0) (s7)
## 30.opt.: output formats: p,n,a,s,v,u,t,m,k,g,* (see docs) (s7)
## 31.opt.: maj. base call at depth>x<mindepth (def.x=mindepth) (s5)
## 32.opt.: keep trimmed reads (def=0). Enter min length. (s2)
## 33.opt.: max stack size (int), def= max(500,mean+2*SD) (s3)
## 34.opt.: minDerep: exclude dereps with <= N copies, def=1 (s3)
## 35.opt.: use hierarchical clustering (def.=0, 1=yes) (s6)
## 36.opt.: repeat masking (def.=1='dust' method, 0=no) (s3,s6)
## 37.opt.: vsearch max threads per job (def.=6; see docs) (s3,s6)
==== optional: list group/clade assignments below this line (see docs) ==================こうなっているはず。
##14までは必須で、それ以降は任意。
集団の情報を入力するには「optional: list group/clade assignments below this line」以下に記述する。
パラメータの各項目について
まだよく理解しきれていないところもある。
そのへんはご勘弁を。
いろいろ自分でいじってみてどうなったかは、あとの章で述べる。
## 1. Working directory
作業ディレクトリを指定する。
ここにアウトプットファイルが出力される。
相対パス、絶対パスのどちらでもいい。
## 2. Loc. of non-demultiplexed files (if not line 18)
インプットファイルのパスを指定。
複数のファイルを使用したいときは、*（ワイルドカード）も使える。
ただし、バーコードなどがついたままのrawデータ（.fastq / .fastq.gz）に限る。
バーコードが取れたあとのファイルを使用したい場合は、ここは空欄にする。
## 3. Loc. of barcode file (if not line 18)
バーコードファイルのパスを指定。
バーコードが取れたあとのファイルを使用したい場合は、ここも空欄にする。
## 4. command (or path) to call vsearch (or usearch)
vsearchのコマンドを指定する。
パスが通っている場合は
vsearch
と書くだけ。
相対パスか絶対パスを記述することもできる。
## 5. command (or path) to call muscle
muscleのコマンドを指定する。
muscleはダウンロードしたままの状態だと、
muscle3.8.31_i86darwin64
みたいな感じになっている。
muscle3.8.31_i86darwin64 → muscle にリネームしていれば、
muscle
だけでよい。
相対パスか絶対パスを記述することもできる。
## 6. Restriction overhang (e.g., C|TGCAG -> TGCAG)
RAD-seqの場合、使用した制限酵素サイトを指定する。
すでに制限酵素サイトがfastqファイルから除かれいるなら空欄でいいはず。
MIG-seqなど、制限酵素を使用していない場合も空欄でいい。
## 7. N processors (parallel)
どれだけパラレルに作業をするか。
２にすると２つ同時進行、３にすると３つ同時進行で作業をする。
多けりゃいいってもんではない。
使用するコンピューターのコア数と同じにするといいとか。
## 8. Mindepth: min coverage for a cluster
リードを採用する最低カバレッジ数。
数字を大きくすれば、質の高くてSNP数の少ないデータに、
小さくすれば、質が低くてSNP数の多いデータになる。
個人的には最低10くらいはほしい気がする。
## 9. NQual: max # sites with qual < 20 (or see line 20)
１つのリードにつき、シーケンスできていないサイト（N）をいくつまで許すか。
short readの場合と、long readの場合で、適正値が違う。
FASTX-toolkitで後ろから削る処理をしていない場合、この処理で大幅に使えるリードが減ってしまう。たぶん。
## 10. Wclust: clustering threshold as a decimal
何%一致した配列をクラスタリングするか。
例えば、90%にしたい場合
.90
という感じに記述する。
## 11. Datatype: rad,gbs,pairgbs,pairddrad,(others:see docs)
rad, ddrad, gbs, pairddrad, pairgbs, mergeが使用できるらしい。
どう違うかはよく知らない。
## 12. MinCov: min samples in a final locus
何個体以上で共通して読めたリードを採用するか。
数字を大きくすれば、質の高くてSNP数の少ないデータに、
小さくすれば、質が低くてSNP数の多いデータになる。
## 13. MaxSH: max inds with shared hetero site
これが未だによくわからない。
ヘテロの個体があまりにも多い場合、その座位はSNPじゃなくてパラログが同じクラスターにまとめられているだけという可能性が高い。
何個体以上でパラログに認定するか、という閾値をここで設定する、かな？たぶん。
## 14. Prefix name for final output (no spaces)
アウトプットファイルの名前。
なんでもよい。
以下はオプションの設定。
だいたい空欄でいいはず。
めぼしい項目だけ解説する。
## 18.opt.: loc. of de-multiplexed data
バーコードが取れたあとの個体ごとに分かれたfastqファイルを使用する場合、ここで指定する。
## 2と同様、ワイルドカードが使える。
## 20.opt.: phred Qscore offset (def= 33)
クオリティフィルターの閾値。
## 30.opt.: output formats: p,n,a,s,v,u,t,m,k,g,* (see docs)
どのアウトプットファイルを出力するか。
とりあえず「*」としておけば全部出力してくれるのでオススメ。
最後に集団の情報を記述する。
例えば、インプットファイルが
aaa001.fastq aaa002.fastq …
bbb001.fastq bbb002.fastq …
ccc001.fastq ccc002.fastq …
※aaa、bbb、cccは集団の名前
みたいなっている場合、次のようにする。
==== optional: list group/clade assignments below this line (see docs) ==================
1 4 aaa*
2 4 bbb*
3 4 ccc*一列目の数字は集団の番号。
集団の名前は入れられず、123とかの数字か、ABCとかのアルファベットのみ使用できるらしい。
二列目の数字はよくわからない。
これ以下の個体数しか得られなかった集団は放棄するってことかな？
三列目にファイル名を記入する。
例のようにワイルドカードを使ってもいいし、カンマ（,）区切りで並べていっても良い。
また、各列はタブで区切る。
pyRADの実行
pythonが立ち上がっている前提。
cd ~/Desktop/MIG/pyrad
python ~/Analysis/pyrad/pyrad/pyRAD.py -p params.txt ##pyRADをインストールしたディレクトリを指定するこれで動くはず。
コンピューターの性能やデータサイズにもよるが、普通のパソコンでやる場合、１〜３日くらいかかる。
ちなみに、pyRADは７つのステップに分かれている。
params.txtの(s1)、(s3,s6)、(all)といった表記はそのパラメータがどのステップに関わっているかを示したものである。
ステップ７に関わるパラメータだけをいじってやり直したいときは
-s 7
というオプションを追加するとステップ７だけやり直すことができる。便利。
パラメータをいじってみた印象
SNP数を増やしたかったり、質を上げたかったりする場合、いじれる項目は、
リードを採用する最低カバレッジ数（## 8）
何個体以上で共通して読めたリードを採用するか（##12）
がメインといったところか。
なので、このへんをいじった感想を書く。今後も追記するかも。しないかも。
あくまでも自分の場合なので、参考程度に。
・リードを採用する最低カバレッジ数（## 8）
3→10にするとSNP数が0.75倍くらいになった。
思ったよりは減らないんだなー、という印象。
3と10だと質が大きく違う（気がする）ので、わりとお得感があった。
・何個体以上で共通して読めたリードを採用するか（## 12）
全個体数の75%→50%にするとSNP数が1.5倍くらいになった。
このパラメータはSNP数にわりと大きく影響するなー、という印象。
・その他
## 13をいじっても大きくSNP数が変動することはなかった。
## 5はあまりいじってない。それよりはFASTX-toolkitでの処理のパラメータいじったほうがいいかも。
FASTX-toolkitで後ろから削ってやると、削らない場合よりもSNP数が多いという結果に。うまく削ってやると1.3倍くらいになった。
