ex401 - nibb-gitc/gitc2025mar-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リファレンスにマッピングさせ、発現量を推定する。

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。パスが通っていなかったり、インストールがうまくいってなければエラーが出る。

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 サイズ。

(発展) 通常はtranscriptのインデクシングで十分であるが、より正確性を高めるには、decoy-awareなインデックス作りが推奨されている。この練習問題は基礎を習得することを目的としているので深入りしない。decoy-aware indexing の詳細と実例については、case1を参照されたい。

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

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

...
[2022-02-26 10:43:48.855] [jointLog] [info] Mapping rate = 36.5575%
...

のようなレポートが表示されて終了する。今回は練習のために、インプットのリード数を減らしているなどの理由で、いくつかエラーメッセージも表示されるが、ここでは気にしない。

Inspect results

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

quant.sf が結果ファイル。中身を確認しよう。(lessコマンドなど)

(example)

Name                    Length  EffectiveLength TPM             NumReads
lcl|ENSMUST00000074761  381     131.000         132.530453      4.000
lcl|ENSMUST00000136312  2205    1955.000        87.535907       39.429
...

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

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

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

以下、R 環境にて。(リモート、ローカル問わない。ログインしているリモート環境でRを起動するには、コマンドラインで、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] 73115
> sum(x$TPM)
[1] 1e+06

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

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

2025-3-4

  • RCCS環境で動作確認。（salmon v1.10.0)

2024-7-30

  • salmon v1.10.0 で動作確認。

2024-2-26

  • decoy-aware indexingについて補足を追加。

2022-8-27

  • salmon v1.8.0で動作確認。v1.7.0の時の出力とほんの少し異なるが、本文はv1.7.0のまま。

2022-2-28

  • Minor updates for GITC 2022 Mar. (salmon v1.7.0; R 3.6.3).

2021-9-11

  • Minor updates for GITC 2021 Sep.

2021-3-9

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

2020-05-24

  • Initial release. Salmon version: 1.2.1