标签归档：Trinity

Genome-guided Trinity for Gene Structure Annotation

发表于由
1

使用genome来引导Trinity进行基因结构注释。

RNA-seq的一个主要用途是识别基因组的转录区，重构转录子结构，同时，鉴定转录子的可变剪切。

现在最新的基于genome的转录子预测方法是将RNA-seq的reads使用剪接比对的方法比对到基因组，然后组装比对结果从而得到转录子的结构。(eg. cufflinks, scripture)。我们将这种方法称为：align-reads then assemble-alignments

Trinity可以进行不需要参考基因组的de novo组装，见：Trinity的安装与使用；也能进行有参考基因组支持的组装：即将RNA-Seq比对到genome、RNA-Seq read的de novo组装 和 转录子比对 结合起来。

1. 步骤

1.1 align-reads

使用GSNAP来将reads比对到基因组。将基因组分成各个被reads覆盖的区。

1.2 assemble-reads

对每个区使用Trinity对相应的reads进行组装。

1.3 align-transcripts

使用PASA软件调用GMAP来将Trinity-assembled transcripts比对到genome.

1.4 assemble-transcript_alignments

使用PASA软件来组装上一步骤的比对结果，得出完整的转录子结构，同时，也能解析可变剪接的转录子结构。该步骤和上一步骤其实是在同一个PASA程序中执行得到的。

2. 需要的软件

Trinity
GSNAP & GMAP
PASA

3. 运行

Below, we describe the steps required for running the genome-guided Trinity-based transcript reconstruction pipeline. 适合于真菌物种，其基因密度较大。

3.1 Align RNA-Seq reads to the genome

$ $TRINITY_HOME/util/alignReads.pl --seqType fq --left reads.left.fq --right reads.right.fq --target genome.fasta --aligner gsnap -- -t 8
$ samtools view gsnap_out/gsnap.coordSorted.bam > gsnap.coordSorted.sam

3.2 Assemble the aligned reads using Trinity

$ % $TRINITY_HOME/util/prep_rnaseq_alignments_for_genome_assisted_assembly.pl --SS_lib_type FR --coord_sorted_SAM gsnap.coordSorted.sam -I 1000000
$ find Dir_* -name "*reads" > read_files.list
$ $TRINITY_HOME/util/GG_write_trinity_cmds.pl --reads_list_file read_files.list --paired --SS --jaccard_clip > trinity_GG.cmds
$ $TRINITY_HOME/Inchworm/bin/ParaFly -c trinity_GG.cmds -CPU 6 -failed_cmds trinity_GG.cmds.failed -v
$ find Dir_*  -name "*inity.fasta" -exec cat {} + | $TRINITY_HOME/util/inchworm_accession_incrementer.pl > Trinity_GG.fasta

3.3 Align and assemble the Trinity-reconstructed transcripts using the PASA pipeline

$ cp $PASA_HOME/pasa_conf/pasa.alignAssembly.Template.txt alignAssembly.config
$ perl -p -i -e 's/MYSQLDB=.*/MYSQLDB=sample_mysql_database/' alignAssembly.config
$ $PASA_HOME/scripts/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g genome.fasta -t Trinity_GG.fasta --ALIGNERS blat,gmap --transcribed_is_aligned_orient --stringent_alignment_overlap 30.0

只根据两个转录组数据分析差异表达基因

发表于由
5

1. 数据

手头有两个转录组数据，分别是某一真菌物种(该物种没有基因组数据)的双核菌丝阶段和子实体阶段的转录组数据。测序平台是Illumina Hiseq2000；插入片段长度200bp，测序的reads长度90bp。得到的数据文件为：

/home/user/RNA-seq/mycelium_reads1.fastq
/home/user/RNA-seq/mycelium_reads2.fastq
/home/user/RNA-seq/fruitingbody_reads1.fastq
/home/user/RNA-seq/fruitingbody_reads2.fastq

2. 数据的预处理

2.1. 去除N的比例大于5%的reads;去除低质量reads（质量值Q≤20的碱基数占整个read的50％以上）；
2.2. 根据FastQC对上述过滤后的reads的质量检测，去除reads首尾各10bp的碱基，得到的预处理数据为：

