前回の記事↓の続きです。
STAR-RSEMによる発現量推定 その１ - Palmsonntagmorgen
STAR-RSEMによる発現量推定 その2 - Palmsonntagmorgen
前回のコマンドを最後まで完了すると、starディレクトリの中にstar/Myers_HUVEC_cell_2x75_200_1.<genes|isoforms>.results のようなファイルがサンプルごとに生成されているのがわかります。
*.genes.resultsと*.isoforms.resultsはそれぞれ遺伝子単位、transcript単位の発現量（count, TPM, FPKM）が出力されます。
transcript単位の出力は基本的には精度が十分でないことが多く、ノイズがまぎれこみやすいので、特に利用目的がないのであれば遺伝子単位のものを利用しましょう。
RSEMの出力
では、ファイルの中をのぞいてみましょう。
$ head star/Myers_HUVEC_cell_2x75_200_1.genes.results
gene_id transcript_id(s) length effective_length expected_count TPM FPKM
ENSG00000000003 ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2146.69 1926.19 1086.00 40.20 34.28
ENSG00000000005 ENST00000373031,ENST00000485971 1046.60 826.17 14.00 1.21 1.03
ENSG00000000419 ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1003.23 782.75 355.00 32.34 27.57
ENSG00000000457 ENST00000367770,ENST00000367771,ENST00000367772,ENST00000423670,ENST00000470238 3672.40 3451.91 112.23 2.32 1.98
ENSG00000000460 ENST00000286031,ENST00000359326,ENST00000413811,ENST00000459772,ENST00000466580,ENST00000472795,ENST00000481744,ENST00000496973,ENST00000498289 2167.45 1946.97 501.77 18.38 15.67
ENSG00000000938 ENST00000374003,ENST00000374004,ENST00000374005,ENST00000399173,ENST00000457296,ENST00000468038,ENST00000475472 2021.00 1800.50 9.00 0.36 0.30
ENSG00000000971 ENST00000359637,ENST00000367429,ENST00000466229,ENST00000470918,ENST00000496761,ENST00000630130 2587.83 2367.38 0.00 0.00 0.00
ENSG00000001036 ENST00000002165,ENST00000367585,ENST00000451668 2285.18 2064.69 1079.67 37.29 31.79
ENSG00000001084 ENST00000229416,ENST00000504353,ENST00000504525,ENST00000505197,ENST00000505294,ENST00000509541,ENST00000510837,ENST00000513939,ENST00000514004,ENST00000514373,ENST00000514933,ENST00000515580,ENST00000616923,ENST00000643939,ENST00000650454 2327.40 2106.94 773.00 26.16 22.30
Ensemblからダウンロードしたgtfファイルを用いているので、gene_idとtranscript_idは当然Ensemblのものになっています。
effective length はpoly(A) tailを除いた転写物の長さだそうです。
expected_count は(isoformごとに振り分けられた)マップリード数、FPKMはリード数を全マップリード数と遺伝子長で正規化したスコアで、TPMはFPKMを更にFPKMの総和で正規化した値になります。
発現量の評価にはTPM、edgeRの入力にはcountを使おう
遺伝子発現量を直接見たい場合、例えば box plot や scatter plot を描く場合には TPMを使っておけば基本的に間違いありません。
一方、edgeRやDESeq2のような発現変動比較ツールにはTPMではなくcountデータを利用することが推奨されます（マニュアル参照）。何故なら、edgeRやDESeq2の中で用いている正規分布や負の二項分布はraw countに対してモデル化するものであり、遺伝子長などで正規化してしまうとこの分布を満たさなくなってしまうおそれがあるからです。全ての正規化はedgeRやDESeq2の中で行われますので、前もって正規化しておく必要はありません。更に言えば、サンプル間比較は遺伝子ごとに行われますので、基本的に遺伝子長の正規化は必要ないのです。
なお、DESeq2は整数値しか扱えませんので、DESeq2を利用する時にはcountデータの少数を切り捨てにしておきましょう。
発現量データをマージ (count)
それでは、４つのサンプルの発現量データをマージしてみましょう。RSEMの中にある rsem-generate-data-matrixを使います。
ここでは遺伝子ごとのファイルを用いています。
$ RSEM/rsem-generate-data-matrix \
star/Myers_H1-hESC_cell_2x75_200_1.genes.results \
star/Myers_H1-hESC_cell_2x75_200_2.genes.results \
star/Myers_HUVEC_cell_2x75_200_1.genes.results \
star/Myers_HUVEC_cell_2x75_200_2.genes.results \
> GeneExpressionMatrix.tsv
得られた結果を見てみましょう。
$ head GeneExpressionMatrix.tsv
"star/Myers_H1-hESC_cell_2x75_200_1.genes.results" "star/Myers_H1-hESC_cell_2x75_200_2.genes.results" "star/Myers_HUVEC_cell_2x75_200_1.genes.results" "star/Myers_HUVEC_cell_2x75_200_2.genes.results"
"ENSG00000000003" 1086.00 2327.00 1086.00 1891.00
"ENSG00000000005" 14.00 34.00 14.00 0.00
"ENSG00000000419" 355.00 755.00 355.00 875.00
"ENSG00000000457" 112.23 212.28 112.23 145.02
"ENSG00000000460" 501.77 878.72 501.77 201.98
"ENSG00000000938" 9.00 15.00 9.00 1.00
"ENSG00000000971" 0.00 0.00 0.00 7301.00
"ENSG00000001036" 1079.67 2036.77 1079.67 1889.43
"ENSG00000001084" 773.00 1263.00 773.00 608.00
行が遺伝子、列がサンプルの行列データになっていることがわかります。タブ区切りのテキストファイルなので、ExcelやR, Pythonなど好きなもので分析することが可能です。
発現量データをマージ (TPM)
上記で出力されているデータはカウントデータです。rsem-generate-data-matrixはカウントデータにしか対応していません（手抜き？）。
TPM や RPKM のマトリックスを出力したい場合は、僕が作成した rsem-generate-data-matrix-modified を使ってください。
rsem-generate-data-matrix-modifiedは、↓の記事で紹介しているscript_rnakato の中にあります。
GitHubからプログラムをダウンロード・インストール - Palmsonntagmorgen
ではやってみましょう。
$ script_rnakato/rsem-generate-data-matrix-modified \
TPM \
star/Myers_H1-hESC_cell_2x75_200_1.genes.results \
star/Myers_H1-hESC_cell_2x75_200_2.genes.results \
star/Myers_HUVEC_cell_2x75_200_1.genes.results \
star/Myers_HUVEC_cell_2x75_200_2.genes.results \
> GeneExpressionMatrix.TPM.tsv
$ head GeneExpressionMatrix.TPM.tsv
Myers_H1-hESC_cell_2x75_200_1.genes.results Myers_H1-hESC_cell_2x75_200_2.genes.results Myers_HUVEC_cell_2x75_200_1.genes.results Myers_HUVEC_cell_2x75_200_2.genes.results
ENSG00000000003 40.20 41.67 40.20 39.30
ENSG00000000005 1.21 1.91 1.21 0.00
ENSG00000000419 32.34 34.32 32.34 39.18
ENSG00000000457 2.32 2.47 2.32 1.88
ENSG00000000460 18.38 15.50 18.38 4.85
ENSG00000000938 0.36 0.30 0.36 0.02
ENSG00000000971 0.00 0.00 0.00 93.63
ENSG00000001036 37.29 34.21 37.29 36.32
ENSG00000001084 26.16 25.42 26.16 11.66
無事TPMのマトリックスを得ることができました。ちなみにgene id の" も必要ないので除去しています。
おわりに
これで無事、fastqファイルから遺伝子発現量マトリックスを生成できるようになりました。
あとはこのtsvファイルをBioconductorなどで解析すればばっちりですね。
細かい点では、gene idを遺伝子シンボルに変更したいなどの要求があると思いますが、機会があれば追記します。
