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序列，有利于于下一步的功能注释。
陈博士您好，我是RNA-SEQ的bioinformatics的新手，主要做数据分析。可能会处理一些来自于Illumina和Ion torrent的可能是病毒方面的数据。现在还没拿到数据都是跟着网上论坛瞎混自己先练练手。
现在有几个问题想请教您:
1. RNA-Seq的数据也要像其它NGS一样预处理么？比如去掉<50bp的或者也要去掉adaptor（这个我想应该是测序的时候自动做的么）。如果需要预处理都有哪些需要或者软件可以推荐谢谢！
2.网上介绍的比较多的是Tuxedo的流程(tophat/bowtie—cufflinks–cuffdiff)。我也跟着论文上说的试了一下，也有结果。今天看到您这个Trinity流程好像也很牛，我之前也看了一些文献介绍edgeR/DESeq。请问您这个流程和我说的比有什么优势和劣势（除了您这个是不需要reference genome的），或者说edgeR/DESeq比cuffdiff好在哪差在哪。
3.当我们得到了结果后，下游分析都有哪些？比如我知道的是我得到差异表达显著的基因(p值小的)，可以去做GO分析或者pathway分析等。请问还有其它么？因为不管是cuffdiff还是edgeR都有很多输出。
1. RNA-Seq的数据一般都要预处理。因为公司给你的数据中一般含有低质量的reads。同时也要去除PCR duplicates，这样计算出的表达量准确些。 公司给的数据中adaptor一般是去除掉的。
2. Tuxedo流程和trinity的流程是差不多的。前者有参考基因组时使用，后者是没有参考基因组时使用。
3. 其实，我们做转录组分析，主要目的就是找和表型有关的关键基因。这些关键基因就在差异表达基因中了。下一步，就是通过GO和Pathway分析确定一个或几个最可能的关键基因，用于实验验证。
您好，我在做真菌时候，也是用的trinity拼接和rsem计算表达量，为什么有超过一半以上的数据表达量（FPKM列）都为0呢？是multiple mapping的原因，还是什么呢？不知道您有没有遇到过这种情况。谢谢
陈博士，请教一个问题啊，我是OS X 10.9.2的系统，安装trinity的时候出错：
sh-3.2# make
Using gnu compiler for Inchworm and Chrysalis
cd Inchworm && (test -e configure || autoreconf) \
&& ./configure –prefix=`pwd` && /Applications/Xcode.app/Contents/Developer/usr/bin/make install
checking for a BSD-compatible install… /usr/bin/install -c
checking whether build environment is sane… yes
checking for gawk… no
checking for mawk… no
checking for nawk… no
checking for awk… awk
checking whether make sets $(MAKE)… yes
checking for g++… g++
checking for C++ compiler default output file name… a.out
checking whether the C++ compiler works… yes
checking whether we are cross compiling… no
checking for suffix of executables…
checking for suffix of object files… o
checking whether we are using the GNU C++ compiler… yes
checking whether g++ accepts -g… yes
checking for style of include used by make… GNU
checking dependency style of g++… gcc3
checking for library containing cos… none required
configure: creating ./config.status
config.status: creating Makefile
config.status: creating src/Makefile
config.status: creating config.h
config.status: config.h is unchanged
config.status: executing depfiles commands
Making install in src
if g++ -DHAVE_CONFIG_H -I. -I. -I.. -pedantic -fopenmp -Wall -Wextra -Wno-long-long -Wno-deprecated -m64 -g -O2 -MT Fasta_entry.o -MD -MP -MF “.deps/Fasta_entry.Tpo” -c -o Fasta_entry.o Fasta_entry.cpp; \
then mv -f “.deps/Fasta_entry.Tpo” “.deps/Fasta_entry.Po”; else rm -f “.deps/Fasta_entry.Tpo”; exit 1; fi
clang: warning: argument unused during compilation: ‘-fopenmp’
if g++ -DHAVE_CONFIG_H -I. -I. -I.. -pedantic -fopenmp -Wall -Wextra -Wno-long-long -Wno-deprecated -m64 -g -O2 -MT IRKE_run.o -MD -MP -MF “.deps/IRKE_run.Tpo” -c -o IRKE_run.o IRKE_run.cpp; \
then mv -f “.deps/IRKE_run.Tpo” “.deps/IRKE_run.Po”; else rm -f “.deps/IRKE_run.Tpo”; exit 1; fi
clang: warning: argument unused during compilation: ‘-fopenmp’
IRKE_run.cpp:9:10: fatal error: ‘omp.h’ file not found
#include
^
1 error generated.
make[2]: *** [IRKE_run.o] Error 1
make[1]: *** [install-recursive] Error 1
make: *** [inchworm_target] Error 2
不知道你有没有过类似的经验，安装那个版本都出现同样的错误。
福哥，您好！我用trinity组装出来的转录本有10万多，这个在做后续的GO 分析的时候就有很大的麻烦。我想问下我是否应该在比对之前 进行组装的转录本去冗余呢？
如果是 这个区冗余的标准应该怎么设置，这个有什么文献可以参考的么？
可以參考這一篇：http://www.plantphysiol.org/content/early/2014/01/07/pp.113.234864.abstract
您好！我按照您的方法安装rsem没有成功，于是按照http://deweylab.biostat.wisc.edu/rsem/README.html#compilation进行安装，我设置的路径是export PATH=”/home/kanlab/rsem/rsem-1.2.22:$PATH”，苦恼的是也没有成功。linux新手，期待您的帮助。
我经常在www.pls173.com网站买魔兽美服代练，服务态度好，代练基本上是立即就开始了，金币发货快，秒发。下次买嗨来这里。