/home/user/RNA-seq/clean_reads/mycelium_reads1.fastq
/home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

3. 使用Trinity和TGICL进行转录组的组装

3.1. 对mycelium数据和fruitingbody数据分别进行转录组的组装

$ pwd
/home/user/RNA-seq/

$ mkdir assembly
$ cd assembly

$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output mycelium_contig --left /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq --right /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output fruitingbody_contig --left /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq --right /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

当然，需要将生成的两个转录组的Trinity.fasta序列按长度进行排序；统一更改的fasta
头名称；更改fasta文件名

$ cp mycelium_contig/Trinity.fasta mycelium_contigs.fasta
$ cp fruitingbody_contig/Trinity.fasta fruitingbody_contigs.fasta

3.2. 使用TGICL将两个转录组序列合并

$ pwd
/home/user/RNA-seq/assembly

$ mkdir all_tissue_contig
$ cd all_tissue_contig

$ cat ../mycelium_contigs.fasta ../fruitingbody_contigs.fasta > all.contigs.fasta
$ tgicl -F all.contigs.fasta

tgicl生成了两个有用的文件： asm_1/contigs  和 all.contigs.fasta.single
tons。其中前者是聚类后的contigs结果；后者是没有聚类的单独的contigs的序列名，需
要分别到../mycelium_contigs.fasta 和 ../fruitingbody_contigs.fasta文
件中提取出相应的序列: all.contigs.fasta.singletons.mycelium.fasta 和 
all.contigs.fasta.singletons.fruitingbody.fasta

$ cat asm_1/contigs all.contigs.fasta.singletons.mycelium.fasta all.contigs.fasta.singletons.fruitingbody.fasta > all_contigs.fasta

在对all_contigs.fasta进行fasta头的重命名，并将序列按长度排序；同时要得到all_
contigs.fasta中的序列和mycelium_contigs.fasta，fruitingbody_contigs.
fasta中序列的对应关系。

$ cp all_contigs.fasta ../

3.3. 至此得出转录组的组装结果：

/home/user/RNA-seq/assembly/all_contigs.fasta

4. 使用cufflinks来分析差异表达基因

4.1 使用tophat来将预处理后的reads比对到转录组序列上

$ pwd
/home/user/RNA-seq/

$ mkdir cufflinks
$ cd cufflinks

$ bowtie2-build ../all_contigs.fasta all_contigs

$ tophat -o tophat_out_mycelium -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ tophat -o tophat_out_fruitingbody -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

4.2 自制一个transcripts.gtf文件
由于cuffdiff运行需要一个参考序列的transcripts.gtf文件: 该文件有9列，使用tab分隔；使exon的范围为整个contig。其格式如：

AA.auricula_all_contig_1     chenlianfu  exon  1  9401  .  .  .  gene_id "A.auricula_all_contig_1"; transcript_id "A.auricula_all_contig_1";
A.auricula_all_contig_10     chenlianfu  exon  1  4464  .  .  .  gene_id "A.auricula_all_contig_10"; transcript_id "A.auricula_all_contig_10";
A.auricula_all_contig_100    chenlianfu  exon  1  3090  .  .  .  gene_id "A.auricula_all_contig_100"; transcript_id "A.auricula_all_contig_100";
A.auricula_all_contig_1000   chenlianfu  exon  1  1768  .  .  .  gene_id "A.auricula_all_contig_1000"; transcript_id "A.auricula_all_contig_1000";
A.auricula_all_contig_10000  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10000"; transcript_id "A.auricula_all_contig_10000"; 
A.auricula_all_contig_10001  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10001"; transcript_id "A.auricula_all_contig_10001";

4.3 使用cuffdiff来分析转录子的表达量和差异表达基因

$ pwd
/home/user/RNA-seq/cufflinks

$ cuffdiff -L mycelium,fruitingbody --library-type fr-unstranded -p 8 -o cuffdiff ./transcriptome.gtf ./tophat_out_mycelium/accepted_hits.bam ./tophat_out_fruitingbody/accepted_hits.bam

