title: 转录组分析流程｜数据处理与De novo组装（一）
tags:
- 转录组组装
- 教程
- 软件
- Trinity
- Rcorrector
- Trimmomatic
categories: 转录组分析
转录组分析流程将分成三部分分别更新，更多内容关注我的个人博客

定义

转录组（transcriptome）广义上指某一生理条件下，细胞内所有转录产物的集合，包括信使RNA、核糖体RNA、转运RNA及非编码RNA；狭义上指所有mRNA的集合。
RNA-Seq (RNA-sequencing)也称为转录组测序，是最新发展起来的利用新一代测序技术进行转录组分析的技术，可以全面快速地获得特定细胞或组织在某一状态下几乎所有转录本的序列信息和表达信息，包括编码蛋白质的mRNA和各种非编码RNA，基因选择性剪接产生的不同转录本的表达丰度等。在分析转录本的结构和表达水平的同时，还发现未知转录本和稀有转录本，从而准确地分析基因表达差异、基因结构变异、筛选分子标记等生命科学的重要问题。

组装软件及用法

Reads前处理可以按照如下思路执行：

  1. 使用Rcorrector进行随机排序错误校正
  2. 删除无法纠正的reads
  3. 使用Trimmomatic去除测序接头和低质量序列
  4. 使用Bowtie2过滤细胞器reads（cpDNA、mtDNA 或两者）。将生成仅包含细胞器reads的文件，可用于组装例如带有Fast-Plast的质体
  5. 运行FastQC以检查reads质量并检测 over-represented reads
  6. 删除 over-represented reads

数据矫正

使用rcorrector软件对数据进行矫正,输入run_rcorrector.pl弹出使用说明

$ run_rcorrector.pl
Usage: perl ./run_rcorrector.pl [OPTIONS]
OPTIONS:
Required parameters:
	-s seq_files: comma separated files for single-end data sets
	-1 seq_files_left: comma separated files for the first mate in the paried-end data sets
	-2 seq_files_right: comma separated files for the second mate in the paired-end data sets
	-i seq_files_interleaved: comma sperated files for interleaved paired-end data sets
Other parameters:
	-k kmer_length (<=32, default: 23)
	-od output_file_directory (default: ./)
	-t number_of_threads (default: 1)
	-maxcorK INT: the maximum number of correction within k-bp window (default: 4)
	-wk FLOAT: the proportion of kmers that are used to estimate weak kmer count threshold, lower for more divergent genome (default: 0.95)
	-ek expected_number_of_kmers: does not affect the correctness of program but affect the memory usage (default: 100000000)
	-stdout: output the corrected reads to stdout (default: not used)
	-verbose: output some correction information to stdout (default: not used)
	-stage INT: start from which stage (default: 0)
		0-start from begining(storing kmers in bloom filter);
		1-start from count kmers showed up in bloom filter;
		2-start from dumping kmer counts into a jf_dump file;
		3-start from error correction.

上面说的很详细，下面在介绍几个常用到的命令

-s seq_files:单端测序的文件选择-s命令
-1 seq_files_left -2 seq_files_right:：双端测序的选择此命令输入左右两端测序的两个文件，目前大多为PE data
-k 设置kmer值，也可以直接不选择按照默认参数
-od 输出文件名称
-t 所用线程数根据自己情况进行选择

当然如果你只有一对数据的话可以单独进行操作，但是当你的转录组数据有很多时，这里推荐进行批量处理，具体的操作如下：

cat 文件名称 | while read f; do run_rcorrector.pl -t 12 -p ${f}"_1.fq.gz" ${f}"_2.fq.gz"; done

cat+自己创建的文件名，如果你的序列为

ABC_1.fq.gz,ABC_2.fq.gz;
BCD_1.fq.gz,BCD_2.fq.gz;
CDE_1.fq.gz,CDE_2.fq.gz;

其中文件的格式应该为

ABC
BCD
CDE

这里如果你的数据是_1.fq.gz以及_2.fq.gz则不需要修改，但是如果后缀和这里不一致，就修改成你自己数据的名称后缀
（注意，处理的所有数据的后缀名应保持一致）
跑完之后会在原有的*.fq.gz生成*.cor.fq.gz

去接头和低质量序列

去接头一般是需要去除掉测序时多余的接头，这一步的话如果测序数据下来之后已经做过了就可以选择可做可不做，说一下具体操作步骤
这里用到的软件为Trimmomatic
软件安装方法：
1、直接进入Trimmomatic官方下载
2、conda install -c bioconda trimmomatic

下载Illumina双端接头序列

curl -O -L http://dib-training.ucdavis.edu.s3.amazonaws.com/mRNAseq-semi-2015-03-04/TruSeq2-PE.fa

或者从官网下载的话自带有一般在adapter文件夹里 ~/Trimmomatic-0.39/adapters

$ trimmomatic
Usage: 
       PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...
   or: 
       SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...
   or: 
       -version

