无参转录组序列组装 — 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

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. 预测基因功能

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,688评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,559评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,749评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,581评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,741评论 3 271
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,684评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,122评论 2 292
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,847评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,441评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,939评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,333评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,783评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,275评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,830评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,444评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,553评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,618评论 2 249

推荐阅读更多精彩内容