生信步骤|转录组mRNA数据的有参组装

转录组的组装能够提供丰富的组学信息，是生物信息学中重要的基础步骤。测序的下机数据经过过滤，比对，排序，组装，最终得到转录本的全部序列信息。下面我们以稻瘟菌测序数据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