至此得出转录子的表达量数据和差异表达分析

/home/user/RNA-seq/cufflinks/cuffdiff/gene_exp.diff

Trinity的安装与使用

发表于由
15

一、 Trinity简介

Trinity，是由 the Broad Institute 开发的转录组de novo组装软件，由三个独立的软件模块组成：Inchworm,Chrysalis和Butterfly。三个软件依次来处理大规模的RNA-seq的reads数据。Trinity的简要工作流程为：
Inchworm: 将RNA-seq的原始reads数据组装成Unique序列；
Chrysalis: 将上一步生成的contigs聚类，然后对每个类构建Bruijn图；
Butterfly: 处理这些Bruijn图，依据图中reads和成对的reads来寻找路径，从而得到具有可变剪接的全长转录子，同时将旁系同源基因的转录子分开。
Trinity发表在 Nature Biotechnology。

二、 Trinity的安装

1. 下载Trinity。
2. 安装Trinity。

$ tar zxvf trinityrnaseq_r2013-02-25.tgz
$ cd trinityrnaseq_r2013-02-25
$ make
仅需要在安装目录下进行make即可。该命令编译了由C++编写的Inchworm和Chrysalis，而
使用Java编写的Butterfly则不需要编译，可以直接使用。

三、Trinity的使用

1. 直接运行安装目录下的程序Trinity.pl来使用该软件，不带参数则给出使用帮助。其典型用法为：

Trinity.pl --seqType fq --JM 50G --left reads_1.fq  --right reads_2.fq --CPU 6

2. Trinity参数

必须的参数：
--seqType     reads的类型：(cfa, cfq, fa, or fq)
--JM          jellyfish使用多少G内存用来进行k-mer的计算，包含‘G’这个字符
--left        左边的reads的文件名
--rigth       右边的reads的文件名
--single      不成对的reads的文件名

可选参数：

