无参转录组序列组装 — Trinity

转录组分析策略

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

  • 有参

    • 比对—定量—差异分析—功能富集分析等
  • 无参

    • 组装—定量—差异分析—功能富集分析等
    image

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

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

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://upload-images.jianshu.io/upload_images/8242255-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. 预测基因功能

推荐阅读更多精彩内容