研究メモブログの目次

【RAD-seq 解析】stacks ― denovo-map.pl 編 ―

ずっと先延ばしにしていたstacks の記事を書きます。後輩からの需要の関係で、まずはリファレンス配列なしのパイプライン解析プログラムdenovo_map.pl の使い方から。

catchenlab.life.illinois.edu

 

 

 

はじめに注意です。denovo_map.pl およびリファレンスマッピング用のref_map.pl はアナコンダさんに入れてもらわずに自分で入れた場合、stacks 本体のフォルダではなく、bin フォルダに入っているので、使う前にstacks だけではなくその中のbin にもパスを通しておく必要があります。

アナコンダ経由以外の方法で入れて「あれ?動かないぞ」となったらbin にパスが通っているか確認してみてください。

アナコンダで入れる方法はいつものように

conda install -c bioconda stacks

です。パスを通したり周辺ソフトウェアを入れるのがめんどくさいので、アナコンダ経由で入れるのが個人的にはおすすめです。

 

 

一番初めにやるべきこと

インプットファイル名をソフトウェアが認識できる形式に変える

インプットファイルがペアエンドの場合、公式ページにも説明があるように、

サンプル名.1.fq.gz、サンプル名.2.fq.gz とすることで、FとRが自動で認識されます。

そのため、まずrename コマンドでフォルダ内のファイルの名前を変えます。失敗すると嫌なので別フォルダに現フォルダの中身をコピーし、名前を変換。

スクリプトは以下。Rename コマンドは何かと応用が利くので覚えておくと便利かも。

#名前の変更 rename s/変える前の名前/変えたい名前/g 変える前のファイル名

rename s/R1.fastq.gz/.1.fq.gz/g *R1.fastq.gz;

rename s/R2.fastq.gz/.2.fq.gz/g *R2.fastq.gz

 

集団情報の作成

#example (Stacks: denovo_map.pl (illinois.edu))

sample_07.01 redlake

sample_07.02 redlake

sample_06 blueriver

sample_09 blueriver

と公式サイトには書いてあるのですが、わたしがやってみたところ、長いアルファベットの集団名ではエラーが出てしまいました。結局自分が実行したときの集団名は数字で表記しました。

集団情報を作るときに便利なコマンド ls -1

これはフォルダ内のサンプルをターミナルにリストとして表示してくれるコマンドです。これでフィルタリング後の生データフォルダ内のサンプルを列で表示し、全てコピー&エクセルにペースト。置換で拡張子を一斉消去して、あとはサンプル名を一つずつに減らす。

最後に、隣の行に集団情報を並べてTab 形式で保存すれば集団情報の出来上がりです。

f:id:nemunemu_nyanko:20210216124750p:plain

実際に自分が使った集団情報データの一部


自分のデータの一部を例として貼っておきます。

 

パラメータの最適化

denovo_map.pl パイプラインで考えなければならないパラメータはnとMです。

M: スタック(仮説対立遺伝子)を仮説座に結合させるために許容されるミスマッチの数

n: カタログ構築時にスタック間で許容されたミスマッチの数(仮説遺伝子座)

これだけ見ても、できるだけ小さいほうがよさそう。

nやMのパラメーターの最適化については、色々な論文が出てますが、たとえば以下が参考になるかと。

https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12775

種特異的な値なので、参考を論文見ても結局は様々な組み合わせで試してみて、より多くのSNP が拾える数値を選ぶのが良いって結論になると思います。

ただし、nの値はMさえ決まれば決まるので(論文要参照)、無限に組み合わせを考える必要はありません。

公式サイトを見ているとほかにも共有率が閾値以下の遺伝子座を除去するオプションとかもあるみたいなのですが、自分はdenovo_map.pl のあとにstacks 内の別のプログラム(populations)やplink というソフトウェアでフィルタリングをかけていたので、ほかのオプションについてはいじっていません。他の気になる人は上の公式サイトを見てみよう。

 

denovo_map.pl の実行

denovo_map.pl -T 16 -o アウトプットフォルダ名 --paired --popmap 集団情報.txt --samples サンプルが入っているフォルダ名 -n (M-1もしくはM+1) -M (自分で決めた値)

 

パラメータの説明

-T 実行スレッド数

-o アウトプットフォルダ指定

--paired ペアエンド処理を指定

--popmap 集団情報を指定

--samples サンプルが入っているフォルダを指定

-n カタログ構築時にスタック間で許容されたミスマッチの数(仮説遺伝子座)

-M スタック(仮説対立遺伝子)を仮説座に結合させるために許容されるミスマッチの数

 

結果のフィルタリングの実行 by populations

ここからはstacks 内の別のプログラム"populations"および、stacks 以外のプログラム"plink"を使って解析します。

catchenlab.life.illinois.edu