Misc：
--SS_lib_type        reads的方向。成对的reads: RF or FR; 不成对的reads
: F or R。在数据具有链特异性的时候，设置此参数，则正义和反义转录子能得到区分。默认
情况下，不设置此参数，reads被当作非链特异性处理。FR: 匹配时，read1在5'端上游, 
和前导链一致, read2在3'下游, 和前导链反向互补. 或者read2在上游, read1在下游反
向互补; RF: read1在5'端上游, 和前导链反向互补, read2在3'端下游, 和前导链一致;
--output             输出结果文件夹。默认情况下生成trinity_out_dir文件夹并
将输出结果保存到此文件夹中。
--CPU                使用的CPU线程数，默认为2
--min_contig_length  报告出的最短的contig长度。默认为200
--jaccard_clip       如果两个转录子之间有UTR区重叠，则这两个转录子很有可能在
de novo组装的时候被拼接成一条序列，称为融合转录子(Fusion Transcript)。如果有
fastq格式的paired reads，并尽可能减少此类组装错误，则选用此参数。值得说明的是：
1. 适合于基因在基因组比较稠密，转录子经常在UTR区域重叠的物种，比如真菌基因组。而对
于脊椎动物和植物，则不推荐使用此参数; 2. 要求fastq格式的paired reads文件(文件
中reads名分别以/1和/2结尾，以利于软件识别)，同时还需要安装bowtie软件用于reads
的比对; 3. 单独使用具有链特异性的RNA-seq数据的时候，能极大地减少UTR重叠区很小的
融合转录子; 4. 此选项耗费运算，若没必要，则不用此参数。
--prep               仅仅准备一些文件(利于I/O）并在kmer计算前停止程序运行
--no_cleanup         保留所有的中间输入文件
--full_cleanup       仅保留Trinity fasta文件，并重命名成${output_dir}.
Trinity.fasta
--cite               显示Trinity文献引证和一些参与的软件工具
--version            报告Trinity版本并推出

Inchworm 和 K-mer 计算相关选项：
--min_kmer_cov      使用Inchworm来计算K-mer数量时候，设置的Kmer的最小值。
默认为1
--inchworm_cpu      Inchworm使用的CPU线程数，默认为6和--CPU设置的值中的
小值。

Chrysalis相关选项：
--max_reads_per_graph   在一个Bruijn图中锚定的最大的reads数目，默认为200
000
--no_run_chrysalis      运行Inchworm完毕，在运行chrysalis之前停止运行
Trinity
--no_run_quantifygraph  在平行化运算quantifygrahp前停止运行Trinity

Butterfly相关选项：
--bfly_opts                    Butterfly额外的参数
--max_number_of_paths_per_node 从node A -> B,最多允许多少条路径。默认
为10
--group_pairs_distance         最大插入片读长度，默认为500
--path_reinforcement_distance  延长转录子路径时候，reads间最小的重叠碱基
数。默认PE:75; SE：25
--no_triplet_lock              不锁定triplet-supported nodes
--bflyHeapSpaceMax             运行Butterfly时java最大的堆积空间，默认
为20G
--bflyHeapSpaceInit            java初始的堆积空间，默认为1G
--bflyGCThreads                java进行无用信息的整理时使用的线程数，默
认由java来决定
--bflyCPU                      运行Butterfly时使用的CPU线程数，默认为2
--bflyCalculateCPU             计算Butterfly所运行的CPU线程数，由公式
 80% * max_memory / maxbflyHeapSpaceMax 得到
--no_run_butterfly             在Chrysalis运行完毕后，停止运行Butterfly

Grid-computing选项：
--grid_computing_module   选定Perl模块，在/Users/bhaas/SVN/trinityr
naseq/trunk/PerlLibAdaptors/。

3. 适合于illumina测序数据的真菌物种转录组组装的Trinity命令为：

Trinity.pl --seqType fq --JM 50G --left reads_1.fq  --right reads_2.fq --SS_lib_type FR --output transcriptome_tissue --CPU 24 --jaccard_clip --inchworm_cpu 24 --group_pairs_distance 500 --bflyCPU 24

4. Trinity生成的结果文件
运行程序结束后，转录组结果为trinity_out_dir/Trinity.fasta。可以使用软件所带的一支程序分析转录组统计信息。

$ $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta
Total trinity transcripts:	30706
Total trinity components:	26628
Contig N50: 554

三. Trinity运行原理与过程

1. 检测java的可运行性，因为buttfly会用到

2. 运行jellyfish，使用其dump命令得到jellyfish.kmers.fa文件

3. Inchworm(Linear contig construction from k-mers)

assembles the RNA-seq data into the unique sequences of transcripts, often generating full-length transcripts for a dominant isoform, but then reports just the unique portions of alternatively spliced transcripts.

4. Chrysalis

clusters the Inchworm contigs into clusters and constructs complete de Bruijn graphs for each cluster. Each cluster represents the full transcriptonal complexity for a given gene (or sets of genes that share sequences in common). Chrysalis then partitions the full read set among these disjoint graphs.

5. Butterfly

then processes the individual graphs in parallel, tracing the paths that reads and pairs of reads take within the graph, ultimately reporting full-length transcripts for alternatively spliced isoforms, and teasing apart transcripts that corresponds to paralogous genes.

四. 注意事项

3.1 Trinity分步运行

当数据量比较大的时候，trinity运行的时间会很长，同时，内存不够等情况出现的时候有可能程序运行崩溃。最好是分步运行。下一步会接着前一步进行下去。

Stage 1: generate the kmer-catalog and run Inchworm: –no_run_chrysalis

Stage 2: Chrysalis clustering of inchworm contigs and mapping reads: –no_run_quantifygraph

Stage 3: Chrysalis deBruijn graph construction: –no_run_butterfly

Stage 4: Run butterfly, generate final Trinity.fasta file. (exclude –no_ options)

3.2 计算资源

Ideally, you will have access to a large-memory server, ideally having ~1G of RAM per 1M reads to be assembled (but often, much less memory may be required).

The assembly from start to finish can take anywhere from ~1/2 hour to 1 hour per million reads (your mileage may vary). 个人记录了一次，使用dell服务器，64GB RAM，24 threads : 53M 的reads，运行了16.5h（平均3.2M/h），内存使用峰值为43G.