![[Share]Twitter](/common/img/lab/sns_tw.png){kind=small}
![[Share]Facebook](/common/img/lab/sns_fb.png){kind=small}
発現量解析をできるだけ自動化する
コマンドラインでは、実行手順を書き込んだファイルを用意することで、多くの作業を自動化できる。自動化のための作業手順
必要なフォルダーを予め作っておく
:> mkdir -p count/RSEM_boiwtie2/count_matrix/DE
:> mkdir -p count/RSEM_boiwtie2/count_matrix/QC
:> mkdir -p count/salmon/count_matrix/DE
:> mkdir -p count/salmon/count_matrix/QC
コマンドを書き込んだスクリプトファイルを、それぞれのフォルダーに置く
カウントファイル作成
:> touch count/RSEM_bowtie2/est_RSEM.sh
:> echo 'perl $TRINITY_HOME/util/align_and_estimate_abundance.pl --thread_count 16 --transcripts ./Trinity.fasta --trinity_mode --seqType fq --prep_reference --samples_file ../BRR_file.txt --est_method RSEM --aln_method bowtie2' >> count/RSEM_bowtie2/est_RSEM.sh
:> touch count/salmon/est_salmon.sh
:> echo 'perl $TRINITY_HOME/util/align_and_estimate_abundance.pl --thread_count 16 --transcripts ./Trinity.fasta --trinity_mode --seqType fq --prep_reference --samples_file ../BRR_file.txt --est_method salmon' >> count/salmon/est_salmon.sh
:> touch count/RSEM_boiwtie2/count_matrix/abundance_matrix.sh
:> echo 'perl $TRINITY_HOME/util/abundance_estimates_to_matrix.pl --est_method RSEM --gene_trans_map ../../../Trinity.fasta.gene_trans_map --name_sample_by_basedir --quant_files ./list_results.txt' >> abundance_matrix.sh
:> touch count/salmon/count_matrix/abundance_matrix.sh
:> echo 'perl $TRINITY_HOME/util/abundance_estimates_to_matrix.pl --est_method salmon --gene_trans_map ../../../Trinity.fasta.gene_trans_map --name_sample_by_basedir --quant_files ./list_sf.txt' >> abundance_matrix.sh
カウントのQC（RSEM）
:> touch count/RSEM_boiwtie2/count_matrix/QC/compare.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../RSEM.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --compare_replicates' >> compare.sh
:> touch count/RSEM_boiwtie2/count_matrix/QC/cpm.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../RSEM.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --CPM --sample_cor_matrix' >> cpm.sh
:> touch count/RSEM_boiwtie2/count_matrix/QC/pca.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../RSEM.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --CPM --center_rows --prin_comp 3' >> pca.sh
カウントのQC（salmon）
:> touch count/salmon/count_matrix/QC/compare.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../salmon.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --compare_replicates' >> compare.sh
:> touch count/salmon/count_matrix/QC/cpm.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../salmon.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --CPM --sample_cor_matrix' >> cpm.sh
:> touch count/salmon/count_matrix/QC/pca.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/PtR --matrix ../salmon.isoform.counts.matrix --samples ../BRR_DE.txt --log2 --min_rowSums 10 --CPM --center_rows --prin_comp 3' >> pca.sh
発現量の統計解析（RSEM）
:> touch count/RSEM_boiwtie2/count_matrix/DE/DE_DESeq2.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method DESeq2 --samples_file ../BRR_DE.txt --matrix ../RSEM.isoform.counts.matrix' >> DE_DESeq2.sh
:> touch count/RSEM_boiwtie2/count_matrix/DE/DE_edgeR.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method edgeR --samples_file ../BRR_DE.txt --matrix ../RSEM.isoform.counts.matrix' >> DE_edgeR.sh
:> touch count/RSEM_boiwtie2/count_matrix/DE/DE_voom.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method voom --samples_file ../BRR_DE.txt --matrix ../RSEM.isoform.counts.matrix' >> DE_voom.sh
:> touch count/RSEM_boiwtie2/count_matrix/DE/diff_P1e-2_C1.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../../RSEM.isoform.TMM.EXPR.matrix --samples ../../BRR_DE.txt -P 1e-2 -C 1' >> diff_P1e-2_C1.sh
この例では、1%水準で2倍の変化があった遺伝子を抽出している。候補が多すぎる場合など、研究の目的に応じてパラメーターを調整する。
発現量の統計解析（salmon）
:> touch count/salmon/count_matrix/DE/DE_DESeq2.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method DESeq2 --samples_file ../BRR_DE.txt --matrix ../salmon.isoform.counts.matrix' >> DE_DESeq2.sh
:> touch count/salmon/count_matrix/DE/DE_edgeR.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method edgeR --samples_file ../BRR_DE.txt --matrix ../salmon.isoform.counts.matrix' >> DE_edgeR.sh
:> touch count/salmon/count_matrix/DE/DE_voom.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --method voom --samples_file ../BRR_DE.txt --matrix ../salmon.isoform.counts.matrix' >> DE_voom.sh
:> touch count/salmon/count_matrix/DE/diff_P1e-2_C1.sh
:> echo '$TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../../salmon.isoform.TMM.EXPR.matrix --samples ../../BRR_DE.txt -P 1e-2 -C 1' >> diff_P1e-2_C1.sh
この例では、1%水準で2倍の変化があった遺伝子を抽出している。候補が多すぎる場合など、研究の目的に応じてパラメーターを調整する。
用意したスクリプトを順次連続で実行するスクリプトを作る
RSEM用
:> touch count/analysis_RSEM.sh
:> echo 'cd ./RSEM_bowtie2 && cp ../../Trinity.fasta ./ && /bin/sh ./est_RSEM.sh && cd ./count_matrix && ls ../*/RSEM.isoforms.results > ./list_sf.txt && /bin/sh ./abundance_matrix.sh && cat ../../BRR_file.txt | awk '{print $1,$2}' > ./BRR_DE.txt && cd ./QC && /bin/sh ./compare.sh && /bin/sh ./cpm.sh && /bin/sh ./pca.sh && cd ../DE && /bin/sh ./DE_edgeR.sh && /bin/sh ./DE_DESeq2.sh && /bin/sh ./DE_voom.sh && cd ./edgeR*dir && /bin/sh ../diff_P1e-2_C1.sh && cd ../DESeq2*dir && /bin/sh ../diff_P1e-2_C1.sh && cd ../voom*dir && /bin/sh ../diff_P1e-2_C1.sh' >> count/analysis_RSEM.sh
salmon用
:> touch count/analysis_salmon.sh
:> echo 'cd ./salmon && cp ../../Trinity.fasta ./ && /bin/sh ./est_salmon.sh && cd ./count_matrix && ls ../*/quant.sf > ./list_sf.txt && /bin/sh ./abundance_matrix.sh && cat ../../BRR_file.txt | awk '{print $1,$2}' > ./BRR_DE.txt && cd ./QC && /bin/sh ./compare.sh && /bin/sh ./cpm.sh && /bin/sh ./pca.sh && cd ../DE && /bin/sh ./DE_edgeR.sh && /bin/sh ./DE_DESeq2.sh && /bin/sh ./DE_voom.sh && cd ./edgeR*dir && /bin/sh ../diff_P1e-2_C1.sh && cd ../DESeq2*dir && /bin/sh ../diff_P1e-2_C1.sh && cd ../voom*dir && /bin/sh ../diff_P1e-2_C1.sh' >> count/analysis_salmon.sh
スクリプファイルの実行とディレクトリの移動を「&&」でつないでいる。
「&&」でつなぐと、処理が完了したら次の処理、と順次実行してくれる。
自動化スクリプトを実行する
上記の準備作業を終えていたら、以後はTrinityのアセンブル終了後に下記の作業を行うことで、定型の発現量解析が完了する。スクリプトの実行が完了したら、edgeR.数字.dir に移動して diff_.sh を実行する。QCの結果やDEの内容を見て、閾値を変更して diff_.sh を作成して詳細な解析を行うと良い。
count フォルダーを Trinity の結果フォルダーにコピーする
:> cd /解析したいTrinity結果フォルダーのパス/count
:> /bin/sh ./analysis_salmon.sh && cd /パス/count && /bin/sh ./analysis_RSEM.sh
スクリプトファイル（.sh という拡張子をつけたファイル）を作成する意味
スクリプトファイルを作成して作業したフォルダー内に残しておくことで、作業ログにもなる。そのフォルダー内に残された解析結果が、どのようなコマンドによって得られたものなのか確実な記録になるので、完璧な再現性を保証する研究ノート代わりにもなる。紙のノートに手書きする際に発生する転記ミスを防ぐ意味もあるので、解析研究での研究ノートの残し方としてお薦めする。
スクリプトは「#」から改行コードまでをメモとして扱い、コマンドには影響しないという仕組みがあるので、上記のスクリプ内に実行した日付や研究の条件などを書き込んでおくと、研究ノートとしての情報量が上がってより有益である。
例
# 2020/05/05
# ナミアゲハのRNA-seq
# Trinity 2.10.0 を使用
目次
- 環境構築
- リードのクオリティチェックとアセンブル
- Transdecoder による遺伝子予測とアノテーション
- リード数のカウント
- カウントマトリクスの作成
- リードカウントのQC解析
- DE解析
- DEG取り出し
- おまけ1: Transdecoderの自動化
- おまけ2: 発現量解析の自動化
- おまけ3: Anacondaコマンドの使い方一覧