转录组分析策略
转录组分析策略根据有无参考转录组可以分为两类：
有参
比对—定量—差异分析—功能富集分析等
无参
组装—定量—差异分析—功能富集分析等
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
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. 预测基因功能