Transdecoderをできるだけ自動化する
コマンドラインでは、実行手順を書き込んだファイル(シェルスクリプト)を用意することで、多くの作業を自動化できる。自動化のための作業手順
必要なフォルダーを予め作っておく
:> mkdir Trandecoder
コマンドを書き込んだスクリプトファイルをフォルダーに置いて、スクリプトを書いておく
最長コード領域の抽出
:> touch transdecoder_long.sh
:> echo 'TransDecoder.LongOrfs -t ../Trinity.fasta' >> transdecoder_long.sh
Blast検索
:> touch transdecoder_blastp.sh
:> echo 'blastp -query ./Trinity.fasta.transdecoder_dir/longest_orfs.pep -db uniprot_sprot -max_target_seqs 1 -outfmt 6 -evalue 1e-10 -num_threads 15 -out blastp.outfmt6' >> transdecoder_blastp.sh
(num_threads の数字は、実際に自分が使うコンピューターのCPU数-1を指定すると最も効率が良い)
Pfam検索
:> touch transdecoder_pfam.sh
:> echo 'hmmscan --cpu 16 -E 1e-10 --domE 1e-10 -o pfam.out --domtblout pfam.out_domtblout /Volumes/DATA_HD/Pfam32/Pfam-A.hmm ./Trinity.fasta.transdecoder_dir/longest_orfs.pep' >> transdecoder_pfam.sh
(cpu の数は、実際に自分が使うコンピューターのCPU数を指定)
(Pfam-A.hmm の場所は、自分が置いた場所のパスを指定)
コード領域の推定
:> touch transdecoder_pred.sh
:> echo 'TransDecoder.Predict -t ../Trinity.fasta --retain_pfam_hits ./pfam.out_domtblout --retain_blastp_hits ./blastp.outfmt6 --cpu 0' >> transdecoder_pred.sh
(cpu を0に設定すると、自動で使える全部のCPUを使う)
推定されたコード領域を使って、Pfam検索
:> touch Pfam_search.sh
:> echo 'hmmscan --cpu 16 -E 1e-10 --domE 1e-10 -o pfam_Final.out --domtblout pfam_Final.out_domtblout /Volumes/DATA_HD/Pfam32/Pfam-A.hmm ./Trinity.fasta.transdecoder.pep' >> Pfam_search.sh
(cpu の数は、実際に自分が使うコンピューターのCPU数を指定)
(Pfam-A.hmm の場所は、自分が置いた場所のパスを指定)
全部を実行するスクリプトを作る
:> touch analyze.sh
:> echo '/bin/sh ./transdecoder_long.sh && /bin/sh ./transdecoder_pfam.sh && touch hmmer.finished && /bin/sh ./transdecoder_blastp.sh && touch blastp.finished && /bin/sh ./transdecoder_pred.sh && touch pred.finished && /bin/sh ./Pfam_search.sh' >> analyze.sh
(どこまで進んだかわかりにくいので、作業完了ごとにtouchコマンドで空ファイルを作っている)
シェルスクリプトでは、コマンドを「&&」でつなぐ事で、完了したら次のコマンドの実行を自動化できる。必要な作業を全部「&&」でつないでおくと、様々な作業を自動化できる。
スクリプトを実行する
Transdecoder フォルダーをTrinityの結果フォルダーにコピーする
:> cd /パス/Transdecoder
:> /bin/sh ./analyze.sh
スクリプトファイル(.sh という拡張子をつけたファイル)を作成する意味
スクリプトファイルを作成して作業したフォルダー内に残しておくことで、作業ログにもなる。そのフォルダー内に残された解析結果が、どのようなコマンドによって得られたものなのか確実な記録になるので、完璧な再現性を保証する研究ノート代わりにもなる。紙のノートに手書きする際に発生する転記ミスを防ぐ意味もあるので、解析研究での研究ノートの残し方としてお薦めする。
スクリプトは「#」から改行コードまでをメモとして扱い、コマンドには影響しないという仕組みがあるので、上記のスクリプ内に実行した日付や研究の条件などを書き込んでおくと、研究ノートとしての情報量が上がってより有益である。
例: 上記 analyze.sh の中に、下記のようにしてメモを書き加えておく
# 2020/05/05
# ナミアゲハのRNA-seq
# Trinity v2.10.0 を使用
目次
- 環境構築
- リードのクオリティチェックとアセンブル
- Transdecoder による遺伝子予測とアノテーション
- リード数のカウント
- カウントマトリクスの作成
- リードカウントのQC解析
- DE解析
- DEG取り出し
- おまけ1: Transdecoderの自動化
- おまけ2: 発現量解析の自動化
- おまけ3: Anacondaコマンドの使い方一覧