Trinity：从RNA-seq数据中重建转录组的高效工具
Trinity：从RNA-seq数据中重建转录组的高效工具
Trinity是一种从RNA-seq数据中有效、稳健地从头重建转录组的新方法。它由三个独立的软件模块组成：Inchworm、Chrysalis和Butterfly，分别对应蝴蝶的三个生长阶段。Trinity能够处理大量RNA-seq读数，将序列数据分割成多个de Bruijn图，每个图表示给定基因或基因座处的转录复杂性，然后独立处理每个图以提取全长剪接同种型并梳理来自类似基因的转录物。
Trinity的工作流程
Trinity的工作流程可以分为三个主要阶段：
虫阶段（Inchworm）：将RNA-seq数据组装成独特的转录本序列，通常为显性同种型产生全长转录本，但随后仅报告可变剪接转录本的独特部分。
蛹阶段（Chrysalis）：将虫阶段的产物contigs聚集成clusters，并为每个cluster构建完整的de Bruijn图。每个簇代表给定基因（或共享序列的一组基因）的完整转录复杂性。然后在这些不相交的图中划分full read set。
蝶阶段（Butterfly）：平行处理各个图，跟踪图中reads和pairs of reads的路径，最终报告可变剪接isoforms的全长转录物，并梳理对应于类似基因的转录物。
Trinity的硬件和配置要求
- Inchworm和Chrysalis阶段是内存密集型的。一个基本的建议是每1M对Illumina reads大约有1G的RAM。较简单的转录组（低等真核生物）比复杂的转录组（如脊椎动物）需要更少的内存。
- Trinity可能需要数百GB的可用磁盘空间，并且可以在运行期间生成数千个中间文件。但是，生成的最终输出文件很少，而且通常相对较小（MB而不是许多GB）。
- 最好有一个临时工作空间，其中有足够的磁盘空间可以在作业执行期间使用。
- 如果不能直接访问高内存机器（通常有256G或512G RAM），考虑在一个外部可用资源上运行Trinity。
Trinity的运行命令
一个典型的用于组装非链特异性RNA-seq数据的Trinity命令将类似于这样，在单个高内存服务器上运行整个过程：
Trinity --seqType fq --max_memory 50G \
--left reads_1.fq.gz --right reads_2.fq.gz --CPU 6
如果你有多组fastq文件，比如对应于多种组织类型或条件等，你可以像下面这样向Trinity指示它们：
Trinity --seqType fq --max_memory 50G \
--left condA_1.fq.gz,condB_1.fq.gz,condC_1.fq.gz \
--right condA_2.fq.gz,condB_2.fq.gz,condC_2.fq.gz \
--CPU 6
或者更好的方法是，创建一个“samples.txt”文件，像下面这样描述数据：
# --samples_file <string> tab-delimited text file indicating biological replicate relationships.
# ex.
# cond_A cond_A_rep1 A_rep1_left.fq A_rep1_right.fq
# cond_A cond_A_rep2 A_rep2_left.fq A_rep2_right.fq
# cond_B cond_B_rep1 B_rep1_left.fq B_rep1_right.fq
# cond_B cond_B_rep2 B_rep2_left.fq B_rep2_right.fq
这个相同的样本文件随后可以与其他下游分析步骤一起使用，包括表达量化和差异表达分析。
Trinity的输出
当Trinity完成后，它将创建一个trinity_out_dir.Trinity.fasta输出文件（或者基于您指定的输出目录的前缀）。Trinity基于共享序列内容将转录物分组成簇。这样的转录本簇被非常宽泛地称为“基因”。这个信息被编码在Trinity的accession中。其中一个转录本的Fasta条目的格式如下：
>TRINITY_DN1000_c115_g5_i1 len=247 path=[31015:0-148 23018:149-246]
AATCTTTTTTGGTATTGGCAGTACTGTGCTCTGGGTAGTGATTAGGGCAAAAGAAGACAC
ACAATAAAGAACCAGGTGTTAGACGTCAGCAAGTCAAGGCCTTGGTTCTCAGCAGACAGA
AGACAGCCCTTCTCAATCCTCATCCCTTCCCTGAACAGACATGTCTTCTGCAAGCTTCTC
CAAGTCAGTTGTTCACAGGAACATCATCAGAATAAATTTGAAATTATGATTAGTATCTGA
TAAAGCA
accession编码Trinity的“gene”和“isoform”信息。在上面的例子中，accessionTRINITY_DN1000_c115_g5_i1表示TRINITY读取簇TRINITY_DN1000_c115，geneg5和isoformi1。因为一个给定的Trinity运行涉及许多读集群，每个读集群都是单独组装的，并且由于“基因”编号在给定的处理读集群中是唯一的，“基因”标识符应该被认为是读取簇和相应基因标识符的集合，在这种情况下应该是TRINITY _ DN1000 _ c115 _ g5。
因此，总之，上面的例子对应于编码“isoform id:TRINITY_DN1000_c115_g5_i1的“基因 id:TRINITY_DN1000_c115_g5。
存储在标头中的Path信息（“Path = [31015:0-14823018:149-246]”）指示在Trinity压缩的de Bruijn图中遍历以构造该转录本的路径。在这种情况下，节点“31015”对应于转录本的序列范围0-148，节点23018对应于转录本序列的序列范围149-246。只有在给定的Trinity基因标识符的上下文中，节点数才是唯一的，因此图形节点可以在同种型之间进行比较，以识别给定基因的每个同种型的唯一和共享序列。
Trinity的基因组引导组装
如果一个基因组序列是可用的，Trinity提供了一种方法，即reads首先比对到基因组，根据基因座进行分区，然后在每个基因座从头转录组装配。在这个用例中，基因组仅被用作将重叠读数分组成簇的底物，然后将其分别输入Trinity进行从头转录组装配。这与典型的基因组引导方法（例如cufflinks）非常不同，其中比对的读数被缝合到转录物结构中，并且基于参考基因组序列重建转录物序列。在这里，基于实际的读取序列重新构建文本。
为什么要这么做？你可能有一个参考基因组，但是你的样本可能来自一个基因组与参考基因组不完全匹配的有机体。基因组引导的从头组装应该捕获RNA-Seq样本中包含的序列变异，其形式是从头重建的转录本。与无基因组从头组装相比，它也可以帮助在情况下，你有旁系同源或其他基因与共享序列，因为基因组是用来分区做任何从头组装之前reads的locus。如果你有一个高度片段化的基因组草图，那么你可能最好执行一个无基因组的从头转录组装配。
用户必须以coordinate-sorted bam文件的形式提供与Trinity的reads齐。使用GSNAP、TopHat、STAR或其他常用的RNA-Seq读对齐工具生成bam文件，并确保它是通过运行“samtools sort”进行坐标排序的。
要运行基因组引导的Trinity并让Trinity执行GSNAP来对齐读数，运行Trinity如下：
Trinity --genome_guided_bam rnaseq.coordSorted.bam \
--genome_guided_max_intron 10000 \
--max_memory 10G --CPU 10
当然，使用最大的内含子长度，最有意义的给你的目标有机体。
如果需要读数的质量修剪，则应在将读数与基因组对齐之前执行，因为Trinity将只使用提供给它的坐标排序的bam文件中存在的读数。
Trinity的输出统计信息
$TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes': 377
Total trinity transcripts: 384
Percent GC: 38.66
########################################
Stats based on ALL transcript contigs:
########################################
Contig N10: 3373
Contig N20: 2605
Contig N30: 2219
Contig N40: 1936
Contig N50: 1703
Median contig length: 772
Average contig: 1047.80
Total assembled bases: 402355
#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################
Contig N10: 3373
Contig N20: 2605
Contig N30: 2216
Contig N40: 1936
Contig N50: 1695
Median contig length: 772
Average contig: 1041.98
Total assembled bases: 392826
转录本丰度估计
为了估计Trinity-reconstructed transcripts的表达水平，我们使用RSEM软件支持的策略。我们首先将原始的rna-seq读取与Trinity转录物对齐，然后运行RSEM以估计映射到每个contig的rna-seq片段的数量。由于各个转录本的丰度在样品之间可能有显著差异，因此必须分别检查每个样品的读数，以获得样品特异性丰度值。
对于alignments，我们使用“bowtie”而不是“tophat”。这有两个原因。首先，因为我们将读数映射到重建的cDNA而不是基因组序列，所以正确对齐的读数不需要在内含子之间留空隙。其次，RSEM软件目前只兼容无间隙对齐。
分别量化每个样本的转录表达
下面的脚本将运行RSEM，它首先使用Bowtie对齐器将RNA-Seq读数与Trinity转录物对齐，然后执行丰度估计。这个过程是
${TRINITY_HOME}/util/align_and_estimate_abundance.pl --seqType fq \
--left RNASEQ_data/Sp_ds.left.fq.gz --right RNASEQ_data/Sp_ds.right.fq.gz \
--transcripts trinity_out_dir/Trinity.fasta \
--output_prefix Sp_ds --est_method RSEM --aln_method bowtie \
--trinity_mode --prep_reference --output_dir Sp_ds.RSEM
完成后，RSEM将生成两个文件‘Sp_ds.isoforms.results’and‘Sp_ds.genes.results’。这些文件包含Trinity transcript和component（Trinity类似于Isoform和gene）rna-seq片段计数和归一化表达值。
> head Sp_ds.RSEM/Sp_ds.isoforms.results
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct
TRINITY_DN0_c0_g1_i1 TRINITY_DN0_c0_g1 1894 1629.48 62.00 432.83 401.52 100.00
TRINITY_DN0_c1_g1_i1 TRINITY_DN0_c1_g1 271 28.49 7.00 2795.30 2593.09 100.00
TRINITY_DN100_c0_g1_i1 TRINITY_DN100_c0_g1 320 62.60 2.00 363.44 337.15 100.00
TRINITY_DN100_c1_g1_i1 TRINITY_DN100_c1_g1 948 683.48 18.00 299.59 277.91 100.00
TRINITY_DN100_c2_g1_i1 TRINITY_DN100_c2_g1 789 524.48 17.00 368.72 342.05 100.00
TRINITY_DN101_c0_g1_i1 TRINITY_DN101_c0_g1 956 691.48 36.00 592.24 549.40 100.00
TRINITY_DN101_c1_g1_i1 TRINITY_DN101_c1_g1 279 33.31 4.00 1365.90 1267.09 100.00
TRINITY_DN101_c2_g1_i1 TRINITY_DN101_c2_g1 250 17.44 0.00 0.00 0.00 0.00
TRINITY_DN102_c0_g1_i1 TRINITY_DN102_c0_g1 4736 4471.48 369.00 938.75 870.84 100.00
生成一个计数矩阵并执行交叉样本标准化
${TRINITY_HOME}/util/abundance_estimates_to_matrix.pl --est_method RSEM \
--out_prefix Trinity_trans \
Sp_ds.RSEM/Sp_ds.isoforms.results \
Sp_hs.RSEM/Sp_hs.isoforms.results \
Sp_log.RSEM/Sp_log.isoforms.results \
Sp_plat.RSEM/Sp_plat.isoforms.results
这应该找到一个名为'Trinity_trans.counts.matrix'的文件，其中包含映射到每个转录本的RNA-Seq片段的计数。此外，还会生成一个包含归一化表达式值的矩阵，称为'Trinity_trans.TMM.EXPR.matrix'。这些表达量的值根据测序深度和转录本长度进行非正规化，然后在大多数转录本没有差异表达的假设下通过TMM标准化进行缩放。
差异表达分析
为了检测差异表达的转录本，使用我们的计数矩阵运行Bioiconductor软件包edgeR：
${TRINITY_HOME}/Analysis/DifferentialExpression/run_DE_analysis.pl \
--matrix Trinity_trans.counts.matrix \
--method edgeR \
--dispersion 0.1 \
--output edgeR