populations は各個体の母集団を解析し、集団遺伝学的統計値を算出してくれるほか、様々な形式のアウトプットファイルを出してくれます。

www.cog-genomics.org

plink はpopulations で得た結果をさらにその後の解析にあった形に計算し直します。

【追記】今回この記事で紹介したコードではログ出力のみの働きしかしてませんでした!ごめんなさい。

 

いずれのソフトウェアもオプションがかなり多様なので、公式サイトで自分の解析にあったオプションを探しましょう。

到底書ききれないので、今回は自分が使っているオプションのみご紹介します。

まず、populations です。

自分が実行しているスクリプトは以下です。

populations -P denovo_map.pl の結果が入っているファイル --popmap 集団情報.txt -p 全集団数 -r 集団内での共有割合(0.7-0.8 を指定) --max-obs-het  0.5 --min-maf 0.01 --vcf  --plink

 

使っているオプションの説明

① データフィルタリングに関するオプション

-p 全集団の内、何集団が共有しているSNP を残すか。

-r 各集団内で何割の個体が共有しているSNP を残すか。

わたしの場合は-p オプションで全集団数を指定して、全集団が共有しているSNP のみを残し、-r オプションでさらに各集団内で特定の割合の個体が維持しているSNP のみを保持するという設定のもと、結果のフィルタリングを行っています。

--max-obs-het 0.5

最大ヘテロ接合度の指定。ハーディー・ワインベルグ平衡を仮定する場合、ヘテロ接合度は0.5 を超えない。

--min-maf 0.01

最小マイナー遺伝子頻度。読み取りエラーの可能性を考慮してわたしは1%に指定していますが、これを指定してしまうとシングルトンが除かれてしまうので、このオプションを使うかは自分が何を見たいのかを考慮して決める。

--write-random-snp データ解析を1つの遺伝子座につき1つのランダムSNPに制限。

※このオプションは次にどんな解析に使うかによって入れたり入れなかったりします。クラスタリング解析に使うアウトプットファイルを出力するときには入れます。

 

②アウトプットファイルに関するオプション

--vcf SNP を.vcf ファイルとして出力する設定。

--phylip SNP を.phy ファイルとして出力する設定。

--plink plink に使えるアウトプットファイルの書き出しを指定。

 

次に、plink でのオプション設定です。

--missing Basic statistics - PLINK 1.9 (cog-genomics.org)

欠損データの数のレポートを出してくれるオプションです。これの結果を見ると一個体あたりどのくらいのSNP がpopulations のフィルタリングで失われたかがわかります。

--allow-extra-chr Complete flag index - PLINK 1.9 (cog-genomics.org)

未認識の染色体上のデータも許容するオプション、と解釈しています。

 

【後輩向け】

わたしはここまでの一連の流れを以下のように一つのスクリプトにまとめて実行しています。

コピペ用スクリプト

denovo_map.pl -T 16 -o アウトプット用フォルダへのパス --paired --popmap 集団情報.txt --samples サンプルファイルへのパス -n 2 -M 1;

※全SNP を拾う場合のpopulations

populations -P denovo_map.pl のアウトプットファイルへのパス --popmap 集団情報.txt --max-obs-het 0.5 --fstats --fst-correction --min-maf 0.01 --vcf -p 全集団数 -r 集団内でのSNP の共有割合 --plink;

plink --file populations.plink --missing --allow-extra-chr;

mv plink* アウトプットフォルダへのパス

※SNP を一遺伝子座あたりひとつランダム抽出する場合のpopulations

※あらかじめ、denovo_map.pl のアウトプット用フォルダ名の中に、わかりやすい名前のランダムサンプリング結果専用フォルダを作っておきましょう(例:one_SNP など)

populations -P denovo_map.pl のアウトプット用ファイルへのパス -O denovo_map.pl のアウトプット用フォルダ名/ランダム抽出結果専用フォルダへのパス --popmap 集団情報.txt -p 全集団数 -r 集団内での共有割合 --max-obs-het 0.5 --min-maf 0.01 --vcf--plink --write-random-snp;

plink --file ランダム抽出結果専用フォルダへのパス/populations.plink --missing --allow-extra-chr;

mv plink* ランダム抽出結果専用フォルダへのパス

 

 

上のスクリプトテキストエディタにコピペして.sh ファイルとして保存し、細字の箇所を自分のデータに合わせて書き換えてもらえれば、あとは「sh スクリプト名.sh」で実行できるかと。

ただし、上手くいかなくても責任は負えません!ごめんね!

あと、扱うデータが大きい場合は時間がかかるのでnohup をつけて実行した方がいいかもしれません。nohup については以下の記事で読めます。

nemunemu-nyanko.hatenablog.com

 

【追記履歴】

20200219 最後のスクリプトの改行、アウトプット及びインプットを修正。集団情報の作り方の修正。

20210222 最後のスクリプトのインプットの修正。