小云三分鐘教會你如何使用Trinity做轉(zhuǎn)錄本組裝（上）
Trinity, 是由Broad?Institute開發(fā)的進(jìn)行de?novo轉(zhuǎn)錄本組裝的工具，它能將RNA-seq?reads拼接成更長的轉(zhuǎn)錄本。Trinity由三個模塊組成，其基本作用如下所示：
Inchworm：將RNA-seq?reads組裝成單個轉(zhuǎn)錄本
Chrysails：將上一步形成的轉(zhuǎn)錄本進(jìn)行聚類，對于每一個聚類形成完整的de?Bruijin?graphs。每一個聚類代表一個基于轉(zhuǎn)錄本的復(fù)雜程度。
Butterfly：并行的處理每一個圖，追蹤圖中reads和配對reads的路徑，最后報道可變剪切亞型的全長轉(zhuǎn)錄本。
快來跟著小云一起學(xué)習(xí)吧！
1.?軟件安裝
Trinity需要在Linux系統(tǒng)上運(yùn)行。
軟件下載鏈接：https://github.com/trinityrnaseq/trinityrnaseq/releases
我們打開可以發(fā)現(xiàn)該軟件開發(fā)的很早，但仍然能保持每年更新的狀態(tài)，是非常難得的！我們選擇最新版本進(jìn)行下載，也就是圖中的trinityrnaseq-v2.15.1.FULL.tar.gz。
下載完成后我們使用如下命令進(jìn)行解壓安裝：
$?tar zxvf trinityrnaseq-v2.15.1.FULL.tar.gz
$?cd trinityrnaseq-v2.15.1
$?make
或者使用conda直接下載安裝：
$?conda install -c bioconda trinity
2.?Trinity使用方法
Trinity的主要參數(shù)如下所示：
# --seqType <string>: fast文件類型（fa或者fq）
# --max_memory <string>:?使用Trinity軟件時允許的最大內(nèi)存
# 當(dāng)數(shù)據(jù)為paired?reads時使用： --left <string>: 左端reads?--right?<string>: 右端reads
# 當(dāng)數(shù)據(jù)為unpaired?reads時使用： --single?<string>： 文件名
# --SS_lib_type <string>: 特異性文件文庫需要添加此參數(shù)，對于paired?reads有兩種參數(shù)（RF、FR）；對于unpaired?reads也有兩種（F、R）
# --CPU?<int>: 指定CPU的使用個數(shù)
跟著小云一起來看一個例子吧！
運(yùn)行一下代碼：
$?Trinity –-seqType fq –-max_memory 10G –-left reads.left.fq.gz –-right reads.right.fq.gz --SS_lib_type RF --CPU 4
運(yùn)行結(jié)束時會自動生成一個trinity_out_dir的文件夾，里面記載著所有輸出文件，如‘trinity_out_dir.Trinity.fasta’,?其中一個fasta條目如下所示：
>TRINITY_DN1000_c115_g5_i1 len=247 path=[31015:0-148 23018:149-246]
?AATCTTTTTTGGTATTGGCAGTACTGTGCTCTGGGTAGTGATTAGGGCAAAAGAAGACAC
ACAATAAAGAACCAGGTGTTAGACGTCAGCAAGTCAAGGCCTTGGTTCTCAGCAGACAGA
?AGACAGCCCTTCTCAATCCTCATCCCTTCCCTGAACAGACATGTCTTCTGCAAGCTTCTC
?CAAGTCAGTTGTTCACAGGAACATCATCAGAATAAATTTGAAATTATGATTAGTATCTGA
?TAAAGCA
‘TRINITY_DN1000_c115‘是聚類的標(biāo)號，’g5‘代表一個基因，’i1‘代表一個轉(zhuǎn)錄本。
因為是算法幫助我們組裝的轉(zhuǎn)錄本，我們怎么知道它組裝的好壞呢？接下來小云將繼續(xù)介紹如何用不同方法查看轉(zhuǎn)錄本的組裝質(zhì)量。
3.?計算全長轉(zhuǎn)錄本數(shù)量
其中一種用來衡量拼接質(zhì)量的方法就是數(shù)全長或者接近全長轉(zhuǎn)錄本的數(shù)量，越多代表組裝的效果越好。對于模式物種，我們可以將轉(zhuǎn)錄本比對到參考集上去，而對于非模式物種，如果存在有高質(zhì)量注釋的近源物種，我們可以使用它們的參考集進(jìn)行衡量。
我們使用Blast+進(jìn)行分析，首先構(gòu)建蛋白質(zhì)數(shù)據(jù)庫
$?makeblastdb -in uniprot_sprot.fasta -dbtype prot
然后進(jìn)行搜索
$?blastx -query Trinity.fasta -db uniprot_sprot.fasta -out blastx.outfmt6 -evalue 1e-10 -num_threads 6 -max_target_seqs 1 -outfmt 6
得到了輸出文件之后，我們可以使用Trinity軟件之中的腳本自動幫我們檢測Trinity轉(zhuǎn)錄本比對到參考組上的百分比
$?TRINITY_HOME/util/analyze_blastPlus_topHit_coverage.pl blastx.outfmt6 Trinity.fasta uniprot_sprot.fasta
跑完之后我們會得到一個’blastx.outfmt6.txt.w_pct_hit_length’文件，通過它我們可以畫出一個表格來更好的展示結(jié)果，如下：
看到圖表中的第二行，其含義為數(shù)據(jù)庫中有268個蛋白質(zhì)與Trinity轉(zhuǎn)錄本匹配，匹配長度為每個轉(zhuǎn)錄本長度的80-90%；數(shù)據(jù)庫中有3510個蛋白質(zhì)可以被幾乎為全長的轉(zhuǎn)錄本表示，比對覆蓋度大于80%。第一行的信息為，有3242個蛋白質(zhì)與轉(zhuǎn)錄本匹配，匹配長度均大于90%。
4.?RNA?Seq?Read?Representation
另外一種方法是為了檢測RNA?seq?reads能否較好的表示Trinity組裝的轉(zhuǎn)錄本，為了全面捕獲reads，我們使用Bowtie2將reads和轉(zhuǎn)錄組對齊，然后計算正確的不正確的和被孤立的reads的數(shù)量
首先，我們?yōu)檗D(zhuǎn)錄本建立索引
$?bowtie2-build Trinity.fasta Trinity.fasta
然后進(jìn)行比對來獲取read的比對信息,?我們以雙端數(shù)據(jù)為例
$?bowtie2 -p 10 -q --no-unal -k 20 -x Trinity.fasta -1 reads_1.fq -2 reads_2.fq 2>align_stats.txt| samtools view -@10 -Sb -o bowtie2.bam
最后可以查看輸出文件’align_stats.txt’
$?cat 2>&1 align_stats.txt
結(jié)果如下所示：
我們可以看到有89.97%的比對率，總的來說對于雙端數(shù)據(jù)，有70~80%的比對率就算合格了。
好了今天的分享就到這里啦~有什么問題可以和小云一起討論哦~下一期小云繼續(xù)為大家?guī)鞹rinity的更多后續(xù)內(nèi)容哦
更多生信小工具可查詢云平臺：http://www.biocloudservice.com/home.html
