RNA-seq转录组数据分析丨一套完整的案例流程

今天分享的学习笔记是一套转录组分析简单流程，适用于初学者入门阅读，从原始测序数据开始，经过质控、序列比对、定量表达、差异表达、功能富集等一系列分析步骤，最终获得基因表达信息。本文所有数据都经过特殊修改，仅供学习参考使用。

转录组是在特定时空条件下细胞中基因转录表达产物，广义的转录组包括信使RNA，核糖体RNA，转运RNA及非编码RNA，狭义上是指所有mRNA的集合，转录组分析能够获得不同基因的表达情况。

1. 数据来源

假设有两个不同组织（PR和SR），每个组织各区三个样本，一共六个样本，利用illumina平台进行转录组测序，得到双端测序数据。数据原始格式为.fq,共有12条测序数据文件（每个样本产生两条）

PR1_2.fq  PR2_2.fq  PR3_2.fq  SR1_2.fq  SR2_2.fq  SR3_2.fq
PR2_1.fq  PR3_1.fq  SR1_1.fq  SR2_1.fq  SR3_1.fq

创建工作文件夹，并新建子步骤目录：00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis

mkdir RNA-Seq_Practice 
cd RNA-Seq_Practice
mkdir 00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis   
#批量创建工作文件夹

2. 测序数据质量评估

利用fastQC软件对获得的fastq序列文件进行质量分析，生成html格式的结果报告，其中含有各项指标，以下用PR1样本为例。

fastqc *.fq        #批量对fq后缀文件运行fastqc程序
输出结果：PR1_1_fastqc.html  
Filename	PR1_1.fq  
File type	Conventional base calls
Encoding	Illumina 1.5
Total Sequences	105300 		 #总序列数
Sequences flagged as poor quality	0
Sequence length	90   			#序列长度
%GC	52     			  #GC碱基含量

质量评估报告，使用浏览器打开：

3. 过滤低质量序列

利用Trimmomatic软件除去序列文件中的接头（adapter），并对碱基进行合适的修改，然后对碱基进行修剪，对低质量的序列进行过滤。

time java -jar trimmomatic-0.39.jar PE  #trimmomatic是依赖java环境运行
-threads 1 
-phred33 PR1_1.fq PR1_2.fq 
-summary ../01_trimmomaticFiltering/PR1.summary 
-baseout ../01_trimmomaticFiltering/PR1 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36

运行程序对序列文件中低质量的序列进行过滤，将输出结果存储到01_trimmomaticFiltering，每一个样本序列会生成如下输出文件，summary文件包含过滤的总结信息。

-rw-r--r-- 1  19M Nov 28 14:29 PR1_1P
-rw-r--r-- 1    0 Nov 28 14:29 PR1_1U
-rw-r--r-- 1  19M Nov 28 14:29 PR1_2P
-rw-r--r-- 1    0 Nov 28 14:29 PR1_2U
-rw-r--r-- 1  282 Nov 28 14:29 PR1.summary

打开其中一个PR1.summary文件进行查看，其中Input Read Pairs表示过滤之前的reads数，Both Surviving Reads表示过滤之后的reads数。

$ cat PR1.summary  #查看总结文件
Input Read Pairs: 102300
Both Surviving Reads: 102300
Both Surviving Read Percent: 100.00
Forward Only Surviving Reads: 0
Forward Only Surviving Read Percent: 0.00
Reverse Only Surviving Reads: 0
Reverse Only Surviving Read Percent: 0.00
Dropped Reads: 0
Dropped Read Percent: 0.00

4. 比对到参考基因组

利用hisat2软件，将fasta序列比对到参考基因组。首先需要构建索引文件，下载或者拷贝参考基因组序列文件Ref和基因组注释文件，通过hisat2-build命令构建索引文件。

cd xx/RNA-Seq_Practice
cp -r xx/Ref .
cd Ref
gunzip -c chr.fa.gz > chr.fa  #解压参考基因组
time hisat2-build -p 1 chr.fa Chr   #建立索引文件