Trinity

最后发布时间 : 2024-06-10 19:12:09 浏览量 :

学习资料

PMC3875132

Trinity是一种从 RNA-seq 数据有效、稳健的从头重建转录组（ de novo reconstruction of transcriptomes）的新方法。Trinity包括三个独立的软件模块。Inchworm(虫), Chrysalis(蛹), and Butterfly(蝶),形象的表示为蝴蝶的三个阶段，依次应用于处理大量 RNA-seq 读数。Trinity 将序列数据分割成许多单独的 de Bruijn 图，每个图表示给定基因或基因座处的转录复杂性，然后独立处理每个图以提取全长剪接同种型(full-length splicing isoforms)并梳理来自类似基因的转录物。简而言之，整个过程是这样的:

虫阶段将 RNA-seq 数据组装成独特的转录本序列(unique sequences of transcripts)，通常为显性同种型( dominant isoform)产生全长转录本，但随后仅报告可变剪接(alternatively spliced)转录本的独特部分。
蛹阶段将虫阶段的产物contigs聚集成clusters，并为每个cluster构建完整的 de Bruijn 图。每个簇代表给定基因(或共享序列的一组基因)的完整转录复杂性( full transcriptonal complexity)。然后蛹在这些不相交的图中划分 full read set。
然后，蝶阶段平行处理各个图，跟踪图中reads和 pairs of reads的路径，最终报告可变剪(alternatively spliced)isoforms的全长转录物( full-length transcripts)，并梳理对应于类似基因的转录物。

Trinity 发表在Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data,转录组组装和下游分析方案已经发表De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity

Trinity --seqType fq --left reads_1.fq --right reads_2.fq --CPU 6 --max_memory 20G

组装的转录本在trinity_out_dir/Trinity.fasta

Trinity Hardware and Configuration Requirements

Inchworm(虫)和Chrysalis(蛹)阶段是内存密集型的。一个基本的建议是每1M 对 Illumina reads大约有1G 的 RAM。较简单的转录组(低等真核生物)比较复杂的转录组(如脊椎动物)需要更少的内存。Trinity 还可能需要数百 GB 的可用磁盘空间，并且可以在运行期间生成数千个中间文件。但是，生成的最终输出文件很少，而且通常相对较小(MB 而不是许多 GB)。最好有一个临时工作空间，其中有足够的磁盘空间可以在作业执行期间使用。

如果您能够在单个高内存多核服务器上运行整个 Trinity 进程，请通过--CPU参数指示要并行运行的butterfly进程的数量，在集群LSF, SGE, UGE, SLURM上处理见这里

在当前的实现中，每100万对reads需要大约半小时到一小时的整个过程。如果不能直接访问高内存机器(通常有256G 或512G RAM)，考虑在一个外部可用资源上运行 Trinity。
请注意，仅 Inchworm 将在6个线程内部设置上限，因为在该设置之外，此步骤的性能不会得到提高)

Running Trinity

一个典型的用于组装非链特异性 RNA-seq 数据的 Trinity 命令将类似于这样，在单个高内存服务器上运行整个过程(目标为每 ~ 1M ~ 76碱基 Illumina 配对读取 ~ 1G RAM，但通常需要更少的内存) :

 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

这个相同的样本文件随后可以与其他下游分析步骤一起使用，包括表达量化和差异表达分析。

请注意，fastq 文件可以像上面所示的那样被 gzip 压缩，在这种情况下，它们应该需要一个’. gz’扩展名。

Options to Consider when Running Trinity

如果转录物来源于重叠 UTR 常见的紧凑基因组，则提供了通过分析整个转录物长度的读取配对的一致性来减轻错误的端到端融合转录物的组装的选项。
如果有特别大的 RNA-Seq 数据集，涉及数亿到数十亿的读数，那么在计算机上标准化是必要的。

#  --no_normalize_reads            :Do *not* run in silico normalization of reads. Defaults to max. read coverage of 200.
#                                       see '--normalize_max_read_cov' under full usage info for tailored settings.
#                                       (note, as of Sept 21, 2016, normalization is on by default)
#
################################################################################ 
####  In silico Read Normalization Options ###
#
#  --normalize_max_read_cov <int>       defaults to 200
#  --normalize_by_read_set              run normalization separate for each pair of fastq files,
#                                       then one final normalization that combines the individual normalized reads.
#                                       Consider using this if RAM limitations are a consideration.
#

尽量减少基因密集型基因组中错误融合的转录本(使用--jaccard_clip)

如果你的转录组 RNA-seq 数据来源于基因密集的致密基因组( a gene-dense compact genome)，例如来自真菌基因组，其中转录物可能经常在 UTR 区域重叠，那么如果您具有配对读数，则可以通过利用--jaccard_clip选项来最小化融合转录物。Trinity 将在几乎没有读配对支持的位置检查读配对和片段转录本的一致性。在脊椎动物和植物的广泛的基因组(expansive genomes)中，这是不必要的，也不推荐。在致密的真菌基因组中，它是强烈推荐的。作为这个分析的一部分，reads与尺蠖重叠(Inchworm contigs)对齐，在尺蠖重叠( Inchworm contigs)中检查reads配对，并且在低配对支持位置剪切重叠。这些剪辑尺蠖重叠然后饲料到蝶蛹下游加工。

将 Trinity 应用到计算网格中，以便对并行步骤进行大规模并行处理机处理

Trinity 有许多并行组件，所有这些组件都可以从单个服务器上的多个 CPU 中受益，但也有例子，比如蝶蛹和蝴蝶，它们可以并行执行数万到数十万个命令，每个命令都有独立的输入和输出。这些原生并行的命令可以在计算场的上下文中进行最有效的计算，将每个命令(或它们的批次)提交给计算网格上的各个节点。

要在计算网格上启动并行作业，您(简单地?)需要有一种方法来启动一组命令到您的计算平台上。这应该是一个脚本，您有(或编写) ，可以采取一个单一的参数: 一个文件，其中包含的命令列表执行。它还应该确保所有命令都成功完成。通过--grid_exec参数指定脚本的路径。这在 Trinity 的用法(在 Trinity--show_full_usage_info下)中描述如下:

Genome-guided Trinity De novo Transcriptome Assembly

如果一个基因组序列是可用的，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 文件中存在的读数。

Output of Trinity Assembly

当 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

Explaining the identifiers: Genes vs. Transcripts

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基因标识符的上下文中，节点数才是唯一的，因此图形节点可以在同种型之间进行比较，以识别给定基因的每个同种型的唯一和共享序列。

De novo assembly of reads using Trinity

为了生成一个reference assembly，我们以后可以用它来分析差异表达式，我们将把不同条件的读数据集合合并到一个针对 Trinity assembly的单一目标中。为此，我们向 Trinity 提供了所有左右 fastq 文件的列表，这些文件分别位于--left和--right参数中，使用逗号分隔，如下所示:

Trinity --seqType fq --SS_lib_type RF  \
           --left 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 \
           --right 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 \
           --CPU 2 --max_memory 1G

在这个数据集上运行 Trinity 可能需要10到15分钟。你会看到它在不同阶段的进展，从Jellyfish 生成 k-mer 目录开始，然后是尺蠖，蝶蛹，最后是蝴蝶。运行一个典型的 Trinity 作业需要每100万个 PE 读取大约1小时和1G 内存。你通常会在一台高内存的机器上运行它，然后让它搅拌几个小时或几天。

运行结束后组装的转录本文件在trinity_out_dir/Trinity.fasta

> head trinity_out_dir/Trinity.fasta
>TRINITY_DN0_c0_g1_i1 len=1894 path=[1872:0-1893] [-1, 1872, -2]
GTATACTGAGGTTTATTGCCTGTAACGGGCAAACTCGAGCAGTTCAATCACGAGGAGATT
ACCAGAAAACACTCGCTATTGCTTTGAAAAAGTTTAGCCTTGAAGATGCTTCAAAATTCA
TTGTATGCGTTTCACAGAGTAGTCGAATTAAGCTGATTACTGAAGAAGAATTTAAACAAA
TTTGTTTTAATTCATCTTCACCGGAACGCGACAGGTTAATTATTGTGCCAAAAGAAAAGC
CTTGTCCATCGTTTGAAGACCTCCGCCGTTCTTGGGAGATTGAATTGGCTCAACCGGCAG
CATTATCATCACAGTCCTCCCTTTCTCCTAAACTTTCCTCTGTTCTTCCCACGAGCACTC
AGAAACGAAGTGTCCGCTCAAATAATGCGAAACCATTTGAATCCTACCAGCGACCTCCTA
GCGAGCTTATTAATTCTAGAATTTCCGATTTTTTCCCCGATCATCAACCAAAGCTACTGG
AAAAAACAATATCAAACTCTCTTCGTAGAAACCTTAGCATACGTACGTCTCAAGGGCACA

Trinity assembly的基础统计信息

 $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

Abundance estimation using RSEM

为了估计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 标准化进行缩放。

Differential Expression Using EdgeR

为了检测差异表达的转录本，使用我们的计数矩阵运行 Bioiconductor 软件包 edgeR:

 ${TRINITY_HOME}/Analysis/DifferentialExpression/run_DE_analysis.pl \
      --matrix Trinity_trans.counts.matrix \
      --method edgeR \
      --dispersion 0.1 \
      --output edgeR

stringTie组装和量化转录本组装结果功能注释