无参转录组序列组装 — Trinity

2019-07-17 本文已影响8人六六_ryx

转录组分析策略

转录组分析策略根据有无参考转录组可以分为两类：

有参
- 比对—定量—差异分析—功能富集分析等
无参
- 组装—定量—差异分析—功能富集分析等
image

无参考转录组序列组装—Trinity概述

其中Trinity 是无参考转录组从头组装转录组的常用软件，且trinity的使用文档非常详细，整合的内容非常完整，包括从组装，比对，定量到差异分析等。因此有大神也推荐Trinity可作为初学者了解熟悉转录组分析流程的入门和进阶学习文档。

官方文档： https://github.com/trinityrnaseq/trinityrnaseq/wiki

Trinity通过有秩序的对大规模的RNA-seq Reads数据进行读取，高效的完成转录组的组装，包含三个独立的软件模块：

Inchworm （虫）（C++）
Chrysalis （蛹）（C++）
Butterfly （蝶）（Java）

image

Trinity组装原理

Trinity组装依据的算法是de Bruijn Graph,即从打断的文库中提取一定长度的K-mer，然后根据k-1错位相似的方法拼接组装的可能路径，最终确定完整的参考组装转录组。

Trinity根据该原理，将主要操作步骤分为3个模块，分别形象的命名为虫，蛹，蝶：

序列延伸 (inchworm) ——虫
- 将 reads切为 k-mers (k bp长度的短片段)
- 利用Overlap关系对k-mers进行延伸 (贪婪算法)
- 输出所有的序列 (“contigs”)
构建 de Bruijn graph (chrysalis)——蛹
- 聚类所有相似区域大于k-1bp的 contigs
- 构图 (区分不同的 “components”)
- 将reads比对回 components，进行验证
解图，列举转录本 (butterfly)——蝶
- 拆分graph 为线性序列
- 使用reads以及 pairs关系消除错误序列

Trinity 组装流程

Trinity组装流程主要包括以下几个步骤：

组装前数据的质控
组装
组装质量的评估
下游分析
- 差异分析
- 富集分析
- 功能分析

STEP0: 准备工作

所有培训数据路径

cd /home/tech/NGS-example/
ll

Trinity数据路径

cd /home/tech/NGS-example/Trinity-example
ll

cp -r /home/tech/NGS-example/Trinity-example/RNASEQ_data/ ./
$ls RNASEQ_data/
Sp_ds.left.fq.gz   Sp_hs.left.fq.gz   Sp_log.left.fq.gz   Sp_plat.left.fq.gz
Sp_ds.right.fq.gz  Sp_hs.right.fq.gz  Sp_log.right.fq.gz  Sp_plat.right.fq.gz

STEP1: 组装前质控

fatqc——trimmamatic等去除低质量的reads等

STEP2: 组装

Trinity 基本命令的含义：

--seqType : [fa/fq] 原始数据的格式，可以是fasta或fastq，可以是.gz或不是
--left/--right：针对双端测序，如果是同一样本不同时期的多个文件，可以中间以逗号分隔，没有空格。
--single: 针对单端测序
--SS_lib_type：[RF]—双端链特异测序；
--output：输出文件，命名必须带Trinity
--CPU：组装使用CPU个数
--max_memory：组装所需内存（eg：10G）

## my_trinity_script.sh
#!/bin/bash
/software/trinityrnaseq-2.2.0/Trinity --seqType fq \
--left RNASEQ_data/Sp_ds.left.fq.gz,RNASEQ_data/Sp_hs.left.fq.gz,RNASEQ_data/Sp_log.left.fq.gz,RNASEQ_data/Sp_plat.left.fq.gz \
--right RNASEQ_data/Sp_ds.right.fq.gz,RNASEQ_data/Sp_hs.right.fq.gz,RNASEQ_data/Sp_log.right.fq.gz,RNASEQ_data/Sp_plat.right.fq.gz \
--CPU  1 \
--max_memory 1G

nohup my_trinity_script.sh>& my_trinity_script.log &

trinity组装结果

会自动生成Trinity_out_dir, 其中Trinity.fasta即组装的初步结果

image

组装结果统计

