RNA-seq预处理分析流程

本文描述了分子细胞科学卓越创新中心生物信息平台所开发的bulk RNA-seq分析流程，核心分析软件包括FastQC、STAR、Kallisto和RNASeQC等，用WDL流程语言编写，在HPC环境运行，经过简单修改也能在云平台运行。

分析流程示意图

针对Bulk RNA-seq，我们开发了快速的分析流程，如下图所示。其中测序质量控制（FastQC)、序列比对（STAR)、转录本定量（Kallisto）互相独立，均只依赖原始文件（FastQ)，因而会并行处理。STAR比对产生BAM文件会作为RNASeQC的输入。

分析模块

本流程主要包括Fastqc、STAR、 Kallisto和RNASeQC四个模块。

Fastqc

FastQC是Babraham Institute开发的二代测序质量控制软件。本流程采用FastQC软件来查看数据的测序质量，在开发者网站上还列出来了不同类型的数据质控的结果供大家参考。

  1. 好的Illumina数据
  2. 差的Illumina数据
  3. 接头污染的数据

对大多数数据包括RNA-seq而言，最重要的结果是单碱基测序质量（Per base sequence quality），原始数据的质量会对后续的分析结果产生影响，分析前会过滤质量差的碱基。

STAR

虽然目前RNA-seq比对软件众多，如比较流行的TopHat2及继任者hisat2，我们选择了STAR，它准确、快速，是众多分析流程的首选，请参考原始文献。STAR也是GATK RNA-seq所选择的工具，用于从RNA-seq数据中分析SNV/Indela.

Kallisto

Kallisto是非常优秀的RNA-seq定量软件，可以用于bulk和single-cell RNA-Seq 高通量数据的定量 。与STAR等依赖基因组序列的比对软件不同，Kallisto直接将测序片段直接比对到cDNA序列，因而也叫假性比对（pseudoalignment）。Kallisto能够在10分钟内完成一套含有30M个读段的双端RNA-seq数据的定量分析，是名副其实的桌面级分析软件。

RNASeQC

RNASeQC是由Broad研究所开发的RNA-seq质控软件，能快速输出比较全面的质控信息。

输入

为了运行这个流程，我们需要原始序列文件和基因组及注释文件。

序列文件

本流程每次处理一个样本，以gzip压缩的FastQ文件作为输入，双端测序为2个文件，单端测序为1个文件。 FASTQ格式文件中每个read由四行描述，如下：

@EAS139:136:FC706VJ:2:21041624:2182 1:Y:18:ATCACG   
CTCGTTCAAATGGACATAAAAGTCGATGCTCATCAGCTTTACAAGAATG
+
CADFFFHHHHGJJJJJJJJICHHIJJJJJIIJGHIGIJJJJJIGIJJJJ

其中第一行是以“@”开头的read测序标识符；第二行是碱基序列；第三行是以“+”开头的占位符；第四行是对应序列的测序质量。

基因组及注释文件

本流程中用到的人和小鼠的基因组信息及注释均来自GENECODE。比如人的基因组及注释信息包括如下三个文件。

  1. 参考序列：GRCh38.primary_assembly.genome.fa
  2. 基因注释：gencode.v34.primary_assembly.annotation.gtf
  3. cDNA序列：gencode.v34.transcripts.fa

JSON配置文件

输入的信息可以方便的通过JSON文件来配置，一个例子如下：

In [1]:
! cat /sibcb2/bioinformatics/WPS/Run_RNAseq/JSON/1302003-C.json

