RNA-seq数据分析
简介
基因表达是功能基因组学研究的一个重要领域。基因表达与基因信息从基因组DNA到功能蛋白产物的流动有关。RNA-seq已成为一种标准的基因表达检测方法，尤其用于检测转录本相对丰度和多样性。已有研究表明，其可靠性可以与其他成熟的方法如定量聚合酶链式反应(qPCR)相媲美。
声明： 文中如有不足请留言讨论！
1 生物基础
吐槽： 系统性总结真的好困难，刚开始阅读文献也很吃力。但是行动是知识特有的果实，慢慢积累吧。
1.1 中心法则
中心法则对于学习生物的人来说再熟悉不过了，基因信息的流动，蛋白的产生等一系列生物学事件都基于中心法则，是基因表达分析的一条主线。遗传信息从双链基因组DNA模板到翻译后修饰蛋白的流动是每个阶段至关重要的分子特征。RNA-seq通常针对成熟的mRNA分子。
图中对于中心法则做了系统的总结，感兴趣可以click here。
1.2 RNA-seq Protocol
典型的RNA-seq流程包括从感兴趣的样本中分离RNA、建库、高通量测序产生数以百万计的reads(长度一般为30-300bp)、比对到参考基因组或转录组，差异表达分析、发现转录本亚型和其它的下游分析。下图很好的概括了RNA-seq数据产生的过程。
RNA-seq的应用：
检测所有基因在特定条件下的表达水平(发育阶段、不同组织、正常与疾病、药物治疗)。 发现新基因，可变剪切，基因突变和基因融合等。
1.3 RNA-seq总的路线图
这张图引自A survey of best practices for RNA-seq data analysis。包括了三个板块，分别是预分析，核心分析和高级分析。很好的概括了实验设计和拿到测序reads后的数据分析，并介绍了不同的分析路线和每一步的数据分析方法，可以说是一个很好的大纲。
更多的细节可以通过阅读文章中的引用进一步了解。
2 数据分析
注意：
- 这里省略了进入相关文件目录的步骤，大家分析时要注意。
- 本次数据是小鼠早期胚胎测序得到的双末端数据，单末端相关参数可以使用–help查看帮助文档。但分析总体流程相同。
我使用的是HISAT2+featureCounts+StringTie流程。
2.1 前期准备
2.1.1 创建目录&安装conda
创建目录
# 首先使用cd命令需要进入自己的目录
cd ~/
# 创建文件夹，用于存储参考基因组
mkdir -p /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/
# 建立软件安装目录
mkdir biosoft
# 此外还需要建立项目分析的目录以及备份文件，方便查找
conda安装
Conda 是一种通用包管理系统，旨在构建和管理任何语言的任何类型的软件。通常与 Anaconda (集成了更多软件包，https://www.anaconda.com/download/#download) 和 Miniconda(只包含基本功能软件包, https://conda.io/miniconda.
html) 一起分发。
# 从官网下载最新版Miniconda3安装包，但速度较慢
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 安装Miniconda3：安装过程只需要输入 yes 或者按 Enter
bash Miniconda3-latest-Linux-x86_64.sh
# miniconda3安装成功，并成功配置好环境变量
source ~/.bashrc
# 镜像设置
# 下面这四行配置清华大学的bioconda的channel地址，国内用户推荐
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
# 配置镜像成功
conda config --set show_channel_urls yes
conda config --set channel_priority flexible
# 查看配置文件
cat ~/.condarc
conda常用命令
conda create -y -n trans python=3 # 创建小环境并成功安装python3
conda info --e # 查看当前conda环境
conda-env list # 列出所有小环境
conda activate trans # 激活小环境
conda deactivate # 退出小环境
2.1.2 常用文件格式简介
1. SRA：NCBI SRA数据库存放格式
SRA（Sequence Read Archive）：SRA是一个数据库，NCBI为了解决高通量数据庞大的存储压力而设计的一种数据压缩方案。
一般使用fastq-dump或fasterq-dump来将其转换为fastq格式的数据，才能做后续分析。
2. FASTQ：高通量数据存储格式
FASTQ格式将序列和Phred质量存储在一个文件中。序列和质量得分皆由单个ASCII字符编码。最早由Sanger Institute开发使用，目前已经是高通量测序结果的标准格式。
在FASTQ文件中，一个序列通常由四行组成：
第一行由@开头，之后为序列的标识符和描述信息；
第二行为序列信息；
第三行以+开头，之后可以再次加上序列的标识和描述信息；
第四行为质量的分信息，与第二行的序列相对应，其长度必须与第二行相同。
Phred质量值简介：
3. FASTA：记录信息的格式
对于每条序列，首行以“>”开头，其之后是注释；
在首行之后，是以单字母标准编码表达的实际序列数据
FASTA的序列表达：
1） 核算编码：A,C,G,T,N,U,R,Y,K,M,S,W,B,D,H,V
2） 氨基酸编码：[A-Z],*（代表终止密码子）
fasta序列格式如下：
4. SAM/BAM：高通量数据比对存放格式
SAM文件是由比对产生的以tab建分割的文件格式，BAM是SAM文件的二进制压缩版本。使用samtools view -S -b -o my.bam my.sam可以将SAM文件转换为BAM文件。
5. BED：基因组浏览器常用格式
BED 文件格式提供了一种灵活的方式来定义的数据行，以用来描述注释信息。BED行有3个必须的列和9个额外可选的列。每行的数据格式要求一致。
必须包含的3列：
(1). chrom - 染色体名字(e.g. chr3，chrY, chr2_random)或scafflold 的名字(e.g. scaffold0671 )。
(2). chromStart - 染色体或scaffold的起始位置，染色体第一个碱基的位置是0。
(3). chromEnd - 染色体或scaffold的结束位置，染色体的末端位置没有包含到显示信息里面。例如，首先得100个碱基的染色体定义为chromStart =0 .chromEnd=100, 碱基的数目是0-99。
9 个额外的可选列:
(4). name - 指定BED行的名字，这个名字标签会展示在基因组浏览器中的bed行的左侧。
(5). score - 0到1000的分值，如果在注释数据的设定中将原始基线设置为１，那么这个分值会决定现示灰度水平（数字越大，灰度越高），下面的这个表格显示GenomeBrowser
(6). strand - 定义链的方向，’’+” 或者”-”。
(7). thickStart - 起始位置（The starting position atwhich the feature is drawn thickly）(例如，基因起始编码位置）。
(8). thickEnd - 终止位置（The ending position at whichthe feature is drawn thickly）（例如：基因终止编码位置）。
(9). itemRGB - 是一个RGB值的形式, R, G, B (eg. 255, 0,0), 如果itemRgb设置为’On”, 这个RBG值将决定数据的显示的颜色。
(10). blockCount - BED行中的block数目，也就是外显子数目。
(11). blockSize - 用逗号分割的外显子的大小, 这个item的数目对应于BlockCount的数目。
(12). blockStarts - 用逗号分割的列表, 所有外显子的起始位置，数目也与blockCount数目对应。
6. GFF/GTF
GFF(General Feature Format)是基于Sanger GFF2的一种格式。GFF有9个必需字段，这些字段必须用制表符分隔。如果字段用空格而不是制表符分隔，则将不能正确显示。GTF (Gene Transfer Format, GTF2.2)是GFF的一种扩展格式。
2.2 软件安装
这里介绍三种软件安装方法。
2.2.1 conda安装软件
# 可以一次安装多个软件
conda install -y -c hcc aspera-cli
conda install -y sra-tools
conda install -y fastqc trim-galore hisat2 multiqc samtools featureCounts
# 运行以下命令，检查软件是否安装成功
hisat2 --help
which ascp # 查看软件安装路径
2.2.2 预编译版本软件安装
# 以aspera为例
wget -c https://d3gcli72yxqn2z.cloudfront.net/connect_3_9_1_171801_ga/bin/ibm-aspera-connect-3.9.1.171801-linux-g2.12-64.tar.gz
# 解压并运行脚本
tar -xzf ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh
# 添加环境变量
export PATH=~/.aspera/connect/bin:\$PATH >>~/.bashrc
# source使变量生效
source ~/.bashrc
2.2.3 源码安装
# 首先使用wget下载
# 安装分三步：
# 配置
./configure --prefix=想要指定的路径
# 编译
make
# 安装
make install
# 具体的安装在下载之后可以查看readme文件，一般有详细的介绍
2.3 数据下载
我这里使用aspera下载sra数据。
首先我们通过阅读文章获取数据的GEO accession，在GEO数据库中找到对应的BioProject，这里建议在ebi数据库下载对应的SRR号的相关信号，并获取aspera下载链接以及md5值等相关信息。
# ------------------------------------------------------
# GEO accession：GSE70605
# ------------------------------------------------------
# 下载得到tsv文件首先需要使用awk命令提取文件内aspera链接，保存为srr.aspra文件。之后运行下面脚本
创建一个脚本,命名为aspera.sh：
```bash
# 找到openssh文件位置，替换以下代码文件路径
find /home/hwb/ -name asperaweb_id_dsa.openssh
#!/bin/bash
cat srr.aspera | while read id
do
ascp -QT -l 300m -P33001 -i ~/miniconda3/envs/aspera/etc/asperaweb_id_dsa.openssh era-fasp@${id} ./
done
# 运行脚本：
./aspera.sh
下载完成之后一定要检查数据完整性，不然分析白做
使用awk命令提取md5值，生成md5文件，运行以下命令检查文件完整性
md5sum -c md5
md5文件如下图：
## 2.4 sra2fastq
将sra文件转换为fastq文件
```bash
创建脚本文件，个人觉得这样很方便。方便下次使用
#!bin/bash
cat srr.list | while read id
do
nohup fastq-dump --split-3 $id -O ./ &
done
2.5 质控（QC）
这一步生成的是HTML文件，可以用浏览器打开查看。其中，各种颜色是各项标准分析结果：绿色代表"PASS"；黄色代表"WARN"；红色代表"FAIL"。
ls ../raw/*.fastq|xargs fastqc -t 10 -o ./
multiqc ./# 整合质控结果
2.6 去接头
去接头完成之后，一定要再次质控，确保数据tidy
同样，写个简单的脚本
#!/bin/bish
# paired sequence
# trim adapter
paste <(ls *1.fastq) <(ls *2.fastq) > config
# mkdir ../clean
cat config | while read id
do
arr=$id
fq1=${arr[0]}
fq2=${arr[1]}
nohup trim_galore --gzip --paired $fq1 $fq2 -o ../clean > ../clean/trim.log &
done
2.7 hisat2比对
我这里使用的是HISAT2，index是在HISAT2官网下载的。所以省去了建index的步骤。hisat2的index有八个文件。这里建议新手将参考基因组的fa文件，注释gtf文件和index文件放在同一文件夹下，不然报错很难找。
这里也是一个小脚本
#!/bin/bash
# hisat2 of paired aligned
# dir = pwd
# 制作文件，第一列为输出文件名，第二列为双末端数据第一个文件，第三列为第二个文件
paste <(ls *fq.gz | cut -d"_" -f 1 |sort | uniq) <(ls *1_val_1.fq.gz) <(ls *2_val_2.fq.gz) > configure
# mkdir ../sam_out
cat configure|while read id;do
arr=(${id})
sample=${arr[0]}
fq1=${arr[1]}
fq2=${arr[2]}
echo $sample $fq1 $fq2
nohup hisat2 -p 4 -k 1 --dta -x /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/genome_tran -1 $fq1 -2 $fq2 -S ../sam
_out/$sample.sam &
done
到此生成了sam文件，比对的结果文件，这个文件一般较大，我们需要将其转换为bam文件，以便下一步分析。当然，可以在这个脚本中再添加命令，直接生成bam文件。
# sam2bam
cd ../sam_out/
ls *sam|while read id;do nohup samtools view -b -F 4 -@ 4 $id | samtools sort -@ 4 -O BAM - > ../bam_out/$(basename $id "sam")mapped.sort.bam &
done
2.8 定量和转录本组装
raw counts：
featuresCounts软件用于统计基因/转录本上mapping的reads数，也就是用于raw count定量。该软件不仅支持基因/转录本的定量，也支持exon, gene bodies, genomic bins, chromsomal locations等区间的定量。此外，还有htseq-count 等。sam和bam文件均可作为输入文件。
# 单样本定量
featureCounts
-T 5 \
-t exon \
-g gene_id -a /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o counts.txt mapping.sam
# 多样本
featureCounts -t exon -g gene_id -a /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o counts.txt ./*.bam
raw count矩阵文件简单处理，就可以用R包DESeq2筛选差异基因，可以参考我的另一篇博客，转录组差异表达分析和火山图可视化。
StrintTie进行转录本组装和定量：
gtf结果文件中有coverage，TPM和FPKM。此外RSEM也可计算FPKM值。
cd ../bam
ls *bam|while read id;do stringtie $id -p 4 -G /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o ../out_gtf/$(basename $id "bam")gtf;done
至此，RNA-seq上游数据分析基本已经完成。当然，有许多的分析流程，比如salmon流程，tophat+cufflinks等。使用之前尽量比对各个流程的优缺点。
数据分析过程坑很多，踩过了才能真正入门啊！（个人观点）
最后，非常感谢Jimmy老师的视频分享！ 由于篇幅太长，本文对于软件的使用没有详细介绍，大家可以去相应软件官网或者使用–help命令详细了解。
##如有侵权请联系作者删除！
欢迎关注我的微信公众号。
参考
[1] 生信技能树jimmy老师一系列课程
[2] A survey of best practices for RNA-seq data analysis
[3] Informatics for RNA Sequencing: A Web Resource for Analysis on the Cloud
168
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
