62.9K Views
December 30, 21
スライド概要
mRNA-Seq のデータを De novo で解析する方法をがちゃがちゃまとめたもの。
社会人になり、アフリカに行き、そして大学に戻ってきたピチピチの博士課程学生(専門:環境科学)である。博士号はまだない。バイオインフォマティクスの諸々についてふわっふわっと説明したい。
mRNA-Seq 解析の流れをざっくりと説明してみた mRNA-Seq 解析 De novo 編 2021/12/06 ⽔産⽣物環境学(九州⼤学) ⾼井優⽣
mRNA-Seq の de novo 解析 de novo 解析を決意するまでの流れ しっかりしたゲノム配列&遺伝⼦情報ファイルがある (いわゆるモデル⽣物である) うむ 前回までの解析⽅法で 頑張ってください (STAR からの edgeR) ない リファレンスに使えそうな 近縁種のゲノム配列もしくは ドラフトゲノム配列がある ある そんなものはない Trinity で Trinity の 頑張りましょう ゲノムガイドアセンブリを 結局こうなる時も 試してみましょう あります とりあえず Trinity が必要です
Trinity の インストール Trinity
Trinity のインストール Trinity(v2.13.2) Trinity の公式 HP に従ってぽちぽちします ubuntu@ubuntu-man$ wget https://github.com/trinityrnaseq/trinityrnaseq/releases/download/Trinityv2.13.2/trinityrnaseq-v2.13.2.FULL.tar.gz ubuntu@ubuntu-man$ tar xvf trinityrnaseq-v2.13.2.FULL.tar.gz ubuntu@ubuntu-man$ cd trinityrnaseq-v2.13.2 ubuntu@ubuntu-man$ make Using gnu compiler for Inchworm and Chrysalis cd Inchworm && make make[1]: Entering directory '/home/ubuntu/src/trinityrnaseq-v2.13.2/Inchworm' mkdir -p build : : ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Performing Unit Tests of Build ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Inchworm: has been Installed Properly Chrysalis: has been Installed Properly QuantifyGraph: has been Installed Properly GraphFromFasta: has been Installed Properly ReadsToTranscripts: has been Installed Properly parafly: has been Installed Properly インストール作業は もうちょっと続きます
Trinity のインストール Trinity(v2.13.2) Trinity の公式 HP に従ってぽちぽちします ubuntu@ubuntu-man$ make plugins cd trinity-plugins && make plugins make[1]: Entering directory '/home/ubuntu/src/trinityrnaseq-v2.13.2/trinity-plugins' cd slclust && make CXX=g++ CC=gcc install : : ## Checking plugin installations: slclust: collectl: has been Installed Properly has been Installed Properly これでインストールは終わり
Trinity の実⾏テスト Trinity(v2.13.2) Trinity の公式 HP に従って実⾏テストをしてみます ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh
Trinity の実⾏テスト Trinity(v2.13.2) Trinity の公式 HP に従って実⾏テストをしてみます ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Error, cannot find jellyfish installed on this system. Be sure to install it. You can get it here: http://www.genome.umd.edu/jellyfish.html at /home/ubuntu/bin/Trinity line 4019.
Trinity の実⾏テスト Trinity(v2.13.2) Trinity の公式 HP に従って実⾏テストをしてみます ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Error, cannot find jellyfish installed on this system. Be sure to install it. You can get it here: http://www.genome.umd.edu/jellyfish.html at /home/ubuntu/bin/Trinity line 4019. エラーになってしまいました jellyfish というソフトウェアがない︕と⾔っているのでぐぐって⼊れます ubuntu@ubuntu-man$ sudo apt install jellyfish
Trinity の実⾏テスト Trinity(v2.13.2) Trinity の公式 HP に従って実⾏テストをしてみます ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Error, cannot find jellyfish installed on this system. Be sure to install it. You can get it here: http://www.genome.umd.edu/jellyfish.html at /home/ubuntu/bin/Trinity line 4019. エラーになってしまいました jellyfish というソフトウェアがない︕と⾔っているのでぐぐって⼊れます ubuntu@ubuntu-man$ sudo apt install jellyfish もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh
Trinity の実⾏テスト Trinity(v2.13.2) Trinity の公式 HP に従って実⾏テストをしてみます ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Error, cannot find jellyfish installed on this system. Be sure to install it. You can get it here: http://www.genome.umd.edu/jellyfish.html at /home/ubuntu/bin/Trinity line 4019. エラーになってしまいました jellyfish というソフトウェアがない︕と⾔っているのでぐぐって⼊れます ubuntu@ubuntu-man$ sudo apt install jellyfish もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : は︖またエラーなんやけど︖ : Error, cannot find path to bowtie2 () or bowtie2-build (), which is now needed as part of Chrysalis' read scaffolding step. If you should choose to not run bowtie, include the --no_bowtie in your Trinity command.
Trinity の実⾏テスト Trinity(v2.13.2) bowtie2 をぐぐって⼊れます ubuntu@ubuntu-man$ linux-x86_64.zip ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.4/bowtie2-2.4.4unzip bowtie2-2.4.4-linux-x86_64.zip ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2 /home/ubuntu/bin/ ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2-build /home/ubuntu/bin/ もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd /home/ubuntu/src/trinityrnaseq-v2.13.2/sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh
Trinity の実⾏テスト Trinity(v2.13.2) bowtie2 をぐぐって⼊れます ubuntu@ubuntu-man$ linux-x86_64.zip ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.4/bowtie2-2.4.4unzip bowtie2-2.4.4-linux-x86_64.zip ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2 /home/ubuntu/bin/ ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2-build /home/ubuntu/bin/ もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd /home/ubuntu/src/trinityrnaseq-v2.13.2/sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Trinity Trinity-v2.13.2 requires salmon to be installed. Get it here: https://combinelab.github.io/salmon/ at /home/ubuntu/bin/Trinity line 4057.
Trinity の実⾏テスト Trinity(v2.13.2) あああああ あああ︕︕︕ bowtie2 をぐぐって⼊れます ubuntu@ubuntu-man$ linux-x86_64.zip ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ ubuntu@ubuntu-man$ wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.4/bowtie2-2.4.4unzip bowtie2-2.4.4-linux-x86_64.zip ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2 /home/ubuntu/bin/ ln -s /home/ubuntu/src/bowtie2-2.4.4-linux-x86_64/bowtie2-build /home/ubuntu/bin/ もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd /home/ubuntu/src/trinityrnaseq-v2.13.2/sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : Trinity Trinity-v2.13.2 requires salmon to be installed. Get it here: https://combinelab.github.io/salmon/ at /home/ubuntu/bin/Trinity line 4057.
Trinity の実⾏テスト Trinity(v2.13.2) salmon をぐぐって⼊れます ubuntu@ubuntu-man$ wget https://github.com/COMBINE-lab/salmon/releases/download/v1.6.0/salmon1.6.0_linux_x86_64.tar.gz ubuntu@ubuntu-man$ tar xvf salmon-1.6.0_linux_x86_64.tar.gz ubuntu@ubuntu-man$ ln -s /home/ubuntu/src/salmon-1.6.0_linux_x86_64/bin/salmon /home/ubuntu/bin/ もう⼀回テストします(これで実⾏できるはず) ubuntu@ubuntu-man$ cd /home/ubuntu/src/trinityrnaseq-v2.13.2/sample_data/test_Trinity_Assembly/ ubuntu@ubuntu-man$ ./runMe.sh : : ##### Done Running Trinity ##### できました
Trinity の実⾏テスト Trinity(v2.13.2) ここまで「jellyfish がない」「bowtie2 がない」「salmon がない」というエラーが⽴て続き 「あああああああ︕なんでなん︕︕︕」とお⾒苦しい姿をお⾒せしてしまいましたが インストールマニュアルより抜粋 しっかりとマニュアルに書いていました(僕の環境では samtools はインストール済みでした) はい、しっかりと マニュアルを読みましょう
リアルガチの De novo 解析 Trinity
リアルガチの de novo 解析とは 使⽤するソフトウェア︓Trinity(v2.13.2) しっかりしたゲノム配列&遺伝⼦情報ファイルがある (いわゆるモデル⽣物である) うむ 前回までの解析⽅法で 頑張ってください (STAR からの edgeR) ない リファレンスに使えそうな 近縁種のゲノム配列もしくは ドラフトゲノム配列がある ある そんなものはない Trinity の ゲノムガイドアセンブリを 頼りにできるものは 試してみましょう ⼰が保有する mRNA-Seq のリードのみ︕ という場合の解析⽅法です Trinity で 頑張りましょう
リアルガチの de novo 解析とは 使⽤するソフトウェア︓Trinity(v2.13.2) mRNA-Seq のリードは断⽚化した mRNA(cDNA)配列なわけで うまいことリードを繋ぎ合わせていけば元々の mRNA 配列を 復元できるのである リードデータ というわけで こういう⾵に 重なってるところを うまいこと使って 復元する 鋳型になった mRNA (cDNA) リードデータから復元した mRNA (cDNA) 配列 鋳型になった mRNA (cDNA) どうにかこうにか頑張って mRNA の配列を復元して それをリファレンスにしちゃえばいいじゃない というのが de novo 解析で、Trinity はこれをやってくれる
リアルガチの de novo 解析とは 使⽤するソフトウェア︓Trinity(v2.13.2) mRNA-Seq のリードは断⽚化した mRNA(cDNA)配列なわけで うまいことリードを繋ぎ合わせていけば元々の mRNA 配列を 復元できるのである コンティグ リードデータ というわけで こういう⾵に 重なってるところを うまいこと使って 復元する 鋳型になった mRNA (cDNA) リードデータから復元した mRNA (cDNA) 配列 鋳型になった mRNA (cDNA) どうにかこうにか頑張って mRNA の配列を復元して それをリファレンスにしちゃえばいいじゃない というのが de novo 解析で、Trinity はこれをやってくれる
サンプルファイルの準備 使⽤するソフトウェア︓Trinity(2.13.2) Trinity のマニュアルに従ってサンプルファイル(タブ区切り*1)を準備する ⼀列⽬に実験条件、⼆列⽬にサンプル名、三列⽬と四列⽬にサンプルのリードファイル名*2を書けばいいらしい ubuntu@ubuntu-man$ nano sample.tsv ubuntu@ubuntu-man$ cat sample.tsv Control Con_1 cleaned_C1_1.fq.gz Control Con_2 cleaned_C2_1.fq.gz Control Con_5 cleaned_C5_1.fq.gz EPS EPS_1 cleaned_E1_1.fq.gz EPS EPS_4 cleaned_E4_1.fq.gz EPS EPS_5 cleaned_E5_1.fq.gz *1: cleaned_C1_2.fq.gz cleaned_C2_2.fq.gz cleaned_C5_2.fq.gz cleaned_E1_2.fq.gz cleaned_E4_2.fq.gz cleaned_E5_2.fq.gz タブ区切りっていうのはカンマとかスペースじゃなくてタブでデータを区切ること *2: リードファイルとサンプルファイルの格納場所が違う場合はサンプルファイルの絶対パスを書くように
Trinity の実⾏ 使⽤するソフトウェア︓Trinity(2.13.2) 以下のオプションを使⽤して Trinity を実⾏します(詳細はマニュアルを参照してください) オプション 指定する内容の⼤まかな説明 --seqType リードファイルの種類(fastq ファイルの場合は fq) --max_memory 解析に使⽤するメモリの最⼤量(でかい⽅が良いです) --samples_file サンプルファイルの場所 --CPU 使⽤するスレッド数 --output 解析結果の出⼒先 (trinity という⽂字列を含んでいた⽅が良いらしいです) こんな感じです ubuntu@ubuntu-man$ Trinity --seqType fq --max_memory 128G --samples_file sample.tsv \ > --CPU 8 --output /home/ubuntu/takai/trinity とりあえずの作業としてはこれで終わりです けっこう時間かかるので気⻑に待ちましょう(使⽤するデータ、解析環境によっては数⽇かかります)
Trinity の実⾏結果 使⽤するソフトウェア︓Trinity(2.13.2) ⻑い⻑い待ち時間が終わると以下 2 つのファイルができます ファイル名 ファイルの中⾝の⼤まかな説明 <output>.Trinity.fasta Trinity が構築したコンティグ配列 <output>.Trinity.fasta.gene_trans_map 遺伝⼦と転写産物の対応表 <output>.Trinity.fasta をリファレンスにしてリードのマッピングを⾏います
コンティグ配列へのマッピング 使⽤するソフトウェア︓Trinity(2.13.2) Trinity はマッピングするためのスクリプトも準備してくれているので、それを使います マッピング⽤スクリプトは Trinity をインストールしたディレクトリの util ディレクトリに⼊っています 使⽤するスクリプトは align_and_estimate_abundance.pl という名前です(詳細はマニュアル参照) オプション 指定する内容の⼤まかな説明 --transcripts コンティグ配列の場所(<output>.Trinity.fasta の場所) --seqType リードファイルの種類 --samples_file サンプルファイルの場所(先ほど作成したサンプルファイルを使いまわします) --est_method 発現量予測の⼿法(RSEM)*マニュアル参照 --aln_method マッピングに使⽤するツール(bowtie2) --output_dir 解析結果の出⼒先 --gene_trans_map <output>.Trinity.fasta.gene_trans_map の場所 --prep_reference マッピングの際にインデックスを作成することの明⽰ --thread_count 使⽤するスレッド数 実⾏例は次のスライドに書きます
コンティグ配列へのマッピング 使⽤するソフトウェア︓Trinity(2.13.2) こんな感じで実⾏します ubuntu@ubuntu-man$ perl /home/ubuntu/src/trinityrnaseq-v2.13.2/util/align_and_estimate_abundance.pl \ > --transcripts /homu/ubuntu/takai/trinity/trinity.Trinity.fasta \ > --seqType fq \ > --samples_file sample.tsv \ > --est_method RSEM \ > --aln_method bowtie2 \ > --output_dir /homu/ubuntu/takai/trinity/RSEM_result \ > --gene_trans_map /homu/ubuntu/takai/trinity/trinity.Trinity.fasta.gene_trans_map \ > --prep_reference \ > --thread_count 8 *環境によっては「は︖RSEM がないんやけど︖」みたいなエラーが吐かれます *その時は素直に RSEM をぐぐってインストールしてください(パスを通すのも忘れずに) 実⾏するとサンプルごとに以下のファイルができます ファイル名 ファイルの中⾝の⼤まかな説明(詳細はマニュアル参照) RSEM.isoforms.results 転写産物単位のマッピング結果 RSEM.genes.results 遺伝⼦単位のマッピング結果 転写産物単位、遺伝⼦単位どちらを解析に使うかはお任せします (転写産物単位で解析してる⼈が多い気がします)
発現量解析⽤のカウントデータ整形 使⽤するソフトウェア︓Trinity(2.13.2) 今の段階ではサンプルそれぞれでカウントデータが出⼒されているので Trinity のデータ整形⽤スクリプトを 使⽤して解析しやすいようにデータを整形します(1 つのファイルにまとめます) 使うスクリプトは Trinity をインストールしたディレクトリの util ディレクトリに⼊っています abundance_estimates_to_matrix.pl という名前のスクリプトです オプション 指定する内容の⼤まかな説明 --est_method 発現量予測に使⽤した⼿法 --gene_trans_map <output>.Trinity.fasta.gene_trans_map の場所 --quant_files それぞれのサンプルのカウントデータファイルをリストしたファイル --name_sample_by_basedir なし(ディレクトリ名をサンプル名として使⽤してね、という意味) とりあえず、サンプルファイルと同じ要領でカウントデータの場所を列挙したテキストファイルを作ります ubuntu@ubuntu-man$ nano RSEM_results_files.txt ubuntu@ubuntu-man$ cat RSEM_results_files.txt Con_1/RSEM.isoforms.results Con_2/RSEM.isoforms.results Con_5/RSEM.isoforms.results EPS_1/RSEM.isoforms.results EPS_4/RSEM.isoforms.results EPS_5/RSEM.isoforms.results
発現量解析⽤のカウントデータ整形 使⽤するソフトウェア︓Trinity(2.13.2) こんな感じで abundance_estimates_to_matrix.pl を実⾏します ubuntu@ubuntu-man$ perl /home/ubuntu/src/trinityrnaseq-v2.13.2/util/abundance_estimates_to_matrix.pl \ > --est_method RSEM \ > --gene_trans_map /home/ubuntu/takai/trinity/trinity.Trinity.fasta.gene_trans_map \ > --quant_files RSEM_results_files.txt \ > --name_sample_by_basedir 実⾏されるとたくさんファイルが出てきますが RSEM.isoform.counts.matrix(遺伝⼦単位で解析を⾏う場合は RSEM.gene.counts.matrix)を使います 中⾝はこんな感じです ubuntu@ubuntu-man$ head RSEM.isoform.counts.matrix Con_1 Con_2 Con_5 TRINITY_DN17980_c1_g1_i1 2.00 5.00 1.00 TRINITY_DN71091_c0_g1_i1 0.00 1.00 0.00 TRINITY_DN1157_c57_g1_i1 0.00 0.00 0.00 TRINITY_DN55094_c1_g1_i4 6.15 18.67 21.33 EPS_1 8.00 2.00 0.00 0.00 EPS_4 4.00 1.00 0.00 4.33 EPS_5 2.00 1.00 0.00 3.01 このファイルを使⽤して edgeR 等の解析ツールで発現量解析を⾏なっていきます(発現量解析編を参照) ただ Trinity ⾃体にも解析を⾏うスクリプトがあるので、それを使ってみても良いと思います(マニュアル参照) とりあえず⼤まかな流れはこれで終わりです 「どのコンティグがどの遺伝⼦なん︖」っていう話については後半でお話しします
リファレンスが あるっちゃある場合の De novo 解析 Trinity
リファレンスがあるっちゃある場合の de novo 解析とは 使⽤するソフトウェア︓マッピング⽤ソフトウェア、Trinity(v2.13.2) しっかりしたゲノム配列&遺伝⼦情報ファイルがある (いわゆるモデル⽣物である) うむ 前回までの解析⽅法で 頑張ってください (STAR からの edgeR) ない リファレンスに使えそうな 近縁種のゲノム配列もしくは ドラフトゲノム配列がある ある Trinity の ゲノムガイドアセンブリを 試してみましょう ドラフトゲノムとか 近縁種のゲノム情報とか Trinity で 他に頼りにできそうなものが 頑張りましょう ある場合の解析⽅法です そんなものはない
リファレンスがあるっちゃある場合の de novo 解析とは 使⽤するソフトウェア︓マッピング⽤ソフトウェア、Trinity(v2.13.2) ⼤雑把に⾔うと mRNA-Seq のリードを ドラフトゲノムとか近縁種へのゲノムに マッピングして そのマッピング情報を元に コンティグを組もうぜ というやり⽅です
リファレンス⽤配列へのマッピング 使⽤するソフトウェア︓HISAT2(2.2.1) まず、HISAT2 を使⽤してリファレンス(近縁種のゲノムとかドラフトゲノムとか)にマッピングをしていきます なぜ STAR ではないのかというと、、、 STAR はメモリ使⽤量が多いため、特にドラフトゲノムではメモリ容量が⾜りなくなることが多いからです 「私のつよつよ PC はメモリ不⾜とか起きないんですけど︖」という⽅は STAR でマッピングして OK です 以下、HISAT2 でのマッピング⽅法について説明します とりあえずインデックスを作成します ubuntu@ubuntu-man$ hisat2-build -p 8 GCA_002091915.1_Lexo_gen_Assembly01_genomic.fna funamushi で、全リードファイルをぶち込んでマッピングします(詳細は HISAT2 のマニュアルを参照してください) ubuntu@ubuntu-man$ hisat2 -x /home/ubuntu/references/funamushi/funamushi \ > -1 cleaned_C1_1.fq.gz,cleaned_C2_1.fq.gz,cleaned_C5_1.fq.gz,cleaned_E1_1.fq.gz,cleaned_E4_1.fq.gz,cleaned_E5_1.fq.gz > -2 cleaned_C1_2.fq.gz,cleaned_C2_2.fq.gz,cleaned_C5_2.fq.gz,cleaned_E1_2.fq.gz,cleaned_E4_2.fq.gz,cleaned_E5_2.fq.gz > -S /home/ubuntu/takai/10_HISAT2/aligned.sam \ > -p 8 \ \
リファレンス⽤配列へのマッピング 使⽤するソフトウェア︓HISAT2(2.2.1) HISAT2 を使⽤してリファレンス(近縁種のゲノムとかドラフトゲノムとか)にマッピングをしていきます なぜ STAR ではないのかというと、、、 STAR はメモリ使⽤量が多いため、特にドラフトゲノムではメモリ容量が⾜りなくなることが多いからです 「私のつよつよ 132875759 PC はメモリ不⾜とか起きないんですけど︖」という⽅は STAR でマッピングして OK です reads; of these: 132875759 (100.00%) were paired; of these: 以下、HISAT2 でのマッピング⽅法について説明します 93855550 (70.63%) aligned concordantly 0 times 38233098 (28.77%) aligned concordantly exactly 1 time とりあえずインデックスを作成します 787111 (0.59%) aligned concordantly >1 times ---- hisat2-build -p 8 GCA_002091915.1_Lexo_gen_Assembly01_genomic.fna funamushi ubuntu@ubuntu-man$ 93855550 pairs aligned concordantly 0 times; of these: 1778352 (1.89%) aligned discordantly 1 time で、マッピングします(詳細はマニュアルを参照してください) ---92077198 pairs 0 times concordantly or discordantly; ubuntu@ubuntu-man$ hisat2 -x aligned /home/ubuntu/references/funamushi/funamushi \ of these: 184154396 mates make up the pairs; of these: > -1 cleaned_C1_1.fq.gz,cleaned_C2_1.fq.gz,cleaned_C5_1.fq.gz,cleaned_E1_1.fq.gz,cleaned_E4_1.fq.gz,cleaned_E5_1.fq.gz 146922079 (79.78%) aligned 0 times > -2 cleaned_C1_2.fq.gz,cleaned_C2_2.fq.gz,cleaned_C5_2.fq.gz,cleaned_E1_2.fq.gz,cleaned_E4_2.fq.gz,cleaned_E5_2.fq.gz 35149894 (19.09%) aligned exactly 1 time > -S /home/ubuntu/lee/10_HISAT2/aligned.sam \ 2082423 (1.13%) aligned >1 times > -p 8 44.71% overall alignment rate 今回のマッピング結果はこんな感じでした(HISAT2 の出⼒結果) 44% のマッピング率ということで微妙な⾹りがぷんぷんしますが とりあえず気にせず進めてみます \ \
ゲノムガイドアセンブリ 使⽤するソフトウェア︓samtools(1.9)、Trinity(2.13.2) マッピングファイルを元に Trinity のゲノムガイドアセンブリを⾏います マニュアルによるとゲノムガイドアセンブリで使⽤するマッピングファイルは BAM 形式で かつソートされてないとだめらしいので samtools で対応します (STAR でマッピングしてソートされた BAM ファイルが既にある場合はこの作業は必要ありません) ubuntu@ubuntu-man$ samtools sort aligned.sam -@ 8 -O bam -o aligned_sort.bam BAM への変換&ソートが終わったら、Trinity でゲノムガイドアセンブリを実⾏します オプション 指定する内容の⼤まかな説明 --genome_guided_bam ドラフト or 近縁種ゲノムへのマッピング結果ファイル --genome_guided_max_intron 許容するイントロンの最⼤⻑ (とりあえずマニュアルと同じ 10000 塩基にしました) --max_memory 解析に使⽤するメモリの最⼤量(でかい⽅が良いです) --CPU 使⽤するスレッド数 --output 解析結果の出⼒先 ubuntu@ubuntu-man$ Trinity --genome_guided_bam aligned_sort.bam \ > --genome_guided_max_intron 10000 --max_memory 128G --CPU 8 \ > --output /home/ubuntu/takai/GG_trinity 時間かかるので気⻑に待ちます
ゲノムガイドアセンブリの実⾏結果 使⽤するソフトウェア︓Trinity(2.13.2) ⻑い⻑い待ち時間が終わると以下 2 つのファイルができます ファイル名 ファイルの中⾝の⼤まかな説明 Trinity-GG.fasta Trinity が構築したコンティグ配列 Trinity-GG.fasta.gene_trans_map 遺伝⼦と転写産物の対応表 このファイルを使⽤してマッピングしたいところではありますが ゲノムガイドアセンブリをした場合は、できあがったコンティグ配列を BUSCO で評価した⽅が良いです というよりも、絶対そうしましょう 詳細は次のスライド「どっちが良いん︖っていう話」に書きます BUSCO でコンティグの評価が「おっけーいい感じ」となったら マッピングからの発現量解析(本資料のスライド23-26)に進んでください
どっちが良いん︖ っていう話 BUSCO
コンティグ配列の評価 使⽤するソフトウェア︓BUSCO(5.2.2) リードだけを使ったコンティグ配列とマッピング情報も使ったコンティグ配列、どっちを使えばいいん︖ という話になるわけですが、簡単な⽅法として BUSCO を⽤いた評価⽅法があります BUSCO はアセンブリされたコンティグの中にコア遺伝⼦(core genes)がどれくらい含まれているかを 計算してくれるソフトウェアです コア遺伝⼦︓同⼀の分類群内で進化的に保存された遺伝⼦のことです anaconda を使⽤して BUSCO をインストールします(anaconda のインストール⽅法はぐぐってください) ubuntu@ubuntu-man$ conda create -n for_busco -c bioconda busco ubuntu@ubuntu-man$ conda activate for_busco こんな感じで使います(詳細は BUSCO のマニュアルを参照してください) ubuntu@ubuntu-man$ busco -m transcriptome -i trinity.Trinity.fasta -l arthropoda_odb10 -o BUSCO -c 8 今回の結果はこんな感じでした(データ︓フナムシ、BUSCO データベース︓arthropoda_odb10) アセンブリの⽅法 コア遺伝⼦の割合(コンティグ中のコア遺伝⼦数/コア遺伝⼦数) リアルガチ (mRNA-Seq リードのみ) 97.7% (990/1013) ゲノムガイドアセンブリ (mRNA-Seq + マッピング情報) 51.3% (520/1013) リアルガチのコンティグを 使った⽅が良さそうですね
コンティグ配列の評価 使⽤するソフトウェア︓BUSCO(5.2.2) リードだけを使ったコンティグ配列とマッピング情報も使ったコンティグ配列、どっちを使えばいいん︖ という話になるわけですが、簡単な⽅法として BUSCO を⽤いた評価⽅法があります BUSCO はアセンブリされたコンティグの中にコア遺伝⼦(core genes)がどれくらい含まれているかを 計算してくれるソフトウェアです コア遺伝⼦︓同⼀の分類群内で進化的に保存された遺伝⼦のことです anaconda を使⽤して BUSCO をインストールします(anaconda のインストール⽅法はぐぐってください) ubuntu@ubuntu-man$ conda create -n for_busco -c bioconda busco ubuntu@ubuntu-man$ conda activate for_busco こんな感じで使います(詳細は BUSCO のマニュアルを参照してください) ubuntu@ubuntu-man$ busco -m transcriptome -i trinity.Trinity.fasta -l arthropoda_odb10 -o BUSCO -c 8 今回の結果はこんな感じでした(データ︓フナムシ、BUSCO データベース︓arthropoda_odb10) アセンブリの⽅法 コア遺伝⼦の割合(コンティグ中のコア遺伝⼦数/コア遺伝⼦数) リアルガチ (mRNA-Seq リードのみ) 97.7% (990/1013) ゲノムガイドアセンブリ (mRNA-Seq + マッピング情報) 51.3% (520/1013) フナムシのドラフトゲノムは スキャフォールドレベルでした なのでこの結果も納得です そもそもマッピングの時点で、、、
コンティグ配列の評価 使⽤するソフトウェア︓BUSCO(5.2.2) 登 リードだけを使ったコンティグ配列とマッピング情報も使ったコンティグ配列、どっちを使えばいいん︖ 録 という話になるわけですが、簡単な⽅法として BUSCO を⽤いた評価⽅法があります さ BUSCO はアセンブリされたコンティグの中にコア遺伝⼦(core genes)がどれくらい含まれているかを れ 計算してくれるソフトウェアです い コア遺伝⼦︓同⼀の分類群内で進化的に保存された遺伝⼦のことです て ろ anaconda を使⽤して BUSCO をインストールします(anaconda のインストール⽅法はぐぐってください) い い ubuntu@ubuntu-man$ conda create -n for_busco -c bioconda busco る ubuntu@ubuntu-man$ conda activate for_busco ろ ゲ あ こんな感じで使います(詳細は BUSCO のマニュアルを参照してください) ノ る-c 8 ubuntu@ubuntu-man$ busco -m transcriptome -i trinity.Trinity.fasta -l arthropoda_odb10 -o BUSCO ム ん に 今回の結果はこんな感じでした(データ︓フナムシ、BUSCO データベース︓arthropoda_odb10) だ た な も アセンブリの⽅法 コア遺伝⼦の割合(コンティグ中のコア遺伝⼦数/コア遺伝⼦数) か あ リアルガチ 97.7% (990/1013) フナムシのドラフトゲノムは (mRNA-Seq リードのみ) スキャフォールドレベルでした い ゲノムガイドアセンブリ (mRNA-Seq + マッピング情報) 51.3% (520/1013) なのでこの結果も納得です そもそもマッピングの時点で、、、
遺伝⼦情報の付与 blast
blast のインストール blast(v2.12.0) Trinity はコンティグ配列を構築してくれて、カウントデータまで出⼒してくれますが それぞれのコンティグがどのような遺伝⼦(転写産物)なのかはわかりません なので BLAST 検索を使⽤してそれぞれのコンティグに遺伝⼦情報の付与を⾏います まずは blast をインストールします ubuntu@ubuntu-man$ wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast2.12.0+-x64-linux.tar.gz ubuntu@ubuntu-man$ tar xvf ncbi-blast-2.12.0+-x64-linux.tar.gz ubuntu@ubuntu-man$ ln -s /home/ubuntu/src/ncbi-blast-2.12.0+/bin/blastx /home/ubuntu/bin/ ubuntu@ubuntu-man$ ln -s /home/ubuntu/src/ncbi-blast-2.12.0+/bin/makeblastdb /home/ubuntu/bin/ BLAST 検索⽤のデータベースをダウンロードします 凄まじく容量が⼤きいので専⽤のディレクトリを作った⽅が良いです ubuntu@ubuntu-man$ pwd /home/ubuntu/blastdb ubuntu@ubuntu-man$ /home/ubuntu/src/ncbi-blast-2.12.0+/bin/update_blastdb.pl nr --decompress これで BLAST 検索を⾏う環境が整いました ネット環境によっては update_blastdb.pl がうまく動かないことがあります その時は https://ftp.ncbi.nlm.nih.gov/blast/db/v5 から nr.xx.tar.gz を wget 等で直接ダウンロード&解凍してください
⽣物種を指定した BLAST 検索 blast(v2.12.0) 今から Trinity のコンティグ配列を BLAST 検索していくわけですが できれば⾃分のサンプルと近い⽣物種で BLAST したいなあ、ということが⼤半だと思います 実際、今回はフナムシのサンプルなので節⾜動物を対象に BLAST したいなあと思っています BLAST の中の get_species_taxids.sh を使うと 対象にする⽣物種の TaxID リストが簡単に取得できます まず -n オプションで調べたい⽣物種、分類群の TaxID を調べます ubuntu@ubuntu-man$ ./get_species_taxids.sh -n Arthropoda Cannot find Entrez EDirect esearch tool, please see installation in https://www.ncbi.nlm.nih.gov/books/NBK179288/ と思ったらエラーになったので、⾔われた通り Entrez Edirect esearch を⼊れてもう⼀度 ubuntu@ubuntu-man$ ./get_species_taxids.sh -n Arthropoda Taxid : 6656 節⾜動物(Arthropoda)の rank : phylum TaxID は 6656 でした division : arthropods scientific name : Arthropoda common name : arthropods 1 matche(s) found. このスクリプトはとても優秀で そこそこ適当に⽣物の名前を⼊れて ばっちり当てはまる検索結果がなくても 「お前が探しとるのはこれか︖」みたいに これっぽいぞリストを出してくれます
⽣物種を指定した BLAST 検索 blast(v2.12.0) 検索した TaxID を使って BLAST 対象の⽣物種の TaxID リストを作ります ubuntu@ubuntu-man$ ./get_species_taxids.sh -t 6656 > /home/ubuntu/blastdb/taxids_Arthropoda_6656.txt *解析⽤ PC の環境によってはセキュリティでこれが実⾏できないことがあります *その場合は同じコマンドを他の PC から実⾏してみてください(もちろんインストールからです) TaxID のリストの中⾝はこんな感じになっています ubuntu@ubuntu-man$ head taxids_Arthropoda_6656.txt 6656 6657 6658 この TaxID リストを使って BLAST 検索を⾏なっていきます
⽣物種を指定した BLAST 検索 blast(v2.12.0) BLAST 検索(blastx)を実⾏します オプション 指定する内容の⼤まかな説明 -query クエリ配列(Trinity のコンティグ配列) -db 使⽤するデータベース(nr) -taxidlist 検索する⽣物種の TaxID リストファイル -outfmt 出⼒形式(マニュアル参照) -evalue E-value の閾値(0.001 がよく使われてる気がします) -out BLAST 検索の結果を出⼒するファイルの名前 -num_threads BLAST 検索に使⽤するスレッド数 ubuntu@ubuntu-man$ blastx -query /home/ubuntu/takai/trinity/trinity.Trinity.fasta \ > -db /home/ubuntu/blastdb/nr -taxidlist /home/ubuntu/blastdb/taxids_Arthropoda_6656.txt \ > -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore stitle" > -evalue 1e-03 \ > -out /home/ubuntu/takai/blast/blastx_Arthropoda.tsv -num_threads 8
⽣物種を指定した BLAST 検索 blast(v2.12.0) BLAST 検索(blastx)を実⾏します これこのまま実⾏したら 1 つのコンティグに対して オプション 指定する内容の⼤まかな説明 複数の検索(ヒット)結果が返ってくるんじゃない︖ Trinity のコンティグ配列ファイル -query -db 使⽤するデータベース(nr) -taxidlist 検索する⽣物種の TaxID リストファイル -outfmt 出⼒形式(マニュアル参照) そうなんです E-value の閾値(0.001 がよく使われてる気がします) -max_target_seqs っていうオプションで BLAST 検索の結果を出⼒するファイルの名前 -out 結果が 1 つだけ返ってくるようにしたらいいんじゃない︖ BLAST 検索に使⽤するスレッド数 -num_threads -evalue そんな⾵に考えている時期が俺にもありました ubuntu@ubuntu-man$ blastx -query /home/ubuntu/takai/trinity/trinity.Trinity.fasta \ > -db /home/ubuntu/blastdb/nr -taxidlist /home/ubuntu/blastdb/taxids_Arthropoda_6656.txt \ > -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore stitle" > -evalue 1e-03 \ > -out /home/ubuntu/takai/blast/blastx_Arthropoda.tsv -num_threads 8
⽣物種を指定した BLAST 検索 blast(v2.12.0) BLAST 検索(blastx)を実⾏します これ本当に⼤きな落とし⽳なのですが オプション 指定する内容の⼤まかな説明 Trinity のコンティグ配列ファイル -max_target_seqs オプションは 使⽤するデータベース(nr) -db 必ずしもベストヒットの結果を返すわけではない 検索する⽣物種の TaxID リストファイル -taxidlist -query -outfmt らしいです 出⼒形式(マニュアル参照) -evalue E-value の閾値(0.001 がよく使われてる気がします) BLAST 検索の結果を出⼒するファイルの名前 ただ結果はヒット順(E-value)にソートされるので BLAST 検索に使⽤するスレッド数 -num_threads とりあえず複数ヒットは気にせずに検索して BLAST blastx 検索が終わったらそのファイルからベストヒットを ubuntu@ubuntu-man$ -query /home/ubuntu/takai/trinity/trinity.Trinity.fasta \ > -db /home/ubuntu/blastdb/nr -taxidlist /home/ubuntu/blastdb/taxids_Arthropoda_6656.txt \ 抜き出してくる、というのが安全です > -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore -out stitle" > -evalue 1e-03 \ > -out /home/ubuntu/takai/blast/blastx_Arthropoda.tsv -num_threads 8
⽣物種を指定した BLAST 検索 とりあえず実⾏ めちゃくちゃ 時間かかります blast(v2.12.0) BLAST 検索(blastx)を実⾏します オプション 指定する内容の⼤まかな説明 -query Trinity のコンティグ配列ファイル -db 使⽤するデータベース(nr) -taxidlist 検索する⽣物種の TaxID リストファイル -outfmt 出⼒形式(マニュアル参照) -evalue E-value の閾値(0.001 がよく使われてる気がします) -out BLAST 検索の結果を出⼒するファイルの名前 -num_threads BLAST 検索に使⽤するスレッド数 ubuntu@ubuntu-man$ blastx -query /home/ubuntu/takai/trinity/trinity.Trinity.fasta \ > -db /home/ubuntu/blastdb/nr -taxidlist /home/ubuntu/blastdb/taxids_Arthropoda_6656.txt \ > -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore stitle" > -evalue 1e-03 \ > -out /home/ubuntu/takai/blast/blastx_Arthropoda.tsv -num_threads 8 ⼀週間、⼆週間、⼀ヶ⽉⻑い⻑い待ち時間です
⽣物種を指定した BLAST 検索 blast(v2.12.0) 「時間もあるし、特に急ぐ必要もないわ」という場合はただ待つだけでも OK です 「そういえば BLAST やってたな」みたいにふと BLAST を思い出す頃に終わると思います ただ「は︖報告書間に合わないんですけど︖」みたいな場合もあると思います そういう場合はとりあえず発現量解析をやって、発現量が変わった遺伝⼦(のコンティグ)だけを BLAST するっていうのも⼀つの⼿段だと思います (seqkit というツールを使うと⽐較的簡単にこの作業を実⾏できます) さらにさらにどこにでもハプニングはあるわけで BLAST 検索中に解析⽤ PC が落ちたりする(間違って消されたりする)こともあるわけです なので僕は BLAST 検索を⾏うコンティグ配列ファイルを seqkit を使って 1000 コンティグずつ くらいに分割して、そのファイルに対して BLAST 検索を実⾏させています こうしておけば途中で万が⼀、BLAST が⽌まったときも、どこまでは検索できてて どこからやり直せば良いかが明確になります(つまり発狂する時間が短くなります)
BLAST 検索結果からベストヒットの抽出
cat, awk, grep
先ほど「ベストヒットの抽出︖簡単だよ︖」みたいに⾔いましたが
この⽅法は僕が考えたわけではなく以下のサイトで⾒つけました
blast+ 使い⽅ best hitの算出 awkとoutfmt7(バイオインフォ道場 [bioinfo-Dojo])
https://bioinfo-dojo.net/2016/03/25/blast_besthit_outfmt7/
そうか、こうすれば良いのかと脱帽したものです
詳細な説明は本家のサイトがとてもわかりやすいのでそちらを参照していただけたらと思います
⼤まかには
・cat <BLAST 検索の結果ファイル> BLAST 検索結果ファイルを表⽰
・awk '/hits found/{getline;print}ʼ 表⽰した内容から hits found という⽂字列を含む⾏の
次の⾏を表⽰
・grep -v "#" 出⼒された⾏から # がついている⾏以外を抽出
この 3 つのコマンドをパイプで繋げることによってベストヒットを抽出しています(脱帽)
ubuntu@ubuntu-man$ cat blastx_Arthropoda.tsv | awk '/hits found/{getline;print}' | grep -v "#" >
blastx_Arthropoda_BH.txt
BLAST 検索が終わったら blast(v2.12.0) ⻑い⻑い待ち時間を経て BLAST 検索が終わり、無事にベストヒットを抽出したら あとは BLAST 検索の結果と発現量解析の結果と照らし合わせて頑張って考察していくのみです 今回はコンティグ配列の遺伝⼦情報を得るために BLAST 検索のみを使⽤しましたが Trinotate というツールを使⽤するともっと幅広く実験結果を観察することができるようになります Trinotate の使い⽅は Trinity のマニュアルでも紹介されているので気になる⽅は読んでみてください De novo での解析は待ち時間に 焦らされ続けることになると思います 僕たちが焦っても解析速度は上がらないので 他の実験を進めつつ気⻑に頑張りましょう