チョウが食草を見分けるしくみを探る | 4. リード数カウント

2020/05/04 公開

リード数をカウントする

Trinity でアセンブルしたcontigを利用して、発現量推定の第一歩であるリード数のカウントを行う。

Trinity のツールについて

Trinity に付属するツール群は親切で、起動すると使い方や設定できるオブションが表示される。
例:
:> $TRINITY_HOME/util/align_and_estimate_abundance.pl
とだけ入力してリターンを押す。
最後の方に実行例があるので、とても参考になる。

リード数カウントの手順

Trinityのアウトプット・フォルダーに移動
フォルダーを作る

:> mkdir ./count

作ったフォルダーに移動

:> $ cd ./count

解析用のスペース区切りファイルを作る

:> ls /sample/fastq/trimmed/*_P.fastq.gz > BRR_file.txt

「/sample/fastq/trimmed」のような、ファイルの置き場所を示すことを「ファイルのパス」という。ファイルのパスは、各自がファイルを保存した場所を指定する。

BRR_file.txt を開いて編集（半角スペースで区切る。スペースが全角にならないように注意）

sample_A sample_A_rep1 A_rep1_R1.fastq.gz A_rep1_R2.fastq.gz
sample_A sample_A_rep2 A_rep2_R1.fastq.gz A_rep2_R2.fastq.gz
sample_B sample_B_rep1 B_rep1_R1.fastq.gz B_rep1_R2.fastq.gz
sample_B sample_B_rep2 B_rep2_R1.fastq.gz B_rep2_R2.fastq.gz

上記の例では、sample_Aについて二回の繰り返し（sample_A_rep1, sample_A_rep2）、sample_Bについても二回の繰り返し（sample_B_rep1, sample_B_rep2）を行った実験であることを示している。各サンプルは、ペアエンド（DNA断片の両末端から塩基配列を読んだ）である（R1とR2）ことを示している。

解析用フォルダーを作る（非アライメント: 高速）

:> mkdir ./salmon

解析用フォルダーに移動

:> cd ./salmon

ファイルをコピーする

:> cp ../../Trinity.fasta ./

スクリプトファイルを作る

:> touch ./est_salmon.sh

スクリプトファイルを開いて編集

perl $TRINITY_HOME/util/align_and_estimate_abundance.pl --thread_count 16 --transcripts ./Trinity.fasta --trinity_mode --seqType fq --prep_reference --samples_file ../BRR_file.txt --est_method salmon

スクリプトを実行

:> /bin/sh ./est_salmon.sh

計算が終わると、サンプル毎にフォルダーが作られていて、その中に「quant.sf」というファイルができている。これがsalmonのカウントファイル本体。
テキストエディター等で開いて結果を確認できる。
Name, Length, EffectiveLength, TPM, NumReads というカラムがあるが、TPMが相対的な発現量（その組織の中でどれくらい発現しているか）を示す。
TPMが高い順に並べ直して、contig名とTPMだけ取り出したい場合は下記のコマンドを実行する。
:> sort -rn -k 4 ./quant.sf | awk '{print $1,$4}' > sort_TPM.txt

sort コマンドの使い方
-r: 逆順、つまり大きい順に並べる
-k: 並べ替えのキーを指定。-k 2 だと2列目をキーとして使うという意味。
-n キーを数字として扱うという意味。

発現量が多い遺伝子を見つけたいという研究目的であれば、ここまでで解析完了。TPMの値がいくつ以上、というような基準で必要なcontigを抜き出して、研究目的に合わせた作業を行う。
組織間や何らかの操作を行なった場合など、異なる条件で標的の遺伝子の発現量が変動しているか確認したい場合は、本マニュアルのこの後の解析を続ける。

解析用フォルダーを作る（アライメント: 高精度）

:> mkdir ./RSEM_bowtie2

解析用フォルダーに移動

:> cd ./RSEM_bowtie2

ファイルをコピーする

:> cp ../../Trinity.fasta ./

スクリプトファイルを作る

:> touch ./est_RSEM.sh

スクリプトファイルを開いて編集

perl $TRINITY_HOME/util/align_and_estimate_abundance.pl --thread_count 16 --transcripts ./Trinity.fasta --trinity_mode --seqType fq --prep_reference --samples_file ../BRR_file.txt --est_method RSEM --aln_method bowtie2

スクリプトを実行

:> /bin/sh ./est_RSEM.sh

計算が終わると、サンプル毎にフォルダーが作られていて、その中に「RSEM.genes.results」と「RSEM.isoforms.results」ができている。どちらもカウントファイルではあるが、以後の解析には「RSEM.isoforms.results」を使う。
テキストエディター等で開いて結果を確認できる。
transcript_id, gene_id, length, effective_length, expected_count, TPM, FPKM, IsoPct というカラムがあるが、TPMが相対的な発現量（その組織の中でどれくらい発現しているか）を示す。
TPMが高い順に並べ直して、contig名とTPMだけ取り出したい場合は下記のコマンドを実行する。
:> sort -rn -k 6 ./RSEM.isoforms.results | awk '{print $1,$6}' > sort_TPM.txt

sort コマンドの使い方
-r: 逆順、つまり大きい順に並べる。
-k: 並べ替えのキーを指定。-k 2 だと2列目をキーとして使うという意味。
-n キーを数字として扱うという意味。

どのソフトを使うか

salmon が高速で、十分な精度になっているので、salmon だけで良いのではないかと思うが、現状では最も信頼されて多く使われているのが RSEM + bowtie2 の組み合わせ。
この他、kallisto を使った方法と RSEM + bowtie を使った方法も可能であるが、行う必要はないだろうと思う。

個人的には、念のため salmon と RSEM + bowtie2 の両方を試しているが、salmon だけで充分である場合がほとんど。

スクリプトファイル（.sh という拡張子をつけたファイル）を作成する意味

スクリプトファイルを作成して作業したフォルダー内に残しておくことで、作業ログにもなる。そのフォルダー内に残された解析結果が、どのようなコマンドによって得られたものなのか確実な記録になるので、完璧な再現性を保証する研究ノート代わりにもなる。紙のノートに手書きする際に発生する転記ミスを防ぐ意味もあるので、解析研究での研究ノートの残し方としてお薦めする。
スクリプトは「#」から改行コードまでをメモとして扱い、コマンドには影響しないという仕組みがあるので、上記のスクリプ内に実行した日付や研究の条件などを書き込んでおくと、研究ノートとしての情報量が上がってより有益である。
例
# 2020/05/05
# ナミアゲハのRNA-seq

Q&A

皆様からいただいた質問と回答など、こちらへ掲載します。

質問やご意見はこちらへ

質問やご意見はこちらへ（投稿フォームに遷移します）

※ボタンのリンク先は、皆で語り合う生命誌研究館の掲示板「みんなの広場」です。ご投稿いただいたご質問やご意見、館員からのお返事は公開されますのでご了承ください。

昆虫食性進化研究室