单端测序
trimmomatic SE -phred33 input.fq.gz output.fq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
双端测序
trimmomatic PE -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
或者将trimmomatic改为java -jar trimmomatic-0.35.jar运行

常用参数：
-threads 线程数
-trimlog 生成日志名，log文件较多建议不开启
-quiet 静默模式

修整步骤：
ILLUMINACLIP：从reads中剪切adapter和其他Illumina特定序列。
SLIDINGWINDOW：执行滑动窗口修剪，一旦窗口内的平均质量低于阈值，则切割。
LEADING：如果低于阈值质量，则在reads起始处剪切碱基
TRAILING：如果低于阈值质量，则在reads末尾处剪切碱基
CROP：将reads从末尾切割为指定长度
HEADCROP：从reads剪切后低于指定长度，则删除
MINLEN：如果reads低于指定长度，则删除
TOPHRED33：将质量得分转换为Phred-33
TOPHRED64：将质量得分转换为Phred-64

若数据有很多对转录组，建议批量操作

cat 文件名 | while read f; do time java -jar /home/你的路径/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 20 -phred33 ${f}"_1.cor.fq.gz" ${f}"_2.cor.fq.gz" ${f}"_1.cor.p.fq.gz" ${f}"_1.cor.u.fq.gz" ${f}"_2.cor.p.fq.gz" ${f}"_2.cor.u.fq.gz" ILLUMINACLIP:/home/lixingze/software/Trimmomatic-0.39/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:5 TRAILING:4 MINLEN:80; done

${f}"_1.cor.fq.gz" ${f}"_2.cor.fq.gz"为输入的文件
${f}"_1.cor.p.fq.gz" ${f}"_1.cor.u.fq.gz" ${f}"_2.cor.p.fq.gz" ${f}"_2.cor.u.fq.gz"为输出的文件
后缀一致，将名称放入文件中 cat 文件 即可,和上面形式一致。

转录组De novo组装

用Trinity进行转录组组装
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PYjf9TWX-1612082766472)(https://s3.ax1x.com/2021/01/24/sbRNvQ.png)]

Trinity包括三个模块，能处理大型的RNA数据
Inchworn: 将RNA-seq数据组装成单个转录本，通常是主要转录亚型的全长转录本
Chrysalis: 这一步将上一步得到contig进行聚类，对于每个聚类构建完整的德布罗意图(de Bruijin graph)每个转录本表示的是给定基因或者一组有着共同序列的基因的全部转录组成。 之后会根据图中不相交的点对全部短读数据进行拆分
Butterfly: 并行处理各个图(graph), 追踪每个图中的短读和配对短读的路径，最后报告可变剪切亚型的全长转录本，并且区分出旁系同源基因的转录本

安装方法

1、conda install -c bioconda trinity
2、trinity官网
因为trinity需要的依赖包很多，这里建议用conda安装,目前最新版本为Trinity-v2.11.0

操作步骤

安装成功后输入Trinity弹出如下

$ Trinity

###############################################################################
#

     ______  ____   ____  ____   ____  ______  __ __
    |      ||    \ |    ||    \ |    ||      ||  |  |
    |      ||  D  ) |  | |  _  | |  | |      ||  |  |
    |_|  |_||    /  |  | |  |  | |  | |_|  |_||  ~  |
      |  |  |    \  |  | |  |  | |  |   |  |  |___, |
      |  |  |  .  \ |  | |  |  | |  |   |  |  |     |
      |__|  |__|\_||____||__|__||____|  |__|  |____/

    Trinity-v2.11.0

