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 ，尋找具有可變剪接的全長轉錄本，同時將旁系基因的轉錄本分開

 https://github.com/trinityrnaseq

Trinity的硬體需求：

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

使用的轉錄組資料為 Schizosaccharomyces pombe 

1. % wget \
2. http://sourceforge.net/projects/trinityrnaseq/files/misc/Trinity
3. NatureProtocolTutorial.tgz/download
4. #解壓後得到如下的檔案
5. tar –xvf TrinityNatureProtocolTutorial.tgz

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

1. % cat *.left.fq > reads.ALL.left.fq
2. % cat *.right.fq > reads.ALL.right.fq
3. #新增環境變數
4. % export PATH=/usr/local/tools:$PATH
5. #一種典型的使用方法入下
6. #其中引數SS_lib_type RF 表示資料是雙端（RF or FR) 單端（F or R）
7. % 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表示的意思如下圖。

1. % $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta

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

1. #首先準備GMAP需要的參考基因組，參考基因組檔案為genome.fa
2. gmap_build -d genome -D ./  
3. #algin 拼接結果，儲存為一個sam檔案
4. gmap -n 0 -D . -d genome ./trinity_out_dir/Trinity.fasta -f samse > trinity_gamp.sam

使用samtools轉換為BAM檔案（binary sam 優點是佔用磁碟空間小，運算速度快，一些對資料的排序或者提取命令需要轉換為BAM檔案)

1. samtools view -Sb trinity_gmap.sam > trinity_gmap.bam
2. #排序，方便後續使用
3. samtools sort trinity_gmap.bam trinity_gmap
4. #建立索引，需要先排序，否則報錯，產生.bai檔案
5. samtools index trinity_gmap.bam

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

1. #準備參考基因組
2. bowtie2-build GENOME_data/genome.fa genome
3. #run tophat 將所有的reads比對到參考基因組上
4. tophat2 -I 300 -i 20 genome \
5.      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 \
6.      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
7. #下面的IGV基因組瀏覽器需要先建立索引
8. samtools index tophat_out/accepted_hits.bam

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

1. 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還準備了一些指令碼進行後續的比如定量，差異表達等一些分析。

1.  #使用Trinity準備好的指令碼先用bowtie
2. #align到拼接好的轉錄組，然後使用RSEM定量
3. #執行這個指令碼後會產生兩個檔案 'Sp_ds.isoforms.results' and 'Sp_ds.genes.results'
4. #包含了Trinity 拼接的轉錄本（isoform) 和基因的raw counts數和標準化後的數值
5.  ${Trinity_home}/util/align_and_estimate_abundance.pl --seqType fq  \
6.       --left RNASEQ_data/Sp_plat.left.fq.gz --right RNASEQ_data/Sp_plat.right.fq.gz \
7.       --transcripts trinity_out_dir/Trinity.fasta \
8.       --output_prefix Sp_plat --est_method RSEM  --aln_method bowtie \
9.       --trinity_mode --prep_reference --output_dir Abundance_quantify/Sp_plat.RSEM
10. #然後再對其他三個樣本進行同樣的操作
11.  
12. #一個樣本間的比較矩陣 ,結果產生一個字尾為 .counts.matrix的檔案
13. #顯示了每個樣本在每個轉錄本（isoform)上的map的數目（raw count)
14. ${Trinity_home}/util/abundance_estimates_to_matrix.pl  --est_method RSEM --out_prefix Trinity_trans \
15.  Abundance_quantify/Sp_ds.RSEM/Sp_ds.isoforms.results \
16.  Abundance_quantify/Sp_hs.RSEM/Sp_hs.isoforms.results \
17.  Abundance_quantify/Sp_log.RSEM/Sp_log.isoforms.results \
18.  Abundance_quantify/Sp_plat.RSEM/Sp_plat.isoforms.results
19.  #另外 Trinity_trans.TMM.EXPR.matrix 是消除了測序深度，基因長度，然後通過TMM方法標準化後的數值(假定其他大多數基因沒有差異表達)

使用 EdgeR 分析差異表達基因

還是通過Trinity安裝包裡自帶的指令碼，不加引數執行會有基本引數的介紹

使用剛才獲得的 Trinity_trans.count.matrix 檔案

1. > ${Trinity_home}/Analysis/DifferentialExpression/run_DE_analysis.pl \
2. >  --matrix Trinity_trans.counts.matrix \
3. >  --method edgeR \
4. >  --dispersion 0.1 \
5. >  --output edgeR

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

logFC = log fold change

logCPM = log counts per million

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

Pre: 【轉載】Samtools安裝及使用