ex401 - nibb-gitc/gitc2021mar-rnaseq GitHub Wiki

ex401

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

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

（受講生の皆さんは、bias5サーバにログインして以下の演習を行います。）

Data

データファイルは、~/gitc/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

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

$ mkdir ex401
$ cd ex401

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

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

salmonが動くかどうか確認

$ salmon -h

helpメッセージが表示されればOK。パスが通っていなかったり、インストールがうまくいってなければエラーが出る。今回の環境ではsalmon v1.4.0が起動するはず。

Build index of reference

リファレンスとなるトランスクリプトームの配列データをインデクシング（索引作成）する。salmon indexサブコマンドを使う。

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

重要・よく使うオプションは

  • -i : salmon index
  • -o : name of output directory
  • -l : library type. A: 自動推定。通常はAで良い。詳細はマニュアル
  • -r : readのファイル名（single-endの場合）
  • -1, -2: readのファイル名（paired-endの場合）
  • -p : 使うCPU coreの数。最近のコンピュータは複数のCPU coreを持っているので多く使えば計算は早く終わる。ただし今回は多数の受講生が同じサーバを利用するため２以上に増やさないこと。

詳細はマニュアルを参照：https://salmon.readthedocs.io/en/latest/salmon.html

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

...
[2021-03-09 16:45:02.747] [jointLog] [info] Mapping rate = 36.559%
...

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

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             131  132.52562    4.000
2 lcl|ENSMUST00000136312   2205            1955   87.53271   39.429
3 lcl|ENSMUST00000004316   1671            1421    0.00000    0.000
4 lcl|ENSMUST00000105465   1665            1415 5644.75627 1840.338
5 lcl|ENSMUST00000165878   1656            1406    0.00000    0.000
6 lcl|ENSMUST00000177779   1674            1424    0.00000    0.000

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

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

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

2021-3-9

  • Minor updates for GITC 2021 Mar. Salmon version: 1.4.0

2020-05-24

  • Initial release. Salmon version: 1.2.1