使用Trinity拼接以及分析差異表達一個小例子

Trinity 將測序數據分為許多獨立的de Brujin graph，理論上每一個圖對應一個表達的基因。

整個流程分為三個步驟：Inchworm, Chrysalis, and Butterfly

Inchworm: 從reads中提取所有的重疊k-mers，根據豐度遞減的順序檢查每個k-mers，然后將重疊的k-mers延長到不能再延長，稱為一個contig

Chrysalis: 將上一部生成的contig聚類，對每個類構建de Brujin graph

Butterfly: 根據構建的de Brujin graph ，尋找具有可變剪接的全長轉錄本，同時將旁系基因的轉錄本分開

Trinity的硬件需求：

Inchworm 和 Chrysails 步驟對內存的需求很大，官方給出的說法是大致為每一百萬對PE reads需要1g內存

使用的轉錄組數據為 Schizosaccharomyces pombe ，共4個樣本（left right 表示雙端測序數據的兩端）

在拼接時，可以將每個樣本都拼接成一個轉錄組，但是更合理的方法是將所有樣本的reads合在一起再進行拼接，所以先將這四個樣本的reads合在一起。

% cat *.left.fq > reads.ALL.left.fq
% cat *.right.fq > reads.ALL.right.fq
#添加環境變量
% export PATH=/usr/local/tools:$PATH
#一種典型的使用方法入下
#其中參數SS_lib_type RF 表示數據是雙端（RF or FR) 單端（F or R）
% Trinity --seqType fq --max_memory 1G --left reads.ALL.left.fq --right reads.ALL.right.fq --SS_lib_type RF --CPU 2

完成后會在當前的工作目錄生成一個 trinity_out_dir 的文件夾，Trinity.fasta為最終拼接結果。

Trinity自帶了一個腳本可以顯示一些結果的基本統計信息，N50表示的意思如下圖。

使用GMAP將拼接結果比對到參考基因組（有參考基因組的情況下）

#首先准備GMAP需要的參考基因組，參考基因組文件為genome.fa
gmap_build -d genome -D ./
#algin 拼接結果，保存為一個sam文件
gmap -n 0 -D . -d genome ./trinity_out_dir/Trinity.fasta -f samse > trinity_gamp.sam

使用samtools轉換為BAM文件（binary sam 優點是占用磁盤空間小，運算速度快，一些對數據的排序或者提取命令需要轉換為BAM文件)

使用tophat 將RNA-seq reads map到參考基因組

#准備參考基因組
bowtie2-build GENOME_data/genome.fa genome
#run tophat 將所有的reads比對到參考基因組上
tophat2 -I 300 -i 20 genome \
RNASEQ_data/Sp_log.left.fq.gz,RNASEQ_data/Sp_hs.left.fq.gz,RNASEQ_data/Sp_ds.left.fq.gz,RNASEQ_data/Sp_plat.left.fq.gz \
RNASEQ_data/Sp_log.right.fq.gz,RNASEQ_data/Sp_hs.right.fq.gz,RNASEQ_data/Sp_ds.right.fq.gz,RNASEQ_data/Sp_plat.right.fq.gz
#下面的IGV基因組瀏覽器需要先建立索引
samtools index tophat_out/accepted_hits.bam

使用基因組瀏覽器IGV （有GUI）查看trinity的拼接結果

igv.sh -g `pwd`/GENOME_data/genome.fa `pwd`/GENOME_data/genes.bed,`pwd`/tophat_out/accepted_hits.bam,`pwd`/trinity_gmap.bam

使用RSEM定量

除了拼接以外，Trinity還准備了一些腳本進行后續的比如定量，差異表達等一些分析。

#使用Trinity准備好的腳本先用bowtie
#align到拼接好的轉錄組，然后使用RSEM定量
#運行這個腳本后會產生兩個文件 'Sp_ds.isoforms.results' and 'Sp_ds.genes.results'
#包含了Trinity 拼接的轉錄本（isoform) 和基因的raw counts數和標准化后的數值
${Trinity_home}/util/align_and_estimate_abundance.pl --seqType fq \
--left RNASEQ_data/Sp_plat.left.fq.gz --right RNASEQ_data/Sp_plat.right.fq.gz \
--transcripts trinity_out_dir/Trinity.fasta \
--output_prefix Sp_plat --est_method RSEM --aln_method bowtie \
--trinity_mode --prep_reference --output_dir Abundance_quantify/Sp_plat.RSEM
#然后再對其他三個樣本進行同樣的操作
#一個樣本間的比較矩陣 ,結果產生一個后綴為 .counts.matrix的文件
#顯示了每個樣本在每個轉錄本（isoform)上的map的數目（raw count)
${Trinity_home}/util/abundance_estimates_to_matrix.pl --est_method RSEM --out_prefix Trinity_trans \
Abundance_quantify/Sp_ds.RSEM/Sp_ds.isoforms.results \
Abundance_quantify/Sp_hs.RSEM/Sp_hs.isoforms.results \
Abundance_quantify/Sp_log.RSEM/Sp_log.isoforms.results \
Abundance_quantify/Sp_plat.RSEM/Sp_plat.isoforms.results
#另外 Trinity_trans.TMM.EXPR.matrix 是消除了測序深度，基因長度，然后通過TMM方法標准化后的數值(假定其他大多數基因沒有差異表達)

使用 EdgeR 分析差異表達基因

還是通過Trinity安裝包里自帶的腳本，不加參數運行會有基本參數的介紹

使用剛才獲得的 Trinity_trans.count.matrix 文件

運行結果 '*.DE_results' 輸出了運行edgeR 分離出來的差異表達的基因

logFC = log fold change

logCPM = log counts per million

#提取FDR<=0.005）
sed '1,1d' edgeR/Trinity_trans.counts.matrix.Sp_log_vs_Sp_plat.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' | wc -l
#畫熱圖，需要進入剛才的/edgeR文件夾作為工作目錄
$TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl \
--matrix ../Trinity_trans.TMM.EXPR.matrix -P 1e-3 -C 2
#-P 為p的閾值，-C 為fold change = 2^2 =4 倍