在上一期的教學中小云講到了如何下載和使用Trinity組裝轉錄本，并且評估轉錄本的質量，在本期中我們將當組裝完成后可以進行那些下游分析來挖掘我們關注的生物學問題。

跟著小云一起來看看吧

1.?轉錄本豐度

在進行其它分析之前，我們通常學對轉錄本和基因的豐度進行定量。在Trinity中提供了調用三種不同軟件的腳本，RSEM，kallisto和salmon。值得注意的是Trinity并沒有同時打包安裝好這三個軟件，需要我們自己下載并配置好相對應的環(huán)境變量。

我們使用一下面的腳本perl腳本進行計算

$?TRINITY_HOME/util/align_and_estimate_abundance.pl

主要參數如下：

# --transcripts <string>: 轉錄本文件

# --seqType <string>: fq|fa

# 若為paired-end數據： --left <string>: 左端reads文件 --right <string>: 右端reads文件

# single-end數據：--single <string>: 原始reads文件

# --est_method <string>: RSEM|kallisto|salmon

比如說我們RSEM來進行定量，相關命令如下：

$ TRINITY_HOME/util/align_and_estimate_abundance.pl --transcript Trinity.fasta --est_method RSEM --aln_method bowtie --seqType fq --left reads_1.fq --right reads_2.fq --trinity_mode --prep_reference --output_dir rsem_outdir

之后會返還兩個輸出文件：

’RSEM.isoforms.results’、’RSEM.genes.results’

基因分型的結果如下所示：

?

至此，我們可以對每個樣本構建轉錄本和基因的表達矩陣，使用下面這個腳本：

$?TRINITY_HOME/util/abundance_estimates_to_matrix.pl

# --est_method <string>: ?RSEM|kallisto|salmon

# --gene_trans_map <string>: 基因和轉錄本的轉換文件

例子如下：

$?TRINITY_HOME/util/abundance_estimates_to_matrix.pl --est_method RSEM --gene_trans_map Trinity.fasta.gene_trans_map --name_sample_by_basedir

我們會得到以下三個文件：

‘RSEM.isoform.count.matrix’: 原始counts數據

‘RSEM.isoform.TPM.not_cross_norm’: TPM表達量（沒有進行跨樣本歸一化）

‘RSEM.isoform.TMM.EXPR.matrix’: TMM歸一化的表達值

?

2.?差異表達分析

得到了表達矩陣之后我們就可以進行差異表達分析了，Trinity中使用的是EdgeR這個R包來做差異表達分析，同樣需要咱們的系統(tǒng)中裝有相對應的R環(huán)境和軟件包。

具體腳本命令和參數如下所示：

$?TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl

# --matrix|m <string>: 原始count矩陣

# --method <string>: edgeR|DESeq2|voom

運行以下腳本：

$?TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix RSEM.isoform.count.matrix --method edgeR?–samples_file samples_described.txt

關于edgeR的參數設置可以查詢edgeR manual。

然后我們會得到三個文件，分別是：

‘${prefix}.sampleA_vs_sampleB.${method}.DE_results’

‘${prefix}.sampleA_vs_sampleB.${method}.MA_n_Volcano.pdf’

${prefix}.sampleA_vs_sampleB.${method}.Rscript’

?

我們一起來看看DE_results文件的頭幾行，記錄著log?fold?change和相對應的統(tǒng)計量（Pvalue，FDR矯正）

同時還有火山圖

對于差異表達的一個比較初始的分析就是去獲取那些前后表達差異極大的轉錄本（最顯著的FDR和fold-change），然后根據不同的表達模式將這些轉錄本劃分聚類。我們可以使用以下腳本進行聚類

$?TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl

# --matrix <string>: TMM.EXPR.matrix

# -P <float>: pvalue閾值

# -C <float>: 最小的fold change數

如小云下面舉的例子：

$?TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../Trinity_trans.TMM .EXPR.matrix -P 1e-3 -C 2

然后我們就能得到如下差異基因和不同樣本的熱圖：

3.?轉錄本功能注釋

在我們得到了拼接的轉錄本以及相應的基因集之后，我們自然而然想知道這些基因富集在什么樣的功能上。同樣的，咱們強大的Trinity提供了一套完整的流程。首先我們通過以下腳本來獲取每個基因所對應的GO term：

$?TRINITY_HOME/util/extract_GO_assignments_from_Trinotate_xls.pl --Trinoate_xls trinotate.xls -G --include_ancestral_terms > go_annotations.txt

然后我們再使用Bioconductor中的包GOseq進行功能富集：

$?TRINITY_HOME/Analysis/DifferentialExpression/run_GOseq.pl --factor_labeling factor_labeling.txt --GO_assignments go_annotations.txt --lengths gene.lengths.txt

其中’factor_labeling.txt’這個文件記載著不同樣本間的差異表達基因集，如下圖所示：

而對于‘gene.lengths.txt’這個文件，我們使用如下命令進行創(chuàng)建：

$?TRINITY_HOME/util/misc/fasta_seq_length.pl Trinity.fasta > Trinity.fasta.seq_lengths

$?TRINITY_HOME/util/misc/TPM_weighted_gene_length.py --gene_trans_map trinity_out_dir/Trinity.fasta.gene_trans_map --trans_lengths Trinity.fasta.seq_lens --TPM_matrix isoforms.TMM.EXPR.matrix > Trinity.gene_lengths.txt

?

在最后我們就會得到一個如下的圖標，其中記錄著差異表達基因富集到的功能通路：

?