Trinity tutorial

Posted by Cooper on February 1, 2019

教程:

Nature protoocl

Trinity tutorial github

Trinity software paper

为什么要de novo?:

Genome-guided assembly VS De novo assembly?

1.缺乏高质量参考基因组,

2.RNA-seq 可以注释转录基因和外显子结构

Genomic DNA assembly: one locus-one contig

Transcriptome : one contig per distinct transcript(isoform) rather than per locus;and different transcripts will have different coverage, reflecting their different expression levels.

安装:

1
conda create -n Trinity Trinity -y

使用:

1
 Trinity --seqType fq --max_memory 50G --left reads_1.fq  --right reads_2.fq --CPU 6 --jaccard_clip

参数说明:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
必须的参数:
--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/。

results:

image-20200229210458849

评估组装质量:

1.RNA Seq Read Representation by Trinity Assembly

1
bowtie2-build Trinity.fasta Trinity.fasta
1
bowtie2 -p 30  -q --no-unal -k 20 -x Trinity.fasta -1 ../C1_1_1.fq -2 ../C1_1_2.fq  2>align_stats.txt| samtools view -@30 -Sb -o bowtie2.bam &
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
30251890 reads; of these:
  30251890 (100.00%) were paired; of these:
    6978680 (23.07%) aligned concordantly 0 times
    1325558 (4.38%) aligned concordantly exactly 1 time
    21947652 (72.55%) aligned concordantly >1 times
    ----
    6978680 pairs aligned concordantly 0 times; of these:
      199258 (2.86%) aligned discordantly 1 time
    ----
    6779422 pairs aligned 0 times concordantly or discordantly; of these:
      13558844 mates make up the pairs; of these:
        1267412 (9.35%) aligned 0 times
        243382 (1.80%) aligned exactly 1 time
        12048050 (88.86%) aligned >1 times
97.91% overall alignment rate

比对率大于80%。

2.Full-length transcript analysis for model and non-model organisms using BLAST+

对全长转录本评估可从三方面:

1.align the assembled transcripts to the reference transcripts and examine the length coverage.

2.or non-model organisms, no such reference transcript set is available. If a high quality annotation exists for a closely related organism, then one might compare the assembled transcripts to that closely related transcriptome to examine full-length coverage.

3.In other cases, a more general analysis to perform is to align the assembled transcripts against all known proteins and to determine the number of unique top matching proteins that align across more than X% of its length.

操作:

1.从 ANiseed 下载 C.savignyi 的 Transcript Models:

image-20200301004658103

为什么是protein?

1
服务器文件位置:~/data/clean_data/trinity_out_dir/protein_data
1
2
3
makeblastdb -in Ciona_savignyi.CSAV2.0.81_pep.2018.fa -dbtype prot

blastx -query trinity_out_dir.Trinity.fasta -db protein_data/Ciona_savignyi.CSAV2.0.81_pep.2018.fa -out blastx.outfmt6 -evalue 1e-20 -num_threads 20 -max_target_seqs 1 -outfmt 6
1
analyze_blastPlus_topHit_coverage.pl blastx.outfmt6 trinity_out_dir.Trinity.fasta protein_data/Ciona_savignyi.CSAV2.0.81_pep.2018.fa 

结果:

1
2
3
4
5
6
7
8
9
10
11
#hit_pct_cov_bin	count_in_bin	>bin_below
100	8770	8770
90	815	9585
80	724	10309
70	711	11020
60	741	11761
50	839	12600
40	910	13510
30	1078	14588
20	1317	15905
10	691	16596

BuSCO

参考:https://busco.ezlab.org/busco_userguide.html

简书:https://www.jianshu.com/p/5041460f7a5d

1
busco -m MODE -i INPUT -o OUTPUT -l LINEAGE

image-20200303192341715

image-20200303181803710

image-20200303192709631

下游分析

基于比对

rsem

1
align_and_estimate_abundance.pl --transcripts trinity_out_dir.Trinity.fasta --seqType fq --samples_file ../samples_trinity.txt --est_method RSEM --output_dir transcript_quant --thread_count 54 --trinity_mode --aln_method bowtie2 --prep_reference