转录组的组装能够提供丰富的组学信息,是生物信息学中重要的基础步骤。测序的下机数据经过过滤,比对,排序,组装,最终得到转录本的全部序列信息。下面我们以稻瘟菌测序数据SRR081556作为例子演示转录组组装的过程。
1.软件准备
转录组拼接所需软件及功能如下:
hisat2
建立参考基因组索引,reads比对。
samtools
进行sam文件向bam文件的格式转化,排序。
stringtie
组装转录本
gffread
转化gtf文件到mRNA fasta文件,至此完成转录本组装。
以上所有软件均可采用conda下载,在此不过多赘述。
2.基因组数据及转录组SRR数据准备
首先我们下载所需要的基因组数据和RNAseq数据:SRR081556
。基因组序列信息和注释信息可以通过NCBI搜索基因组70-15名称找到详情页,复制对应链接即可下载。RNAseq数据可通过prefetch命令搜索SRR号下载,得到的SRR081556文件夹通过fastq-dump命令进行格式转化,产生fastq文件。--split-3
会自动识别下载的SRR是单端测序还是双端测序并进行相应的格式转化,可以说非常人性化了!
$ prefetch SRR081556 #下载转录组数据
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/495/GCF_000002495.2_MG8/GCF_000002495.2_MG8_genomic.fna.gz #下载参考基因组序列,而后mv重命名为70-15.BAC.fa
$ wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/495/GCF_000002495.2_MG8/GCF_000002495.2_MG8_genomic.gff.gz #下载参考基因组注释信息
#转化SRR数据格式为fastq
$ fastq-dump SRR081556 --split-3 --gzip -O ./
由于我们下载数据的时候并不知道作者上传的是rawdata还是cleandata,因此在拿到单端或者双端测序数据之后,首先要做的就是数据过滤。这里采用fastp软件进行过滤。
$ fastp -i SRR081556_1.fastq.gz -o SRR081556_1.fastq.cleandata -I SRR081556_2.fastq.gz -O SRR081556_2.fastq.cleandata -w 8 #-w设置进程数8。默认2,max=16,或者是--thread 8
3.hisat2对测序数据进行mapping
在使用Hisat2软件进行比对之前,需要获得参考基因组的index文件。我们先使用hisat2-build
为参考基因组构建索引。
$ hisat2-build 70-15.BAC.fa 70-15.BAC.fa
构建基因组索引后可以正式进行双端测序的reads mapping:
$ hisat2 -p 8 --dta --no-mixed -x 70-15.BAC.fa -1 SRR081556_1.fastq.cleandata -2 SRR081556_2.fastq.cleandata --no-unal -S SRR081556.sam 2>SRR081556.summary.txt
命令中所用到的参数如下:
-x
参考基因组索引文件的前缀
-1
reads1,可以是压缩格式.gz .bz2
-2
reads2
-S
输出文件名(sam文件)
-P
线程数
-t
打印加载索引文件和对齐读取所需的时间
--dta-cufflinks
出来的结果更适合cufflinks处理 (主要用于基因表达量的计算和差异表达基因的寻找。本文不涉及转录本定量,故而不使用此参数)
--no-unal
不记录没比对上的reads
--un-conc
在输出文件写出与参考基因组不一致的reads对
--no-mixed
不做非配对的reads比对
4.samtools将比对结果二进制转化,排序
二进制转化将有效压缩比对文件的空间,加快后继程序的运行效率。参数-bS
输入待转化的sam文件。
$ samtools view -bS SRR081556.sam -o SRR081556.bam
$ samtools sort SRR081556.bam SRR081556.sorted.bam #注意此步骤产生的结果文件名称是SRR081556.sorted.bam.bam
5.stringtie生成转录本注释信息
Stringtie是一个基因和转录本组装定量的软件,stringtie的输入文件有两个:经过排序的bam文件和参考基因组的注释文件。而后stringtie生成转录本注释信息。此处的SRR081556.sorted.bam.bam
即为刚刚排序好的bam文件。-o
代表输出文件,-p
代表使用的核心数,-G
代表参考基因组的注释信息。
$ stringtie SRR081556.sorted.bam.bam -o SRR081556.gtf -p 20 -G GCF_000002495.2_MG8_genomic.gff
6.gffread生成转录本序列
我们将得到的gtf文件转化为fasta文件。gffread会读取刚生成的SRR081556.gtf
文件,-g
输入参考基因组序列,-w
选择输出文件,保存于SRR081556.fa
文件中。从而得到RNAseq的全部转录本。
$ gffread -w SRR081556.fa -g 70-15.BAC.fa SRR081556.gtf