![[Share]Twitter](/common/img/lab/sns_tw.png){kind=small}
![[Share]Facebook](/common/img/lab/sns_fb.png){kind=small}
2023/03/28 誤表記を修正「1のマイナス3乗＝0.001」→「1 × 10のマイナス3乗 ＝ 0.001」
DE（発現量比較）解析
統計解析用プログラム言語「R (アール)」のパッケージを使って、発現量の変動を統計学的に解析する。作業手順
一つ上のフォルダー（count_matrix）に移動する（ひとつ前の手順で cout_matrix/QC に移動しているので、そこにいるという前提で解説している）
:> cd ../
count_matrix フォルダーの中で、新しいフォルダーを作る
:> mkdir ./DE
フォルダーに移動
:> cd ./DE
解析に使用する BRR_DE.txt は、ひとつ前の手順のQCと共通なので、作成を繰り返す必要はない。
スクリプトファイルを作る
:> touch ./DE_edgeR.sh
スクリプトファイルを開いて編集
:>$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method edgeR --samples_file ../BRR_DE.txt --matrix ../salmon.isoform.counts.matrix
（RSEMの場合は、salmonのところを書き換える）
スクリプトを実行
:> /bin/sh ./DE_edgeR.sh
統計解析のツールは、edgeR の他に DESeq2 と voom が使える。上記スクリプトの「edgeR」の部分を「DESeq2」「voom」に変更する。
スクリプトファイルを作る
:> touch ./diff_P1e-2_C1.sh
スクリプトファイルを開いて編集
$TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../../salmon.isoform.TMM.EXPR.matrix --samples ../../BRR_DE.txt -P 1e-2 -C 1
-P: p値（有意差検定の値）の閾値。1e-2 の場合は0.01なので、1%水準で「発現が上昇している」という仮説が支持されるかを示している。
-C: 2の何乗の倍率で発現が上昇しているか。この場合2の1乗＝2倍。
この二つのオプションの組み合わせで、「二つのサンプル間で2倍以上発現が変化しているという仮説が1％水準で支持される遺伝子を検出する」という意味になる。5%水準であれば「-P 0.05」、0.1%水準であれば「1e-3」と設定する。4倍の変化で検出したい場合は「-C 2（2の2乗）」、8倍の変化で検出したい場合は「-C 3（2の3乗）」とする。
どれくらいの厳しさで比較するかは、結果を見ながら適宜調整して各自の研究目的に最適なものを見つけること。
DE解析結果のフォルダーに移動
:> cd ./edgeR.数字.dir
（結果ファイルがあるフォルダーで動かすプログラム）
スクリプトを実行
:> /bin/sh ../diff_P1e-2_C1.sh
以後、DESeq2 と voom についても同様に、DESeq2.数字.dir または voom.数字.dir に移動して、「/bin/sh ../diff.sh」を実行する。
結果
ファイル名の最後に「サンプル名-UP.subset」と付いているファイルが、比較対象に対して設定した条件で発現が上昇している遺伝子のリスト。この結果がどれくらい信頼できるかについては、edgeRとvoomの結果についてはFDR（False Discovery Rate）、DESeq2の結果についてはpadj（補正したp値）を見て判断する。padjは、解析パッケージによってはadjPとい表記されている場合もあるが、意味は同じ。FDRやpadjの値は小さいほど良い。FDRやpadjは「1e-3」のような表記がされているが、「e-」はマイナス何乗なのかを表している。「1e-3」の場合、「1 x 10のマイナス3乗＝0.001」という意味。「e-」の後ろの数字が大きいほど、「0.~」の0の数が多くなるので値としては小さい。
スクリプトファイルを作成して作業したフォルダー内に残しておくことで、作業ログにもなる。そのフォルダー内に残された解析結果が、どのようなコマンドによって得られたものなのか確実な記録になるので、完璧な再現性を保証する研究ノート代わりにもなる。紙のノートに手書きする際に発生する転記ミスを防ぐ意味もあるので、解析研究での研究ノートの残し方としてお薦めする。
スクリプトは「#」から改行コードまでをメモとして扱い、コマンドには影響しないという仕組みがあるので、上記のスクリプ内に実行した日付や研究の条件などを書き込んでおくと、研究ノートとしての情報量が上がってより有益である。
例
# 2020/05/05
# ナミアゲハのRNA-seq
どの解析パッケージを使うか
全部試せ手間も計算時間もたいしたことない。
この条件だけ試せばどの研究にもOK！なんて設定はない。
DESeq2, edgR, voom では微妙に結果が異なる。全体的に、DESeq2ではedgeRやvoomより抽出される遺伝子の数が多くなる傾向があるが、どちらが優れているということはない。研究目的や条件によって、検出される遺伝子の数が多すぎたり少なすぎたりすることもあるので、3種類のパッケージとp値や倍率を変えながらじっくり検証して、自分の研究目的に合う結果を丁寧に探索すること。
目次
- 環境構築
- リードのクオリティチェックとアセンブル
- Transdecoder による遺伝子予測とアノテーション
- リード数のカウント
- カウントマトリクスの作成
- リードカウントのQC解析
- DE解析
- DEG取り出し
- おまけ1: Transdecoderの自動化
- おまけ2: 発現量解析の自動化
- おまけ3: Anacondaコマンドの使い方一覧