转录组分析流程｜数据处理与De novo组装（一）
title: 转录组分析流程｜数据处理与De novo组装（一）
tags:
- 转录组组装
- 教程
- 软件
- Trinity
- Rcorrector
- Trimmomatic
categories: 转录组分析
转录组分析流程将分成三部分分别更新，更多内容关注我的个人博客
定义
转录组（transcriptome）广义上指某一生理条件下，细胞内所有转录产物的集合，包括信使RNA、核糖体RNA、转运RNA及非编码RNA；狭义上指所有mRNA的集合。
RNA-Seq (RNA-sequencing)也称为转录组测序，是最新发展起来的利用新一代测序技术进行转录组分析的技术，可以全面快速地获得特定细胞或组织在某一状态下几乎所有转录本的序列信息和表达信息，包括编码蛋白质的mRNA和各种非编码RNA，基因选择性剪接产生的不同转录本的表达丰度等。在分析转录本的结构和表达水平的同时，还发现未知转录本和稀有转录本，从而准确地分析基因表达差异、基因结构变异、筛选分子标记等生命科学的重要问题。
组装软件及用法
Reads前处理可以按照如下思路执行：
- 使用Rcorrector进行随机排序错误校正
- 删除无法纠正的reads
- 使用Trimmomatic去除测序接头和低质量序列
- 使用Bowtie2过滤细胞器reads（cpDNA、mtDNA 或两者）。将生成仅包含细胞器reads的文件，可用于组装例如带有Fast-Plast的质体
- 运行FastQC以检查reads质量并检测 over-represented reads
- 删除 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，公众号：生信技术
转录组分析流程｜数据处理与De novo组装（一）相关推荐
- Trinity简介（1）--用于无参考基因组的转录组de novo组装
一. Trinity简介 Trinity,是由 the Broad Institute 开发的转录组de novo组装软件,由三个独立的软件模块组成: Inchworm,Chrysalis和Butte ...
- 二代测序技术相匹配的de novo组装工具
.目前,Velvet,ABySS,SOAPdenovo,VCAKE,SPAdes等[4-6]多种与二代测序技术相匹配的de novo组装工具应运而生,而如何在众多组装工具中,根据序列属性和具体要求来选 ...
- 空间转录组分析流程（使用Seurat对空间数据集进行分析）
空间转录组分析流程(使用Seurat对空间数据集进行分析) 因为每次打开这个网页都非常慢,所以我讲这个网页进行一个翻译,方便学习. 使用Seurat对空间数据集进行分析,可视化和集成 1.介绍 本教程 ...
- 生物信息学之rnaseq转录组分析流程--转换文件中的ensemble id到gene名
生物信息学之rnaseq转录组分析--转换文件中的ensemble id到gene名 如何解决转录组分析中count之后遇到ensemble id的问题 一个将ensemble id转换成gene名的 ...
- 转录组分析流程｜TransDecoder预测转录本的开放阅读框（二）
使用TransDecoder预测CDS TransDecoder按照其官网的说明,主要用于识别转录本序列中的潜在的编码区域,也就是预测CDS.转录本可以由RNA-Seq数据通过Trinity组装来的, ...
- 转录组分析流程｜基于salmon转录组批量定量流程（三）
TransDecoder那一步最终得到了*.cds序列,之后就需要用到salmon进行下面操作 salmon介绍 Salmon是不基于比对计数而直接对基因进行定量的工具,适用于转录组.宏基因组等的分析 ...
- 转录组分析流程：表达差异分析之edgeR
edgeR edgeR是非常经典的转录组表达差异分析软件. 样本量:72个转录组样本 library(edgeR) library(HTSFilter)fc <- read.table('cou ...
- 转录组分析流程：比对（有参）及统计Counts矩阵
1. 质控 fastqc * multiqc * trimmomatic_run.sh #去掉前9个碱基 trimmomatic_run.sh #! /bin/bash #history: # Gos ...
- 一文搞定细菌基因组De Novo测序分析
本文转自基因的生物信息学分析,链接 https://mp.weixin.qq.com/s/xWOlv5WVJ7LwTuRQDXmGzg 以一个细菌的测序数据为例子,介绍细菌基因组测序分析流程.本次实验 ...
最新文章
- 常用数据结构讲解与案例分析
- 相关计算机专业的英语文献,英文文献及翻译计算机专业.doc
- C++中substr()函数用法详解
- 这 3 个字是未来发展关键，不重视的企业，正在被淘汰
- 力扣算法题—045跳跃游戏二
- 古剑奇谭二服务器维护,《古剑奇谭二》10月4日例行维护更新公告
- 区块链项目-Lisk
- 地...地...地震了
- QToolButton设置背景无效的思考
- 【JEECG技术文档】JEECG 组织机构导入V3.7
- mysql 脚本安装工具_mysql 非安装版的一个自动安装脚本及工具(更新版)
- python写一个数据库的界面_Python GUI教程（十四）：在PyQt5中使用数据库
- 使用计算机在什么上传输,MODEM的作用是使计算机数据能在什么上传输
- 物联网卡产品的应用和拓展
- 计算机主机异常经常蓝屏,电脑频繁蓝屏怎么办
- vue中播放flv流视频
- java 登录界面加验证码_java 做登陆窗口，带有用户名和密码输入框和验证码。求修改...
- 为什么要学习Linux？
- 项目进度经常延误，该怎么破？
- 学习笔记：WEB安全防护
热门文章
- scrcpy 使用教程：将安卓设备投屏到 PC 端
- JAVA进阶学习书籍
- mysql记录锁（record lock）,间隙锁（gap lock）,Next-key锁（Next-key lock）亲测
- Linux学习——yum学习和光盘yum源搭建
- “科学”流言榜：空气炸锅做菜会致癌？锅表示“不服”
- Python 使用win32相关的库实现简单自动操作电脑QQ--（2，打开好友对话框以及其他操作）
- 代码开源！！行人检测与行人重识别结合 person search
- 9、IAR中断向量表与中断服务函数的编写
- Ubuntu 16.04 用户登陆界面 循环输入密码 进不去图形界面
- Python爬虫入门【3】：美空网数据爬取