仓库源文，站点原文

title: "专题：RNA-Seq 上游分析学习" tags:

转录组测序的研究对象为特定细胞在某一功能状态下所能转录出来的所有 RNA 的总和，包括 mRNA 和非编码 RNA，相对于传统的芯片杂交平台，转录组测序无需预先针对已知序列设计探针，即可对任意物种的整体转录活动进行检测，提供更准确的数字化信号，更高的检测通量以及更广泛的检测范围，是目前深入研究转录组复杂性的强大工具，基于高通量测序平台的转录组测序技术能够全面获得物种特定组织或器官的转录本信息，从而进行基因表达水平研究、新转录本发现研究、转录本结构变异研究等。

<!--more-->

RNA-seq 概述

RNA-seq 是研究转录组应用最广泛，也是最重要的技术之一，RNA-seq 分析内容包括序列比对、转录本拼接、表达定量、差异分析、融合基因检测、可变剪接、RNA 编辑和突变检测等，具体流程和常用工具如下图所示，通常的分析不一定需要走完全部流程，按需进行，某些步骤可以跳过、简化等。

RNA-seq 中最常用的分析方法就是找出差异表达基因 (Differential gene expression, DEG)，在实验室中，标准流程就分为三步：

版本信息

一些名词解释

大致流程

软件安装

只需要使用 conda 就可以安装所有需要的软件，主要使用的软件有以下一些

先创建 conda 虚拟环境，安装所需要的软件，可以自行手动安装，也可以直接导入我的 conda 环境：

git clone https://github.com/catuncle4/NGS-analysis.git
cd NGS-analysis/environment
cp .condarc ~/
conda env create --file rna-seq-env.yml   # RNA-seq

创建完成后激活环境就可以使用了：

conda activate rna-seq

如果想自己安装也可以，用下列命令就可以：

conda install bioconda::sra-tools
conda install bioconda::fastqc
conda install bioconda::trimmomatic
conda install bioconda::samtools
conda install bioconda::hisat2
conda install bioconda::subread
conda install bioconda::multiqc

数据获取与预处理

测序数据下载

先在 NCBI 的 SRA 数据库搜索感兴趣的物种

选择符合自己要求的文章，找到下面 Runs 这里，点击 SRR 开头的编号

查看数据是否符合要求

文件是_1.fq.gz、_2.fq.gz 这种是双端测序数据，我们需要这种双端测序的数据来进行 RNA-seq 分析

如果数据是双端测序的，那么就复制 SRR 编号，使用 sratools 下载，例如：

prefetch SRR8956151

批量下载需要先建立一个 txt 文件，将 SRR 编号写进去，例如：

SRR5830630
SRR5830631
SRR5830632
SRR5830633
SRR5830634

然后使用下面的命令下载

prefetch --option-file SRR_Acc_List.txt

由于数据比较大，可以使用 nohup 命令挂在后台下载

nohup prefetch -O . $(<SRR_Acc_List.txt) &

刚刚下载好的数据是 sra 格式的，使用 sratools 将其拆分

fastq-dump --gzip --split-3 SRR25909836.sra

如果数据比较多，就写一个 bash 脚本

#!/bin/bash
mkdir SRR
cat SRR_Acc_List.txt | while read id; do mv -f ${id}/${id}.sra ./SRR; done
cd SRR
for i in *sra
do
    fastq-dump --gzip --split-3 ${i} -f ../fastqgz
done

参考基因组及注释文件

植物的我一般在 Ensembl Plants 下载，用 wget 或 curl 都可以，内存不大

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.57.gff3.gz

然后解压

gzip -d Oryza_sativa.IRGSP-1.0.57.gff3.gz
gzip -d Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

为了方便，我把两个文件分别重命名为 oryza_sativa.fa 和 oryza_sativa.gff3

mv Oryza_sativa.IRGSP-1.0.dna.toplevel.fa oryza_sativa.fa
mv Oryza_sativa.IRGSP-1.0.57.gff3 oryza_sativa.gff3

FastQc 质控

常用参数

fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

进行质量检测

这里我们使用 fastqc [文件名] 即可

fastqc SRR25909836_1.fastq.gz

当然，数据比较多的时候还是挂在后台批处理然后等着就行

nohup fastqc SRR*.fastq.gz &
multiqc ./fastqc_report # 将多个质量检测报告合并

程序运行完成后会输出一堆 html 文件和 zip 压缩包，html 是网页版报告，zip 是本地宝报告，下载到本地用浏览器打开就可以看到质量检测报告了

左侧 Summary 部分就是整个报告的目录，整个报告分成若干个部分

我们一般比较关心的是下面几个部分

Basic Statistics

每个位置的碱基的测序质量

Per base sequence quality

一般要求此图中，所有位置的 10%分位数大于 20，否则切除 20 以下的碱基，从而保证后续分析的正确性

Per tile sequence quality

每个 tile 测序的测序质量

检查 reads 中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色表示低于平均偏离度，偏离度小，质量好；越红表示偏离平均质量越多，质量也越差，如果出现质量问题可能是短暂的，如有气泡产生，也可能是长期的，如在某一小孔中存在残骸，问题不大，可以看到我的这个数据几乎没有瑕疵

Per sequence quality scores

每条序列质量得分的分布情况

假如我测的 1 条序列长度为 101bp，那么这 101 个位置每个位置 Q 值的平均值就是这条 reads 的质量值，我这个数据的质量不错，reads 大都集中在高分上，纵轴数值越大，该序列测序错误的可能就越小

Per base sequence content

统计 reads 每个位置 ATCG 四种碱基的分布

Per sequence GC content

序列平均 GC 含量分布