/software/trinityrnaseq-2.2.0/util/TrinityStats.pl trinity_out_dir/Trinity.fasta > ./trinity_out_dir/assembly_report.txt
![image](https://img.haomeiwen.com/i8242255/23e540035140a649.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

STEP3: 组装后质量评估

组装质量评估标准

Unigene数量
N50
比对率--比对率大于80%
注释比率--物种近缘性良好CDS序列相对完整60%以上
核心蛋白比对率--真核生物中存在一些高度保守区域所编码的蛋白（2748/80%以上）

组装完整性

N50长度，可以初步评估，但不是越长越好，对应物种考虑
通过统计Unigene对近缘物种基因覆盖度分布

组装准确性

组装长度会影响定量的准确性

后续定量准确性

组装冗余度

影响定量准确性的最大因素：冗余
冗余：组装出的Unigene的数量大大超过基因数；
冗余的来源：（1）可变剪切；（2）测序错误引入的“新” 转录本；
冗余的最大影响：导致多重比对，给定量带来困难
减少冗余策略
（1）去除低质量的reads;（2）是否有外援污染;（3）使用Normalization参数，降低高丰度基因的reads数据，同时提高组装效率；（4）后续序列聚类或过滤
去冗余方法

（1）筛选同一基因的最长转录本作为unigene

（2）软件：TGICL、CAP3通过聚类筛选unigene

# 1.提取最长转录本
/software/trinityrnaseq-2.2.0/util/misc/get_longest_isoform_seq_per_trinity_gene.pl \ 
./trinity_out_dir/Trinity.fasta > ./trinity_out_dir/unigene1.fasta

# 2.软件聚类去冗余
cd-hit-est -i ./trinity_out_dir/Trinity.fasta -o  output-cdhit -T 1 -M 1000

# unigene长度分布
perl /software/trinityrnaseq-2.2.0/util/misc/fasta_seq_length.pl ./trinity_out_dir/Trinity.fasta > ./trinity_out_dir/length.txt
#R画图
data <- read.table("length.txt",header=T)
data[,2][which(as.numeric(data[,2])>=2000)]<-2000

library(ggplot2)
pdf("length_distribution.pdf",height=7,width=10)
ggplot(as.data.frame(data), aes(x = as.numeric(data[,2])))+geom_histogram(binwidth =100)+
xlab("Transcripts Length Interval")+
ylab("Number ofTranscripts")+
labs(title="Transcripts Length Distribution")+
scale_x_continuous(breaks=seq(100,2000,by=100),
labels=c("100","200","300","400","500","600","700","800","900","1000","1100","1200","1
300","1400","1500","1600","1700","1800","1900",">=2000"))
dev.off()

Trinity组装后下游分析

STEP4: 定量：比对和丰度计算

利用RSEM进行丰度估计，

image

# 比对reads评估表达量（每个样品运行一次）
# Sp_log
/software/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts ./trinity_out_dir/unigene1.fasta --seqType fq \
--left RNASEQ_data/Sp_log.left.fq.gz --right RNASEQ_data/Sp_log.right.fq.gz \
--est_method RSEM --aln_method bowtie --trinity_mode --prep_reference \
--output_dir rsem_Sp_log_outdir
# Sp_ds
/software/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts ./trinity_out_dir/unigene1.fasta --seqType fq \
--left RNASEQ_data/Sp_ds.left.fq.gz --right RNASEQ_data/Sp_ds.right.fq.gz \
--est_method RSEM --aln_method bowtie --trinity_mode --prep_reference \
--output_dir rsem_Sp_ds_outdir

# Sp_hs
/software/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts ./trinity_out_dir/unigene1.fasta --seqType fq \
--left RNASEQ_data/Sp_hs.left.fq.gz --right RNASEQ_data/Sp_hs.right.fq.gz \
--est_method RSEM --aln_method bowtie --trinity_mode --prep_reference \
--output_dir rsem_Sp_hs_outdir

# Sp_plat
/software/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts ./trinity_out_dir/unigene1.fasta --seqType fq \
--left RNASEQ_data/Sp_plat.left.fq.gz --right RNASEQ_data/Sp_plat.right.fq.gz \
--est_method RSEM --aln_method bowtie --trinity_mode --prep_reference \
--output_dir rsem_Sp_plat_outdir

参数含义：
- --aln_method: 默认是bowtie
- --est_method:丰度评估方法，默认是RSEM，更准确，也可选择eXpress, 更快和所需RAM少
- --trinity_mode: 参考文件来自trinity，如果不是需要一个gene-isoform 比对文件 --gene_trans_map
输出结果
image

查看mapping结果

# ds
perl /software/trinityrnaseq-2.2.0/util/SAM_nameSorted_to_uniq_count_stats.pl rsem_Sp_ds_outdir/bowtie.bam > rsem_Sp_ds_outdir/mapping.out

STEP5: 差异分析

创建数量矩阵

/software/trinityrnaseq-2.2.0/util/abundance_estimates_to_matrix.pl \
--est_method RSEM \
--out_prefix Trinity_trans \
rsem_Sp_ds_outdir/ds_RSEM.isoforms.results \
rsem_Sp_hs_outdir/hs_RSEM.isoforms.results \
rsem_Sp_log_outdir/log_RSEM.isoforms.results \
rsem_Sp_plat_outdir/plat_RSEM.isoforms.results

无生物学重复进行差异表达分析

/software/trinityrnaseq-2.2.0/Analysis/DifferentialExpression/run_DE_analysis.pl \
--matrix trinity_trans.matrix/Trinity_trans.counts.matrix \
--dispersion 0.1 --method edgeR \
--output edgeR

查看差异数

sed '1,1d' edgeR/Trinity_trans.counts.matrix.ds_RSEM_vs_hs_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/ds_hs_DE_num

sed '1,1d' edgeR/Trinity_trans.counts.matrix.ds_RSEM_vs_log_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/ds_log_DE_num

sed '1,1d' edgeR/Trinity_trans.counts.matrix.ds_RSEM_vs_plat_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/ds_plat_DE_num

sed '1,1d' edgeR/Trinity_trans.counts.matrix.hs_RSEM_vs_log_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/hs_log_DE_num

sed '1,1d' edgeR/Trinity_trans.counts.matrix.hs_RSEM_vs_plat_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/hs_plat_DE_num

sed '1,1d' edgeR/Trinity_trans.counts.matrix.log_RSEM_vs_plat_RSEM.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' > edgeR/log_plat_DE_num

STEP6: 聚类/网络分析

提取差异表达的序列绘制热图

# 提取差异表达的序列绘制热图
cd edgeR/
/software/trinityrnaseq-2.2.0/Analysis/DifferentialExpression/analyze_diff_expr.pl \
--matrix ../trinity_trans.matrix/Trinity_trans.TMM.EXPR.matrix -P 1e-3 -C 2

# 根据聚类图提取子类
/software/trinityrnaseq-2.2.0/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl \
--Ptree 60 -R diffExpr.P1e-3_C2.matrix.RData

根据聚类图提取子类

# 提取差异表达的序列绘制热图
cd edgeR/
/software/trinityrnaseq-2.2.0/Analysis/DifferentialExpression/analyze_diff_expr.pl \
--matrix ../trinity_trans.matrix/Trinity_trans.TMM.EXPR.matrix -P 1e-3 -C 2

# 根据聚类图提取子类
/software/trinityrnaseq-2.2.0/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl \
--Ptree 60 -R diffExpr.P1e-3_C2.matrix.RData

STEP7: 功能注释和富集分析

使用TransDecoder-对trinity结果进行注释

1. 预测ORF (open-reading-frame) ( DNA -> protein sequence)

/software/TransDecoder/TransDecoder.LongOrfs -t trinity_out_dir/unigene1.fasta
# 
/software/TransDecoder/TransDecoder.Predict -t trinity_out_dir/unigene1.fasta

2. 预测基因功能

无参转录组序列组装 — Trinity

转录组分析策略

无参考转录组序列组装—Trinity概述

Trinity组装原理

Trinity 组装流程

STEP0: 准备工作

STEP1: 组装前质控

STEP2: 组装

STEP3: 组装后质量评估

Trinity组装后下游分析

STEP4: 定量：比对和丰度计算

STEP5: 差异分析

STEP6: 聚类/网络分析

根据聚类图提取子类

STEP7: 功能注释和富集分析

猜你喜欢

热点阅读