ex14 - nibb-gitc/gitc2020jun-rnaseq GitHub Wiki

ex14

この練習問題では、Transcript-based RNA-seq pipelineとして、salmonソフトを使った alignment-free のread count estimation のパイプラインを学ぶ。

マウス Mus musculus のRNA-seqを行った。ライブラリは1種類のみで、single end (片側 Read1のみ) 75bpシークエンスを行った。これらのリードをマウスmRNAリファレンスにマッピングさせ、発現量を推定する。

Data

データファイルは、~/data/SS 以下。

Input reads

  • IlluminaReads1.fq : FASTQ format. 200,000 reads included.

Reference

  • minimouse_mRNA.fa : FASTA format. 5,000 transcripts included. (本来マウスの遺伝子数は2万個以上だが、練習用に一部を抽出してある）

Setup

Setup environment

ex14 ディレクトリをつくり、以下の解析はその下で作業しよう。

$ mkdir ex14
$ cd ex14

使うファイルをdataディレクトリからコピーしておこう。

$ cp ~/data/SS/minimouse_mRNA.fa ./
$ cp ~/data/SS/IlluminaReads1.fq ./

salmonが動くかどうか確認

$ salmon -h

helpメッセージが表示されればOK。パスが通っていなかったり、インストールがうまくいってなければエラーが出る。

Build index of reference

salmon index -t minimouse_mRNA.fa -i minimouse_mRNA.fa.salmon.idx -k 31

optionの意味はヘルプで調べられる。

salmon index -h

重要なoptionの説明

  • -t: インデックスを作成するリファレンス配列のファイル名を指定する。
  • -i: 作成するインデックスの名前。好きな名前で良い。
  • -k: kmer サイズ。

Quantification with salmon

salmonでabundance estimation

salmon quant -i minimouse_mRNA.fa.salmon.idx -l A -o salmon_out -p 2 -r IlluminaReads1.fq

optionの意味はヘルプで調べられる。

salmon quant --help-reads

コマンドを実行するとしばらくして、

...
[2020-05-24 13:49:00.750] [jointLog] [info] Mapping rate = 36.638%
...

のようなレポートが表示されて終了する。マッピング率など有用な情報なので、テキストファイルにコピー＆ペーストして保存しておくと良い。（今回は練習のために、インプットのリード数を減らしているなどの理由で、いくつかエラーメッセージも表示されるが、ここでは気にしない）

Inspect results

計算が終わったら、出力ディレクトリ（今回の場合、-oで指定した　salmon_out）にどのようなファイルが生成されたか確認する。 (''ls -l''など)

quant.sf が結果ファイル。中身を確認しよう。

(example)

Name                    Length  EffectiveLength TPM             NumReads
lcl|ENSMUST00000074761  381     132.000         135.948330      4.000
lcl|ENSMUST00000136312  2205    1956.000        90.483290       39.451
...

format の説明は、開発者のHPを参照すること。

https://salmon.readthedocs.io/en/latest/file_formats.html#fileformats

quant.sf　の内容をもう少し詳しく調べてみよう。

以下、R 環境にて。

> x <- read.delim("quant.sf")
> head(x)
                    Name Length EffectiveLength        TPM NumReads
1 lcl|ENSMUST00000074761    381             132  135.94833    4.000
2 lcl|ENSMUST00000136312   2205            1956   90.48329   39.451
3 lcl|ENSMUST00000004316   1671            1422    0.00000    0.000
4 lcl|ENSMUST00000105465   1665            1416 5830.57819 1840.332
5 lcl|ENSMUST00000165878   1656            1407    0.00000    0.000
6 lcl|ENSMUST00000177779   1674            1425    0.00000    0.000

> sum(x$NumReads)
[1] 73276
> sum(x$TPM)
[1] 1e+06

total のマップされたreads数が73276だとわかる。インプットのリード数が200000だったので、73276 / 200000 = 0.366 つまりmapping rate = 36.6% 上記のログファイルの"[info] Mapping rate = 36.638%"と整合性が取れている。

TPM の定義から、全ての合計が100万になるのは当然。（TPMの定義を思い出してみよう。）

Links

References

  • Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14: 417–419. doi:10.1038/nmeth.4197

Revision history

2020-05-24

  • Initial release. Salmon version: 1.2.1