蓝红色线应该比较接近才比较好，当红色的线出现双峰，很有可能是混入了其他物种的 DNA 序列，比如我这张图

Per base N content

每个位置无法检测的值的比例，当测序仪无法确定是何种碱基时，用 N 表示

正常情况下，N 的比例是很小的，所以图上常常看到一条直线，但放大 Y 轴之后会发现还是有 N 的存在，这不算问题，当 Y 轴在 0%-100%的范围内也能看到“鼓包”时，说明测序系统出了问题

Sequence Length Distribution

序列测序长度分布

每次测序仪测出来的长度在理论上应该是完全相等的，但是总会有一些偏差，当测序的长度不同时，如果很严重，则表明测序仪在此次测序过程中产生的数据不可信，比如我的这个图中，150pb 是最主要的，其他的几乎没有，所以数据的质量还是比较高的

Sequence Duplication Levels

统计 reads 重复水平

测序本身就会产生重复 reads, 测序深度越高，reads 重复数越大；如果重复出现峰值，就提示可能存在偏差（如建库过程中的 PCR duplication）

fastqc 抽取 reads 文件前 200,000 条 reads 统计其重复情况，重复数目大于等于 10 的 reads 被合并统计，这也是为什么我们看到上图的中间那里略有上扬，大于 75bp 的 reads 只取 50bp 进行比较，由于 reads 越长错误率越高，所以其重复程度仍有可能被低估

Overrepresented sequences

过度重复出现的序列的统计信息，上图中没有

Adapter Content

衡量的是序列中两端 adapter 的情况

如果在当时 fastqc 分析的时候-a 选项没有内容，则默认使用图例中的四种通用 adapter 序列进行统计

上图中 adapter 都已经去除，如果有 adapter 序列没有去除干净的情况，在后续分析的时候需要先使用 cutadapt 等软件进行去接头

Trimmomatic 过滤低质量序列

常用参数

过滤低质量序列

我使用的是下面的命令，需要根据自己的文件进行调整

trimmomatic PE -threads 1 -phred33 SRR25909836_1.fastq.gz SRR25909836_2.fastq.gz -summary SRR25909836.summary -baseout SRR25909836.fastq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36

解释一下这些参数

-phred33 之后的两个是正向和反向的测序文件

会看到五个输出文件 SRR25909836_1P.fastq.gz、SRR25909836_2P.fastq.gz、SRR25909836_1U.fastq.gz、SRR25909836_2U.fastq.gz、oryza_sativa.summary

hisat2 比对

hisat2 可以快速准确地将测序得到的 RNA 片段（reads）比对到参考基因组，从而确定这些 RNA 片段在基因组上的精确位置，进一步可以用于基因表达量定量，剪接位点的检测等多种 RNA-Seq 分析任务

建立索引

hisat2 需要一个 index 索引才能进行比对，hisat2 提供了一些 index，但很少，只有人类、小鼠等基因组的，可以在下面的 ftp 地址中进行下载

ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data

由于这里我做的是水稻的，所以就需要自己建立索引，使用的是水稻的参考基因组序列，前面已经下载好了，使用下列命令建立索引

hisat2-build -p 4 oryza_sativa.fa oryza_sativa

建立完成之后就可以看到系统中多了一个 oryza_sativa 文件夹，cd 进去就可以看到 8 个以 ht2 为拓展名的文件，要将参考基因组序列文件也放到这个文件夹

常用参数

进行比对

我使用的命令如下（使用经过过滤的测序数据）：

hisat2 -x oryza_sativa/oryza_sativa -p 5 -1 SRR25909836_1P.fastq.gz -2 RR25909836_2P.fastq.gz -S SRR25909836.sam

注意-x 后跟索引文件，不加拓展名，保证 ht2 文件和 fa 文件的文件名一致即可，这里由于前面过滤后的序列是没有拓展名的，所以会提示 Warning: Unsupported file format，不影响结果

运行完毕后便得到 sam 文件，还会输出一段信息

z -S oryza_sativa.sam
21675765 reads; of these:
  21675765 (100.00%) were paired; of these:
    2281601 (10.53%) aligned concordantly 0 times
    18927559 (87.32%) aligned concordantly exactly 1 time
    466605 (2.15%) aligned concordantly >1 times
    ----
    2281601 pairs aligned concordantly 0 times; of these:
      248174 (10.88%) aligned discordantly 1 time
    ----
    2033427 pairs aligned 0 times concordantly or discordantly; of these:
      4066854 mates make up the pairs; of these:
        2387493 (58.71%) aligned 0 times
        1619885 (39.83%) aligned exactly 1 time
        59476 (1.46%) aligned >1 times
94.49% overall alignment rate

这些输出记录包含以下信息：

在 RNA-Seq 分析中，比对成功率是一个重要的质量控制指标， 94.49%的比对成功率表明，绝大部分读取序列都能够成功地比对到基因组上，这表示 RNA-Seq 实验和测序质量都相对较好

samtools 排序压缩

常用命令

排序压缩

我使用的是以下命令

samtools sort -n -@ 5 SRR25909836.sam -o SRR25909836.bam

运行完成后会得到一个 bam 文件

featureCounts 生成基因计数表

常用参数

计数统计

我使用的是以下命令

featureCounts -T 5 -t exon -g Name -a oryza_sativa.gff3 -o counts -p SRR25909836.bam

oryza_sativa.gff3 就是最初下载的注释文件，如果要统计多个文件的话，在-p 后面跟上就可以，会生成 counts、counts.summary 两个文件，

counts.summary 文件是计数统计情况

counts 文件是基因的具体信息

我这里只有一组数据，所以数量统计也只有一列，通常做 RNA-Seq 时是需要多组数据进行分析的，将 bam 文件跟在-p 参数后面即可。

参考