top of page

5. リードカウントのマトリクスを作成する

2021.06.10 初稿

​

 

こちらのページを参考にし、文章を抜粋したものが以下。

発現量の評価にはTPM、edgeRの入力にはcountを使おう

遺伝子発現量を直接見たい場合、例えば box plot や scatter plot を描く場合には TPMを使っておけば基本的に間違いありません。 一方、edgeRやDESeq2のような発現変動比較ツールにはTPMではなくcountデータを利用することが推奨されます（マニュアル参照）。何故なら、edgeRやDESeq2の中で用いている正規分布や負の二項分布はraw countに対してモデル化するものであり、遺伝子長などで正規化してしまうとこの分布を満たさなくなってしまうおそれがあるからです。全ての正規化はedgeRやDESeq2の中で行われますので、前もって正規化しておく必要はありません。更に言えば、サンプル間比較は遺伝子ごとに行われますので、基本的に遺伝子長の正規化は必要ないのです。

 

ということなのでリードカウントマトリクスを作る必要がある。

expected_countの列のみを抽出するため、RSEMの中にあるrsem-generate-data-matrixを使用する。

遺伝研スパコンの中にはsingularityで使えるRSEMが揃っている。

一番最新のバージョンは /usr/local/biotools/r/rsem:1.3.3--pl526ha52163a_0なので、これを使用する。

singularity exec -e /usr/local/biotools/r/rsem:1.3.3--pl526ha52163a_0 \

rsem-generate-data-matrix \

sample1_RSEM.isoforms.results \

sample2_RSEM.isoforms.results \

sample3_RSEM.isoforms.results \

> expected_count_matrix.tsv

 

上記ではsample1~3のalign_and_estimateで得たマトリクスのファイル名をそのまま載せているが、本当はサンプル名をリスト（.txt）にしたものをcat等で読み込む操作をはさみたい。サンプルがたくさん（10以上）ある時もあるので。しかし、リストにして読み込む方法を探ったけれどパッとは見つからなかった。今後良い方法を見つけたら追記します。

 

 

TPMのマトリクスを作成するには、trinityrnaseqの中に含まれているabundance_estimates_to_matrix.plを用いる。

--quant_filesオプションでTPM列を持っているファイル名を指定する必要があるため、事前にファイル名をリスト化したファイル.txtを作っておく。

$ ls ../*/RSEM.isoforms.results > ./list_sample.txt

$ less list_sample.txt

~/RNAseq/sample1_rep1/RSEM.isoforms.results

~/RNAseq/sample1_rep2/RSEM.isoforms.results

~/RNAseq/sample1_rep3/RSEM.isoforms.results

 

 

---------------------------------------------

"abundance_estimates.sh"の中身

 

TRINITY_HOME=/usr/local/bin/trinityrnaseq

singularity exec -e ${HOME}/trinityrnaseq.v2.11.0.simg \

${TRINITY_HOME}/util/abundance_estimates_to_matrix.pl \

--est_method RSEM \

--gene_trans_map ${HOME}/RNAseq/trinity_out/sample.Trinity.fasta.gene_trans_map \ 

--name_sample_by_basedir \ 

--quant_files ${HOME}/RNAseq/align_and_estimate/count_matrix/list_sample.txt 

 

---------------------------------------------

 

このスクリプトを走らせると以下のようにたくさんのファイルが生成されるが、TPMの列を抽出してまとめたものはRSEM.gene.TPM.not_cross_norm。

 RSEM.gene.TMM.EXPR.matrix

 RSEM.gene.TPM.not_cross_norm

 RSEM.gene.TPM.not_cross_norm.TMM_info.txt

 RSEM.gene.TPM.not_cross_norm.runTMM.R

 RSEM.gene.counts.matrix

（geneをisoformに置換したものも同数生成される。）

​

​

参考にしたページ

Estimating transcript abundance （Trinity Github）：https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Transcript-Quantification

​

bottom of page