使用Trinity拼接以及分析差异表达一个小例子

Trinity 将测序数据分为许多独立的de Brujin graph，理论上每一个图对应一个表达的基因。

整个流程分为三个步骤：Inchworm, Chrysalis, and Butterfly

Inchworm: 从reads中提取所有的重叠k-mers，根据丰度递减的顺序检查每个k-mers，然后将重叠的k-mers延长到不能再延长，称为一个contig

Chrysalis: 将上一部生成的contig聚类，对每个类构建de Brujin graph

Butterfly: 根据构建的de Brujin graph ，寻找具有可变剪接的全长转录本，同时将旁系基因的转录本分开

Trinity的硬件需求：

Inchworm 和 Chrysails 步骤对内存的需求很大，官方给出的说法是大致为每一百万对PE reads需要1g内存

使用的转录组数据为 Schizosaccharomyces pombe ，共4个样本（left right 表示双端测序数据的两端）

在拼接时，可以将每个样本都拼接成一个转录组，但是更合理的方法是将所有样本的reads合在一起再进行拼接，所以先将这四个样本的reads合在一起。

% cat *.left.fq > reads.ALL.left.fq
% cat *.right.fq > reads.ALL.right.fq
#添加环境变量
% export PATH=/usr/local/tools:$PATH
#一种典型的使用方法入下
#其中参数SS_lib_type RF 表示数据是双端（RF or FR) 单端（F or R）
% Trinity --seqType fq --max_memory 1G --left reads.ALL.left.fq --right reads.ALL.right.fq --SS_lib_type RF --CPU 2

完成后会在当前的工作目录生成一个 trinity_out_dir 的文件夹，Trinity.fasta为最终拼接结果。

Trinity自带了一个脚本可以显示一些结果的基本统计信息，N50表示的意思如下图。

使用GMAP将拼接结果比对到参考基因组（有参考基因组的情况下）

#首先准备GMAP需要的参考基因组，参考基因组文件为genome.fa
gmap_build -d genome -D ./
#algin 拼接结果，保存为一个sam文件
gmap -n 0 -D . -d genome ./trinity_out_dir/Trinity.fasta -f samse > trinity_gamp.sam

使用samtools转换为BAM文件（binary sam 优点是占用磁盘空间小，运算速度快，一些对数据的排序或者提取命令需要转换为BAM文件)

使用tophat 将RNA-seq reads map到参考基因组

#准备参考基因组
bowtie2-build GENOME_data/genome.fa genome
#run tophat 将所有的reads比对到参考基因组上
tophat2 -I 300 -i 20 genome
RNASEQ_data/Sp_log.left.fq.gz,RNASEQ_data/Sp_hs.left.fq.gz,RNASEQ_data/Sp_ds.left.fq.gz,RNASEQ_data/Sp_plat.left.fq.gz
RNASEQ_data/Sp_log.right.fq.gz,RNASEQ_data/Sp_hs.right.fq.gz,RNASEQ_data/Sp_ds.right.fq.gz,RNASEQ_data/Sp_plat.right.fq.gz
#下面的IGV基因组浏览器需要先建立索引
samtools index tophat_out/accepted_hits.bam

使用基因组浏览器IGV （有GUI）查看trinity的拼接结果

igv.sh -g `pwd`/GENOME_data/genome.fa `pwd`/GENOME_data/genes.bed,`pwd`/tophat_out/accepted_hits.bam,`pwd`/trinity_gmap.bam

使用RSEM定量

除了拼接以外，Trinity还准备了一些脚本进行后续的比如定量，差异表达等一些分析。

#使用Trinity准备好的脚本先用bowtie
#align到拼接好的转录组，然后使用RSEM定量
#运行这个脚本后会产生两个文件 'Sp_ds.isoforms.results' and 'Sp_ds.genes.results'
#包含了Trinity 拼接的转录本（isoform) 和基因的raw counts数和标准化后的数值
${Trinity_home}/util/align_and_estimate_abundance.pl --seqType fq
--left RNASEQ_data/Sp_plat.left.fq.gz --right RNASEQ_data/Sp_plat.right.fq.gz
--transcripts trinity_out_dir/Trinity.fasta
--output_prefix Sp_plat --est_method RSEM --aln_method bowtie
--trinity_mode --prep_reference --output_dir Abundance_quantify/Sp_plat.RSEM
#然后再对其他三个样本进行同样的操作
#一个样本间的比较矩阵 ,结果产生一个后缀为 .counts.matrix的文件
#显示了每个样本在每个转录本（isoform)上的map的数目（raw count)
${Trinity_home}/util/abundance_estimates_to_matrix.pl --est_method RSEM --out_prefix Trinity_trans
Abundance_quantify/Sp_ds.RSEM/Sp_ds.isoforms.results
Abundance_quantify/Sp_hs.RSEM/Sp_hs.isoforms.results
Abundance_quantify/Sp_log.RSEM/Sp_log.isoforms.results
Abundance_quantify/Sp_plat.RSEM/Sp_plat.isoforms.results
#另外 Trinity_trans.TMM.EXPR.matrix 是消除了测序深度，基因长度，然后通过TMM方法标准化后的数值(假定其他大多数基因没有差异表达)

使用 EdgeR 分析差异表达基因

还是通过Trinity安装包里自带的脚本，不加参数运行会有基本参数的介绍

使用刚才获得的 Trinity_trans.count.matrix 文件

运行结果 '*.DE_results' 输出了运行edgeR 分离出来的差异表达的基因

logFC = log fold change

logCPM = log counts per million

#提取FDR<=0.005）
sed '1,1d' edgeR/Trinity_trans.counts.matrix.Sp_log_vs_Sp_plat.edgeR.DE_results | awk '{ if ($5 <= 0.05) print;}' | wc -l
#画热图，需要进入刚才的/edgeR文件夹作为工作目录
$TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl
--matrix ../Trinity_trans.TMM.EXPR.matrix -P 1e-3 -C 2
#-P 为p的阈值，-C 为fold change = 2^2 =4 倍