{
	"RNAseq.masterFolder":"/sibcb2/bioinformatics/WPS/Run_RNAseq/workflow",
	"RNAseq.output_basename":"1302003-C",
	"RNAseq.read1":"/sibcb1/bioinformatics/shijiantao/WPS/Raw_Platform/Fastq_RNA/1302003-C_R1.fq.gz",
	"RNAseq.read2":"/sibcb1/bioinformatics/shijiantao/WPS/Raw_Platform/Fastq_RNA/1302003-C_R2.fq.gz",
	"RNAseq.genomeIndexSTAR":"/sibcb2/bioinformatics/iGenome/STAR/GENCODE/human_hg38/STAR_Idx",
	"RNAseq.gtf_file":"/sibcb2/bioinformatics/iGenome/STAR/GENCODE/human_hg38/raw/gencode.v34.primary_assembly.annotation.gtf",
	"RNAseq.IndexKallisto":"/sibcb2/bioinformatics/iGenome/Kallisto/GENECODE_hg38/transcriptome.idx",
	"RNAseq.thread":"8",
	"RNAseq.collapsed_gtf":"/sibcb2/bioinformatics/iGenome/STAR/GENCODE/human_hg38/RNA-SeQC/gencode.v34_collapsed.gtf"
}

输出

本流程会在结果目录中生成四个结果文件夹分别为bamCoverage、kallisto、report、STAR。结构如下:

In [6]:
! tree /sibcb2/bioinformatics/WPS/Run_RNAseq/workflow/1302003-C

/sibcb2/bioinformatics/WPS/Run_RNAseq/workflow/1302003-C
|-- STAR
|   |-- 1302003-C.bam
|   |-- 1302003-C.bam.bai
|   |-- 1302003-C_STAR_Log.txt
|   `-- 1302003-C_STAR_QC.txt
|-- bamCoverage
|   `-- 1302003-C.bw
|-- kallisto
|   |-- 1302003-C.bam
|   |-- 1302003-C.bam.bai
|   |-- 1302003-C.h5
|   |-- 1302003-C.tsv
|   `-- run_info.json
`-- report
    |-- 1302003-C.bam.metrics.tsv
    |-- 1302003-C_R1_fastqc.html
    |-- 1302003-C_R1_fastqc.zip
    |-- 1302003-C_R2_fastqc.html
    `-- 1302003-C_R2_fastqc.zip

4 directories, 15 files

质量控制信息

STAR可以产生丰富的质量控制信息供我们参考：

In [3]:
! cat /sibcb2/bioinformatics/WPS/Run_RNAseq/workflow/1302003-C/STAR/1302003-C_STAR_QC.txt

                                 Started job on |	May 20 19:56:39
                             Started mapping on |	May 20 20:20:42
                                    Finished on |	May 20 20:40:35
       Mapping speed, Million of reads per hour |	86.88

                          Number of input reads |	28790824
                      Average input read length |	300
                                    UNIQUE READS:
                   Uniquely mapped reads number |	27404635
                        Uniquely mapped reads % |	95.19%
                          Average mapped length |	298.37
                       Number of splices: Total |	25657037
            Number of splices: Annotated (sjdb) |	25656468
                       Number of splices: GT/AG |	25373678
                       Number of splices: GC/AG |	222458
                       Number of splices: AT/AC |	26511
               Number of splices: Non-canonical |	34390
                      Mismatch rate per base, % |	0.33%
                         Deletion rate per base |	0.01%
                        Deletion average length |	1.97
                        Insertion rate per base |	0.01%
                       Insertion average length |	1.81
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |	1066478
             % of reads mapped to multiple loci |	3.70%
        Number of reads mapped to too many loci |	2391
             % of reads mapped to too many loci |	0.01%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |	0.00%
                 % of reads unmapped: too short |	0.94%
                     % of reads unmapped: other |	0.16%
                                  CHIMERIC READS:
                       Number of chimeric reads |	0
                            % of chimeric reads |	0.00%

一致性验证

为了验证kallisto统计的表达量是否正确，我们用featureCounts软件来统计基因的表达量（输入为STAR生成的bam文件），然后统计这两款软件的统计结果的相关性，经过验证两个结果相关性很高，结果见下图，横纵坐标的值为表达值的对数值（以10为底）：

Drawing

In [ ]: