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

时间：2021-12-09 15:13:18
使用Trinity拼接以及分析差异表达一个小例子
 2017-06-12 09:42:47     293     0     0

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 ，寻找具有可变剪接的全长转录本，同时将旁系基因的转录本分开

https://github.com/trinityrnaseq

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

Trinity的硬件需求：

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

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

  1. % wget \
  2. http://sourceforge.net/projects/trinityrnaseq/files/misc/Trinity
  3. NatureProtocolTutorial.tgz/download
  4. #解压后得到如下的文件
  5. tar –xvf TrinityNatureProtocolTutorial.tgz

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

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

  1. % cat *.left.fq > reads.ALL.left.fq
  2. % cat *.right.fq > reads.ALL.right.fq
  3. #添加环境变量
  4. % export PATH=/usr/local/tools:$PATH
  5. #一种典型的使用方法入下
  6. #其中参数SS_lib_type RF 表示数据是双端（RF or FR) 单端（F or R）
  7. % 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拼接以及分析差异表达一个小例子

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

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

  1. % $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta

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

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

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

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

  1. samtools view -Sb trinity_gmap.sam > trinity_gmap.bam
  2. #排序，方便后续使用
  3. samtools sort trinity_gmap.bam trinity_gmap
  4. #建立索引，需要先排序，否则报错，产生.bai文件
  5. samtools index trinity_gmap.bam

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

  1. #准备参考基因组
  2. bowtie2-build GENOME_data/genome.fa genome
  3. #run tophat 将所有的reads比对到参考基因组上
  4. tophat2 -I 300 -i 20 genome \
  5. 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 \
  6. 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
  7. #下面的IGV基因组浏览器需要先建立索引
  8. samtools index tophat_out/accepted_hits.bam

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

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

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

使用RSEM定量

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

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

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

使用 EdgeR 分析差异表达基因

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

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

  1. > ${Trinity_home}/Analysis/DifferentialExpression/run_DE_analysis.pl \
  2. >  --matrix Trinity_trans.counts.matrix \
  3. >  --method edgeR \
  4. >  --dispersion 0.1 \
  5. >  --output edgeR

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

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

logFC = log fold change

logCPM = log counts per million

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

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

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

相关文章