新司机带你学RNA-Seq数据分析
注：阅读原文就是作者的简书原文连接，阅读效果更加！
前言
这次给大家带来的是16年发表在NATURE PROTOCOLS上面的一篇处理RNA-seq数据的文章：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
和它的12年发表于同一个杂志的兄弟文章：Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks都是NATURE PROTOCOLS上阅读量最大的文章。
NATURE PROTOCOLS
当然还有很多其它的介绍不止生信的方法文章，大家有时间可以去探究。
同时我这里就不再赘述RNA-seq的具体原理，有需要了解的请移步：一个简略的RNA-seq演示
至于软件的安装到官网下载，解压后将bin/添加进路径即可，这里不再做讲解。
注：所有操作皆在LINUX&R上完成，默认基本处理软件已经安装
本体介绍
大佬的文章
事实上作者团队一直致力于开发出更好的解决数据处理的软件，就目前12年推出的Tohat2和Cufflinks已经不是那么的令人满意了，所以他们又开发了 HISAT, StringTie and Ballgown三件套完成大家对于一个RNA-seq所有的幻想。#但实际上目前还存在着其他的性能很优秀软件也可以满足需求，例如STAR等等，但作为菜鸟来说只有先一步一步的学习了。
好了，我们先对这篇文章给出的处理方法有个大概的了解。
pipeline：
alignment of the reads to the genome
assembly the alignments into full-length transcripts
quantification of the differencesin expression levels of each gene and transcript
calculation of the differences in expression for all genes among the different conditions
An overview of the 'new Tuxedo' protocol
这个protocol首先从原始RAN-seq数据入手，输出数据包括基因list，转录本，及每个样本的表达量，能够表现差异表达基因的表格并完成显著性的计算。
第一步是使用HISAT将读段匹配到参考基因组上,使用者可以提供注释文件，但HISAT依旧会检测注释文件没有列出来的剪切位点
第二步，比对上的reads将会被呈递给StringTie进行转录本组装，StringTie单独的对每个样本进行组装，在组装的过程中顺带估算每个基因及isoform的表达水平
第三步，所有的转录本都被呈递给StringTie的merge函数进行merge，这一步是必须的，因为有些样本的转录本可能仅仅被部分reads覆盖，无法被第二步的StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本，方便下一步的对比
merge的数据再一次被呈递给StringTie，StringTie可以利用merge的数据重新估算转录本的丰度，还能额外的提供转录本reads数量的数据给下一步的ballgown
最后一步：Ballgown从上一步获得所有转录本及其丰度，根据实验条件进行分类统计
注意事项Warning：
HISAT, StringTie, Ballgown并不是唯一处理RNA-seq的方法，也并不能处理所有的RNA-seq数据。
例如质量控制或者去除测序污染/去除接头/去除低质量数据（可以使用FastQC以及FASTXtoolkit做到）
这个pipeline无法处理第三代测序的数据
Ballgwon可以使用stringtie/cufflinks/RSEM产出的数据，但是可能无法接受其它程序产出的结果
默认参数适用于小至三个，大至小数百的样本处理，对于特殊需求还需要参考manual
数据处理设计
Read alignment with HISAT
总体上来说HISAT利用了BWA和Bowtie的算法，同时解决了mRNA中不存在内含子难以比对的问题，比对上代主流RNA-seq比对工具能快50倍，同时需求更少的内存<8G（这就意味着你可以在PC上跑数据），20个样本，每个样本一亿reads的估计，用一台电脑一天之内能够跑完。使用者可以提供精确的基因注释来提高在已知基因区域的准确性，但这是可选项。
事实上，针对RNA-seq数据的align目前有很多工具，16年12月12日出版的NATURE METHODS上的一篇文章，利用分别使用人和一种叫malaria的寄生虫的数据，模拟出三种复杂度的数据集（T1、T2和T3）。对于复杂度的衡量，定义了三个尺度：
突变率
indel比例
错误率
T1复杂度最低，T3最高。
随后使用了目前主流的比对工具进行了比对。
因为RNA-Seq测序的特性，天然的会有一部分数据延伸到内含子区，这部分跨越外显子和内含子的reads就称为『junction reads』，所以RNA-Seq比对软件需要针对此进行优化，而文章做benchmark也考虑到这点。分两个水平做比较：base and read level，这点和一般的DNA aligner一样
junction-level
度量软件的结果时，用了两个值：recall： 所有base中正确比对的比例。Tophat2在T1 base&read的表现，recall是85%-95%，T2降到80%，T3更是雪崩式降至20%以下。这部分表现好的是Novoalign和MapSplice2。
precision：所有比对上的base中，正确比对的比例。对于T1和T2，大部分软件这个值都较高。
结果上图：
接下来看看调参（自行设置参数，而不是使用软件默认参数）比对，直接说结论吧，不管是否调参，表现robust的是CLC，GSNAP、Novoalign、STAR以及MapSplice2。针对复杂度高的T3数据，Tophat2调参和默认参数得到的比对率相差67%还多。
另外还有运行时间的比较，这轮表现好的是HISAT/HISAT2，其实也能看出，数据复杂度越高，比对耗时越长。
Transcript assembly and quantification with StringTie
RNA-seq的分析依赖于精准的对于基因isoform的重建以及对于基因相对丰度的预测。继承于Cufflinks，StringTie相对于之前开发的工具更为准确，需求内存和耗时也更少。
使用者一样可以使用注释文件来帮助StringTie运行，对于低丰度的数据比较有帮助。此时StringTie依旧会对非注释区域进行转录本的组装，所以注释文件在这里是可选选项。
转录本组装完成后，使用者可以利用gffcompare(StringTie工具包含)工具来得知有多少转录本和注释文件相同，又有多少新的转录本:
#input: gff or gtf formattranscripts.gtf
command line example:
$ gffcompare -G -r chrX.gtf transcripts.gtf
#-r : reference
#-G :compare all transcripts
#-o :prefixoutput file
# 1: gffcmp.annotated.gtf
# 显示StringTie组装的转录本与注释文件内的转录本的差别（会给每个转录本加入一个class code，我理解为一个标识，释义如下图）
output dile
# 2: gffcmp.stats# 根据注释文件显示组装结果的准确性和阳性预测率class codes
由于在实验中，我们会处理多个RNA-seq数据，单个样本里面的基因和iosforms与其它样本的数据很少相同，但是为了进行比较就需要它们以一个相同的形式进行组装，作者通过StringTie的merge工具解决了这个问题，将所有样本中发现的基因进行merge。下图的例子可以帮助理解。
example
如图所示，四个样本中组装得到四个转录本，最后merge得到A/B两个转录本。1/2样本与参考基因组相同merge得到A转录本，3/4样本相互吻合但与参考基因组不一样，所以merge得到B转录本。
在merge步骤之后，需要StringTie再运行一遍去重新估算merge得到的转录本的丰度。Chevy理解：这就是相当于重新定义了一个annotation file进行二次重新组装，相对于第一次组织可以提高准确度和敏感性。比已有注释文件的优势在于可以发现更多的isoforms。
Differential expression analysis with Ballgown
组装和定量完转录本之后自然需要进行基因表达差异分析，统计建模和可视化。R和Bioconductor提供了一揽子的工具处理这些任务，包括raw data的plot作图，数据的标准化，下游的统计建。此处使用Ballgown包作为承上启下的桥梁。
在R里面使用Ballgown处理需要得到：
表型数据。关于样本的信息
表达数据。标准化和未标准化的关于外显子，junction，转录本，基因的表达数量
基因信息。有关外显子，junction，转录本，基因的坐标以及注释信息
而大多数的差异表达分析都会包括下面几个步骤：
数据可视化和检查
差异表达的统计分析
多重检验校正
下游检查和数据summary
Ballgown包可以完成以上的几个步骤并且可以联合R语言的其它操作进行分析
Ballgown使用：
数据的读入：需要将StringTie输出的数据结合表型数据，这里需要保证表型数据的identifier和StringTie输出数据的一致，不然会报错
预测丰度的检查：以FPKM（fragments per kilobase of transcript per million mapped reads）为单位的丰度预测将会根据library size进行标准化。#极端差异此处需要引起注意
使用线性模型进行差异表达分析，由于FPKM对于转录本解读过于曲解，所以这里需要使用log转化处理数据，随后再使用线性模型进行差异分析。
Ballgown可以对time-course和fixed-condition数据进行差异分析，但是无法避免批次效应带来的误差。#使用stattest功能可以指定任何已知的cofounder
Ballgown可以帮助你在基因，转录本，外显子，junction水平上进行差异分析
结果会以表格形式展现，如果样本够多还有p值和q值
结果数据可以用来进行下一步的gene set分析（例如GSEA）等等
实战示例
准备阶段
实际上本轮操作就是按照文章给的示例数据走了一遍流程：
文章中，为了减少计算时间，方便读者重复，作者截取了一批实验数据中能够匹配到chromosomeX上的数据作为示例数据，chromosomeX是一个基因相对较为丰富的染色体，可以占到基因组的5%左右。
首先是下载数据随后解压：
$nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
$tar zxvf chrX_data.tar.gz其中包括如下数据：
文件数据
sample文件夹下有12个fastq格式的paired-end RNA-seq文件，从两地人群中（英格兰岛住民和约鲁巴住民）各取六个样，六个样又分为男女性别各三个（最少的生物学重复）。
sample
随后介绍一个骚命令：
HISAT2可以直接从NCBI下载sra格式的文件并作为输入文件格式
下面以ERR188245测序数据为例$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran --sra-acc ERR188245 –S ERR188245_chrX.sam
可以说是相当的帮人偷懒了
index文件夹下包含着HISAT2对于染色体X的index文件
index
当然HISAT2已经为各位老爷准备好了常见基因组的index文件和genes文件。
点我
genome文件夹下包含这染色体X的序列文件（GRCh38 build 81）
genome
genes文件夹下则包含着针对基因组的注释文件（来自于RefSeq数据库）
annotation
geuvadis_phenodata.csv和mergelist.txt则是用户着要自己创建表明数据关系的文件，这里作者提供出来方便使用。
如果需要建立官网上没有基因组的index，就需要需要使用Hisat2包里面的python脚本:
extract_splice_sites.py和extract_exons.py，从注释文件里面抽取出剪切位点和外显子信息$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon先提取剪切位点和外显子数据 $ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran 随后建立HISAT2 index
正式开始
1.将reads比对上genome
# 样本比对操作
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ~/RNA-seq/chrx/chrX_data/samples/ERRERR188044_chrX_1.fastq.gz -2 ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam &
# 数据太多我就写了个小脚本处理,
in sample dir/
for i in *1.fastq.gz;
do
i=${i%1.fastq.gz*};
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran \
-1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S ~/RNA-seq/chrx/chrX_data/result/${i}align.sam \
2>~/RNA-seq/chrx/chrX_data/result/${i}align.log &
done运行结果如下：
比对结果
2.将sam文件转化为bam文件
# 转化操作
samtools sort @ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
# 批处理,
in result dir/
for i in *.sam;
do i=${i%_align.sam*};
nohup samtools sort -@ 8 -o ${i}.bam ${i}_align.sam &
done结果如下：
sam to bam
3.组装转录本并定量表达基因
# 单文件操作
stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf ERR188044_chrX.bam
# 批处理,
in result dir/
for i in *.bam;
do
i=${i%.bam*};
nohup stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ${i}.gtf ${i}.bam &
done结果如下:
strigtie
4.将所有转录本合并
warning: 此处的mergelist.txt是自己创建的，需要包含之前output.gtf文件的路径
cat ~/RNA-seq/chrx/chrX_data/mergelist.txt
ERR188044_chrX.gtfERR188104_chrX.gtfERR188234_chrX.gtfERR188245_chrX.gtfERR188257_chrX.gtfERR188273_chrX.gtfERR188337_chrX.gtfERR188383_chrX.gtfERR188401_chrX.gtfERR188428_chrX.gtfERR188454_chrX.gtfERR204916_chrX.gtf
#因为就在结果目录下面操作，所以直接显示文件名即可
stringtie --merge -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o stringtie_merged.gtf ~/RNA-seq/chrx/chrX_data/mergelist.txt
#output文件即为
stringtie_merged.gtf5.检测相对于注释基因组，转录本的组装情况
gffcompare -r ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -G -o merged stringtie_merged.gtf
#输出文件信息可见上面的Transcript assembly and quantification with StringTie介绍6.重新组装转录本并估算基因表达丰度，并为ballgown创建读入文件
# 单文件操作
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188044_chrX/ERR188044_chrX.gtf ERR188044_chrX.bam
# 批处理, in result dir/
for i in *_chrX.bam;
do
i=${i%_chrX.bam*};
nohup stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/${i}/${i}_chrX.gtf ${i}_chrX.bam &
done结果如下：
table
7.R包的安装与载入
R语言需要安装ballgown包和一些接下来处理会用到的包（RSkittleBrewer/genefilter/dplyr/devtools）
事实上我还没搞清楚LINUX里面的R包安装问题，老是报error，所以这里我就用Windows下的R进行操作。
# 首先是几个包的载入>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)8.数据的读入
# 随后是数据的读入，CSV文件，我把所有文件放在桌面，上一步得到的ballgown文件夹直接拷到桌面
> setwd("C:/Users/DELL/Desktop")
> read.csv("geuvadis_phenodata.csv")
ids sex population
1 ERR188044 male YRI
2 ERR188104 male YRI
3 ERR188234 female YRI
4 ERR188245 female GBR
5 ERR188257 male GBR
6 ERR188273 female YRI
7 ERR188337 female GBR
8 ERR188383 male GBR
9 ERR188401 male GBR
10 ERR188428 female GBR
11 ERR188454 male YRI
12 ERR204916 female YRI
> pheno_data<-read.csv("geuvadis_phenodata.csv")
> bg_chrX=ballgown(dataDir = "C:/Users/DELL/Desktop/ballgown",samplePattern = "ERR",pData = pheno_data)
# dataDir指的是数据储存的地方，
# samplePattern则依据样本的名字来，
# pheno_data则指明了样本数据的关系，
# 这个里面第一列样本名需要和ballgown下面的文件夹的样本名一样，不然会报错10.滤掉低丰度的基因
# 这里滤掉了样本间差异少于一个转录本的数据
> bg_chrX_filt=subset(bg_chrX,"rowVars(texpr(bg_chrX))>1",genomesubset=TRUE)11.确认组间有差异的转录本
> result_transcripts=stattest(bg_chrX_filt,feature = "transcript",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_transcripts
feature id fc pval qval
1 transcript
1 0.941367186 7.353184e-01 9.468599e-012 transcript
2 1.207949354 8.668638e-01 9.744055e-013 transcript
3 1.013100019 9.920824e-01 9.986659e-014 transcript
4 0.387671042 5.233364e-01 9.189906e-015 transcript
5 0.615797877 3.324164e-01 9.117015e-01......12.确认组间有差异的基因
> result_genes=stattest(bg_chrX_filt,feature = "gene",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_genes
feature id fc pval qval1 gene
MSTRG.1 1.114999374 7.305975e-01 9.218157e-012 gene MSTRG.10 1.749747142 2.767538e-01 7.922755e-013 gene
MSTRG.1000 1.399358968 6.039065e-01 8.945639e-014 gene MSTRG.1002 0.937668743 8.825216e-01 9.816878e-015 gene
MSTRG.1003 1.044840944 6.155345e-01 8.977043e-016 gene MSTRG.1004 1.220626116 4.160214e-01 8.388688e-01.....13.对结果 result_transcripts增加基因名
> result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneI Ds(bg_chrX_filt),result_transcripts)
> result_transcripts
geneNames geneIDs feature id fc pval qval
1 . MSTRG.4 transcript 1 0.941367186 7.353184e-01 9.468599e-01
2 PLCXD1 MSTRG.4 transcript 2 1.207949354 8.668638e-01 9.744055e-01
3 . MSTRG.4 transcript 3 1.013100019 9.920824e-01 9.986659e-01
4 . MSTRG.4 transcript 4 0.387671042 5.233364e-01 9.189906e-01
5 . MSTRG.5 transcript 5 0.615797877 3.324164e-01 9.117015e-01
6 PLCXD1 MSTRG.4 transcript 6 0.630018786 2.938396e-01 9.002815e-01......14.按照P值排序（从小到大）
> result_transcripts=arrange(result_transcripts,pval)
> result_genes=arrange(result_genes,pval)15.把结果写入csv文件
> write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
> write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)16.筛选出P值小于0.05的transcripts和genes
> subset(result_transcripts,result_transcripts$qval<0.05)
geneNames geneIDs feature id fc pval qval
1 . MSTRG.531 transcript 1657 0.030820591 1.427864e-10 2.082377e-07
2 XIST MSTRG.531 transcript 1656 0.003014576 1.860927e-10 2.082377e-073 . MSTRG.531 transcript 1655 0.016144997 3.762791e-10 2.807042e-074 . MSTRG.531 transcript 1658 0.028308029 6.599039e-08 3.692163e-055 TSIX MSTRG.530 transcript 1654 0.078461818 1.690716e-06 7.567643e-046 . MSTRG.613 transcript 1848 7.391342987 1.249275e-05 4.659796e-037 . MSTRG.141 transcript 421 3.204058932 2.729898e-05 8.727874e-038 . MSTRG.618 transcript 1852 9.136857610 4.804244e-05 1.343987e-029 . MSTRG.777 transcript 2338 0.127674288 1.000751e-04 2.488533e-0210 KDM6A MSTRG.258 transcript 736 0.054212867 1.173824e-04 2.627018e-0211 PNPLA4 MSTRG.62 transcript 186 0.592778584 2.083496e-04 4.238966e-0212 RPS4X MSTRG.519 transcript 1611 0.597532121 2.493976e-04 4.651265e-02> subset(result_genes,result_genes$qval<0.05)
feature id fc pval qval
gene MSTRG.531 0.002396124 2.469125e-11 2.523446e-082
gene MSTRG.141 3.165966412 1.666206e-06 8.514310e-043
gene MSTRG.530 0.082749314 7.018439e-06 2.390948e-034
gene MSTRG.613 7.314423295 1.128589e-05 2.883544e-035
gene MSTRG.618 9.050399157 5.022017e-05 1.026500e-026
gene MSTRG.519 0.637382385 6.953432e-05 1.184401e-027
gene MSTRG.376 0.620693674 1.134561e-04 1.656458e-028
gene MSTRG.62 0.600686117 1.638615e-04 2.093330e-029
gene MSTRG.777 0.159178288 2.177206e-04 2.472338e-0210
gene MSTRG.622 7.870447732 3.737270e-04 3.527295e-0211
gene MSTRG.778 0.127712244 3.796501e-04 3.527295e-0212
gene MSTRG.229 1.412379185 4.200674e-04 3.577574e-0213
gene MSTRG.48 4.158987103 5.279898e-04 4.150812e-0217.数据可视化之颜色设定
# 赋予调色板五个指定颜色
> tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow')
> palette(tropical)# 当然rainbow()函数也可以完成这个任务
> palette(rainbow(5))18.以FPKM为参考值作图，以性别作为区分条件
# 提取FPKM值
> fpkm = texpr(bg_chrX,meas="FPKM")
#方便作图将其log转换，+1是为了避免出现log2(0)的情况
> fpkm = log2(fpkm+1)
# 作图
> boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')结果如下：
FPKM, male in blue, females in orange
19.就单个转录本的查看在样品中的分布
> ballgown::transcriptNames(bg_chrX)[12]
12 "NR_027232"
> ballgown::geneNames(bg_chrX)[12]
12 "LINC00685"
# 绘制箱线图
> plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
+ main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
+ ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",
+ ylab='log2(FPKM+1)')
> points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
+ col=as.numeric(pheno_data$sex))结果如下：
transcript
20.查看某一基因位置上所有的转录本
# plotTranscripts函数可以根据指定基因的id画出在特定区段的转录本
# 可以通过sample函数指定看在某个样本中的表达情况，这里选用id=1750, sample=ERR188234
> plotTranscripts(ballgown::geneIDs(bg_chrX)[1750], bg_chrX, main=c('Gene MSTRG.575 in sample ERR188234'), sample=c('ERR188234'))结果如下：
21.以性别为区分查看基因表达情况
# 这里以id=575的基因为例（对应上一步作图）
> plotMeans('MSTRG.575', bg_chrX_filt,groupvar="sex",legend=FALSE)结果如下：
MSTRG.575
临终讨论:
关于时间的使用：针对chrX的数据量的比对和组装，在PC上可以40min内搞定，可以说是比较快的了。作者论述如果将12个样本的全数据下下来做分析的话，8核&8G内存的配置花费12.5h可以搞定比对，花费8h可以搞定组装和表达分析。
RNA_seq的比对情况一览
转录组组装情况一览
其实这篇文章说到底还是在讲方法论，它只负责上游的数据处理到可以分析的一步，下游的分析和可视化还需要依赖自身的实验设计和研究目的。
官方给出的方案是：stringtie处理得到的数据是直接呈递给ballgown进行差异分析和可视化作图，但是我还不是很清楚输出文件是否可以无差别的被其它软件使用。（我也是菜鸟）
RNA-seq的文章其实很多，也不是非要走这个pipeline。另外我觉得ballgown的作图其实也没有很elegant，org
参考文献
Pertea M, Kim D, Pertea G M, et al. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown[J]. Nature protocols, 2016, 11(9): 1650-1667.