#
#
# Required:
#
#  --seqType <string>      :type of reads: ('fa' or 'fq')
#
#  --max_memory <string>      :suggested max memory to use by Trinity where limiting can be enabled. (jellyfish, sorting, etc)
#                            provided in Gb of RAM, ie.  '--max_memory 10G'
#
#  If paired reads:
#      --left  <string>    :left reads, one or more file names (separated by commas, no spaces)
#      --right <string>    :right reads, one or more file names (separated by commas, no spaces)
#
#  Or, if unpaired reads:
#      --single <string>   :single reads, one or more file names, comma-delimited (note, if single file contains pairs, can use flag: --run_as_paired )
#
#  Or,
#      --samples_file <string>         tab-delimited text file indicating biological replicate relationships.
#                                   ex.
#                                        cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq
#                                        cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq
#                                        cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq
#                                        cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq
#
#                      # if single-end instead of paired-end, then leave the 4th column above empty.
#
####################################
##  Misc:  #########################
#
#  --include_supertranscripts      :yield supertranscripts fasta and gtf files as outputs.
#
#  --SS_lib_type <string>          :Strand-specific RNA-Seq read orientation.
#                                   if paired: RF or FR,
#                                   if single: F or R.   (dUTP method = RF)
#                                   See web documentation.
#
#  --CPU <int>                     :number of CPUs to use, default: 2
#  --min_contig_length <int>       :minimum assembled contig length to report
#                                   (def=200)
#
#  --long_reads <string>           :fasta file containing error-corrected or circular consensus (CCS) pac bio reads
#                                   (** note: experimental parameter **, this functionality continues to be under development)
#
#  --genome_guided_bam <string>    :genome guided mode, provide path to coordinate-sorted bam file.
#                                   (see genome-guided param section under --show_full_usage_info)
#
#  --jaccard_clip                  :option, set if you have paired reads and
#                                   you expect high gene density with UTR
#                                   overlap (use FASTQ input file format
#                                   for reads).
#                                   (note: jaccard_clip is an expensive
#                                   operation, so avoid using it unless
#                                   necessary due to finding excessive fusion
#                                   transcripts w/o it.)
#
#  --trimmomatic                   :run Trimmomatic to quality trim reads
#                                        see '--quality_trimming_params' under full usage info for tailored settings.
#
#  --output <string>               :name of directory for output (will be
#                                   created if it doesn't already exist)
#                                   default( your current working directory: "/home/lixingze/trinity_out_dir" 
#                                    note: must include 'trinity' in the name as a safety precaution! )
#  
#  --full_cleanup                  :only retain the Trinity fasta file, rename as ${output_dir}.Trinity.fasta
#
#  --cite                          :show the Trinity literature citation
#
#  --verbose                       :provide additional job status info during the run.
#
#  --version                       :reports Trinity version (Trinity-v2.11.0) and exits.
#
#  --show_full_usage_info          :show the many many more options available for running Trinity (expert usage).
#
#
###############################################################################
#
#  *Note, a typical Trinity command might be:
#
#        Trinity --seqType fq --max_memory 50G --left reads_1.fq  --right reads_2.fq --CPU 6
#
#            (if you have multiple samples, use --samples_file ... see above for details)
#
#    and for Genome-guided Trinity, provide a coordinate-sorted bam:
#
#        Trinity --genome_guided_bam rnaseq_alignments.csorted.bam --max_memory 50G
#                --genome_guided_max_intron 10000 --CPU 6
#
#     see: /home/lixingze/software/trinityrnaseq-v2.11.0/sample_data/test_Trinity_Assembly/
#          for sample data and 'runMe.sh' for example Trinity execution
#
#     For more details, visit: http://trinityrnaseq.github.io
#
###############################################################################

上面给的信息很详细运行只需要下面一条命令，接下来就是等待结果～

几个常用命令介绍
Trinity --seqType fq --samples_file 文件名称 --CPU 20 --max_memory 70G --output 输出文件夹

Trinity 
--seqType reads类型fq or fa
--samples_file 如果文件过多，可以将信息写入samples文件中
文件格式已给出：
cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq
cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq
cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq
cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq
--CPU 线程数
--max_memory 最大内存数根据自己服务器内存自己定
--output 输出文件名
--left : left reads，多个文件以逗号隔开
--right : right reads，多个文件以逗号隔开

查看trinity的完整的全部参数 Trinity --show_full_usage_info

数据解读

组装完成之后会在输出文件夹中输出一个Trinity.fasta这个就是最终组装结果！

>TRINITY_DN59_c0_g1_i10 len=4860 path=[0:0-614 2:615-757 4:758-815 5:816-882 7:883-919 9:920-1052 10:1053-1080 11:1081-1408 13:1409-1434 15:1435-1614 16:1615-1638 18:1639-1773 20:1774-1960 22:1961-2333 24:2334-2374 26:2375-2404 27:2405-2991 29:2992-4149 30:4150-4212 32:4213-4318 33:4319-4526 34:4527-4559 36:4560-4859]
CTTTGCCTCTCCTCAGAAGAGATGACAAGTGCATCAAGGAGTAGTTTTGCATCAACTACACATCATGCTTCGATGTTGCAGATAAGAGATTTCTCAGCATACCATATGGAAAATGTTTCCTACTATGCCTCCATGCATACTCAAAACGCCTCCTGAATTGGGGATAAATAGTTCTCGTATAGAATATAGGTTAATTATACAACAAATAAAGAACAAAGCTGTTAAACTGATTTTATGCATGAACACACTTTTTGA略

其中TRINITY_DN59_c0是Trinity聚类编号，g1是基因编号，i10是转录组亚型编号

取最长的片段

···
trinityrnaseqv2.11.0/util/misc/get_longest_isoform_seq_per_trinity_gene.pl 文件名
···

转录组后续分析

组装完成之后，就需要对其进行注释等后续操作，将会在下一篇中提到。
教程二

出自：lxz9.com，公众号：生信技术

# github # linux
Logo
华为开发者空间

华为开发者空间，是为全球开发者打造的专属开发空间，汇聚了华为优质开发资源及工具，致力于让每一位开发者拥有一台云主机，基于华为根生态开发、创新。

更多推荐

cover

一文解读GaussDB(DWS)监控运维诊断优化能力

avatar 华为开发者空间
cover

GaussDB技术解读——GaussDB架构介绍（四）

avatar 华为开发者空间
cover

华为云数据仓库专业级开发者认证重磅发布

avatar 华为开发者空间