1. Trinity進行轉錄組組裝
Trinity進行轉錄組組裝的典型命令如下:
$ /opt/biosoft/trinityrnaseq_r20131110/Trinity.pl --seqType fq --JM 50G\ --left sample1_1.clean.fastq sample2_1.clean.fastq\ --right sample1_2.clean.fastq sample2_2.clean.fastq\ --jaccard_clip --CPU 6 --SS_lib_type FR
–JM后的參數設定與轉錄組的大小有關，在內存足夠的情況下，設定大點能節約時間；
–left 和 –right后可以接多個樣平的數據，並用空格隔開，值得注意的是，left reads name以/1結尾，rigth reads name以/2結尾；
–jaccard_clip 適合於基因稠密的真菌物種；
–SS_lib_type 適合於鏈特異性測序
大數據量(>300M pairs)的RNA-seq數據，最好使用TRINITY_RNASEQ_ROOT/util/normalize_by_kmer_coverage.pl對reads進行處理后再使用trinity進行組裝，以降低內存消耗和大量時間。
也可以設置–min_kmer_cov 2，丟棄uniquely occurring kmer, 從而降低內存消耗。
參考文獻：
1. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011 May 15;29(7):644-52. doi: 10.1038/nbt.1883. PubMed PMID: 21572440.
2. Borodina T, Adjaye J, Sultan M. A strand-specific library preparation protocol for RNA sequencing. Methods Enzymol. 2011;500:79-98. PubMed PMID: 21943893.
2. Trinity輸出結果的統計
Trinity默認的輸出結果為：trinity_out_dir/Trinity.fasta。
該fasta格式文件中序列名例如：
>comp6749_c0_seq1 len=328 path=[471:0-83 388:84-208 679:209-327] >comp6749_c0_seq2 len=328 path=[304:0-83 388:84-208 679:209-327] >comp6749_c0_seq3 len=245 path=[901:0-125 679:126-244]
可以看到，trinity生成的結果為components,而一個components可能有多個seq。這相當於一個gene能有多個transcripts。
可以使用trinity自帶的程序TrinityStats.pl對components和transcripts的數目，大小和N50等進行統計。
$ $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta Total trinity transcripts: 40138 Total trinity components: 31067 Percent GC: 61.31
3. 將reads比對到轉錄組,並進行可視化
TRINITY_RNASEQ_ROOT/util/alignReads.pl能調用bowtie將reads map到轉錄組，並可以設置鏈特異性參數。
$ TRINITY_RNASEQ_ROOT/util/alignReads.pl --left left.fq --right right.fq --seqType fq\ --target Trinity.fasta --aligner bowtie --retain_intermediate_files
結果中生成coordSorted和nameSorted的sam和bam文件。如果設置了鏈特異性參數，則額外生成+鏈和-鏈的比對結果文件。
TRINITY_RNASEQ_ROOT/util/SAM_nameSorted_to_uniq_count_stats.pl用於統計比對結果
$ $TRINITY_HOME/util/SAM_nameSorted_to_uniq_count_stats.pl bowtie_out.nameSorted.sam.+.sam #read_type count pct proper_pairs 21194964 93.22 both read pairs align to a single contig and point toward each other. left_only 836213 3.68 only the left (/1) read is reported in an alignment right_only 687576 3.02 only the right (/2) read is reported in an alignment improper_pairs 16640 0.07 both left and right reads align, but to separate contigs, or to a single contig in the wrong expected relative orientations.
可以將Trinity.fasta導入到IGV中作為genome，上載bam文件，從而可視化比對結果。
4. 使用RSEM進行表達量計算
首先，需要下載最新版本的RSEM，安裝並將程序加入到$PATH中。
$ wget http://deweylab.biostat.wisc.edu/rsem/src/rsem-1.2.8.tar.gz $ tar zxf rsem-1.2.8.tar.gz $ cd rsem-1.2.8 $ make $ echo "PATH=$PWD:\$PATH" >> ~/.bashrc
使用$TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl可以調用RSEM，從而計算表達量。如果是鏈特異性測序，則加入–SS_lib_type參數。
$TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl --transcripts Trinity.fasta \
--seqType fq --left left.reads.fq --right right.reads.fq --SS_lib_type FR \
--prefix RSEM --thread_count 4 -- --bowtie-phred64-quals --no-bam-output
將rsem-calculate-expression程序的參數–bowtie-phred64-quals和–no-bam-output加入到run_RSEM_align_n_estimate.pl中，則如上所示。這兩個參數分別代表fastq的質量格式是phred64，不輸出bam文件(節約大量時間)。
若運行出現問題，點擊：RSEM的README文件。
結果生成兩個abundance estimation information文件:
RSEM.isoforms.results : EM read counts per Trinity transcript
RSEM.genes.results : EM read counts on a per-Trinity-component (aka… gene) basis, ‘gene’ used loosely here.
可以根據得到的結果，去除掉IsoPct低於1%的transcripts。可以依據RSEM.isoforms.results使用TRINITY_RNASEQ_ROOT/util/filter_fasta_by_rsem_values.pl過濾掉trinity組裝結果中的lowly supported transcripts。
但不推薦過濾掉這些序列。
5. 鑒定差異表達transcripts
Trinity可以使用Bioconductor package中的edgeR或DESeq來鑒定差異表達trancripts。因此，需要安裝R和相關的一些包。
source("http://bioconductor.org/biocLite.R")
biocLite('edgeR')
biocLite('DESeq')
biocLite('ctc')
biocLite('Biobase')
install.packages('gplots’)
install.packages(‘ape’)
5.1 使用上一節中的RSEM來分別對每個樣品的每個生物學重復進行表達量計算
5.2 將每個樣的RSEM的結果進行合並
$ $TRINITY_HOME/util/RSEM_util/merge_RSEM_frag_counts_single_table.pl \ sampleA.RSEM.isoform.results sampleB.RSEM.isoform.results ... \ > transcripts.counts.matrix $ TRINITY_HOME/util/RSEM_util/merge_RSEM_frag_counts_single_table.pl \ sampleA.RSEM.gene.results sampleB.RSEM.gene.results ... \ > genes.counts.matrix
然后修改生成的兩個matrix文件的column headers（代表着樣品和重復的名字），有利於下游的分析。如果要分析transcripts水平的差異表達，則使用transcripts.counts.matrix文件；若要分析gene水平的差異表達，則使用genes.counts.matrix。
5.3 無生物學重復進行差異表達分析
$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl用於調用edgeR或DESeq進行差異表達基因分析。直接輸入該命令查看其用法。
Trinty推薦使用edgeR進行差異表達分析。
$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl \ --matrix counts.matrix --method edgeR
注意輸入的matrix是counts的數據，而不要是FPKM的數據。
5.4 有生物學重復進行差異表達分析
首先，要建立文件samples_described.txt，內容為：
conditionA condA-rep1 conditionA condA-rep2 conditionB condB-rep1 conditionB condB-rep2 conditionC condC-rep1 conditionC condC-rep2
condA-rep1, condA-rep2, condB-rep1… 等對應着counts.matrix文件中的column names。
命令如下：
$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl \ --matrix SP2.rnaseq.counts.matrix --method edgeR \ --samples_file samples_described.txt
結果文件中 logFC 是 log2 Fold Change; logCPM 是 log2-counts-per-million。
值得注意的是：程序默認去除counts數都少於10的transcripts或genes，不對其進行差異分析。所以有差異分析的genes或transcripts數目低於原始的數目。
5.5 提取差異表達基因，對其進行聚類分析
5.5.1 表達量的 NORMALIZED
使用TMM方法將counts轉換為FPKM。
首先從1個樣平的RSEM結果中提取長度數據：
$ cut -f 1,3,4 sampleA.RSEM.isoforms.results > feature_lengths.txt
然后使用TMM方法將counts數據轉換為FPKM數據:
$ $TRINITY_HOME/Analysis/DifferentialExpression/run_TMM_normalization_write_FPKM_matrix.pl \ --matrix counts.matrix --lengths feature_lengths.txt
5.5.2 提取差異表達轉錄子
注意的是，這一步要在edgeR的結果文件中運行程序：
$ $TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl \ --matrix matrix.TMM_normalized.FPKM -P 0.001 -C 2
默認下選擇FDR值低於0.001，log2fold-change的絕對值>=2為差異表達基因。
程序輸出差異表達基因FPKM、log2FC、FDR等值 和 聚類圖 Heat Map.
5.5.3 根據聚類圖提取子類
根據聚類結果，可以自動或手動確定子類。
自動確定子類：
$ $TRINITY_HOME/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl \ --Ptree 20 -R file.all.RData
上例中從數的20%處來自動划分子類。
手動確定子類：
$ R
> load("all.RData") # check for your corresponding .RData file name to use here, replace all.RData accordingly
> source("$TRINITY_HOME/Analysis/DifferentialExpression/R/manually_define_clusters.R")
> manually_define_clusters(hc_genes, centered_data)
然后左鍵點擊選擇子類，右鍵結束選擇
6. 提取蛋白編碼區
使用transdecoder從trinity的轉錄子中提取coding region。最新版的transdecoder貌似有點問題。
$ $TRINITY_HOME/trinity-plugins/transdecoder/transcripts_to_best_scoring_ORFs.pl \ -t transcripts.fasta -m 100
默認下允許的最小的protein長度為100.
提取出了coding region，得出對應的protein序列，有利於於下一步的功能注釋。
原文出自：http://www.chenlianfu.com/?p=2026