使用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 ,尋找具有可變剪接的全長轉錄本,同時將旁系基因的轉錄本分開
https://github.com/trinityrnaseq
Trinity的硬體需求:
Inchworm 和 Chrysails 步驟對記憶體的需求很大,官方給出的說法是大致為每一百萬對PE reads需要1g記憶體
使用的轉錄組資料為 Schizosaccharomyces pombe
- % wget \
- http://sourceforge.net/projects/trinityrnaseq/files/misc/Trinity
- NatureProtocolTutorial.tgz/download
- #解壓後得到如下的檔案
- tar –xvf TrinityNatureProtocolTutorial.tgz
在拼接時,可以將每個樣本都拼接成一個轉錄組,但是更合理的方法是將所有樣本的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表示的意思如下圖。
- % $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta
使用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檔案)
- samtools view -Sb trinity_gmap.sam > trinity_gmap.bam
- #排序,方便後續使用
- samtools sort trinity_gmap.bam trinity_gmap
- #建立索引,需要先排序,否則報錯,產生.bai檔案
- samtools index trinity_gmap.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 檔案
- > ${Trinity_home}/Analysis/DifferentialExpression/run_DE_analysis.pl \
- > --matrix Trinity_trans.counts.matrix \
- > --method edgeR \
- > --dispersion 0.1 \
- > --output edgeR
執行結果 '*.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 倍
Pre: 【轉載】Samtools安裝及使用