文章题目：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown

利用Hisat、Stringtie、Ballgown进行RNA测序中转录本定量

来源：Nature Protocols DOI：10.1038/nprot.2016.095

第一部分：词句经典

1. High-throughput sequencing of mRNA (RNA-seq) has become the standard method for measuring and comparing the levels of gene expression in a wide variety of species and conditions.
   【一句很通用的话介绍了一下NGS对于RNA-seq的意义】
2. RNA-seq experiments capture the total mRNA from a collection of cells and then sequence that RNA in order to determine which genes were active, or expressed, in those cells.
   【说明了RNA-seq的作用 p.s. 本文做的都是mRNA的研究】
3. generate enormous numbers of raw sequencing reads, typically numbering in the tens of millions,
4. three is the minimum number of replicates for valid statistical results

第二部分：文章结构

摘要前言：
RNA测序能够在一次试验中捕获成千上万基因的表达量，能产生原始reads数十M，通过检测每个基因测出的reads数目可以估算基因丰度，当有两个或多个重复时能够判断差异表达基因。另外利用RNA-seq可以检测注释文件中没有的变异基因，也可以寻找每个基因不同的调控表达途径。

分析RNA-seq一般需要四步:

1. 将reads比对到参考基因组上（这篇文章是以有参转录组为例）；
2. 将比对上的alignments组装成全长转录本；
3. 基因与转录本的表达定量；
4. 计算不同实验条件下，所有基因的表达差异

之前有一套公认的好用又快的处理转录组的工具Tophat2和Cufflinks是由约翰霍普金斯大学团队开发，后来被被创作团队淘汰，当大家不理解的时候，他们发声：别用之前的啦，不好用！不够快！～我们有了新的三件套！于是新版三件套工具就这么出来了，还命名为“Tuxedo ”～～真是应了一句话无敌是多么寂寞

Hisat相对于Tophat2比对和发现可变剪切位点速度更快，消耗内存更少；【可以提供gff基因注释文件，当然如果没有程序也会自动检测剪切位点】

StringTie负责拼接转录本，构建isoforms用于估算基因表达量；【先接收Histat传来的alignments进行转录本拼接，每组数据之间是独立的，拼接的过程中估算每个gene和每个isoform的表达量；拼接完以后，所有的拼接好的序列进行merge，这一步是必须的，因为某些样本中的转录本只是被reads部分覆盖，导致的结果就是仅有那些被覆盖的区域得到了拼接，利用merge可以降低这一误差。】merge完成的转录本重新返回给StringTie，计算转录本丰度，它使用的算法和一开始拼接的一样，但如果在merge过程中转录本的丰度虽然不变，但结构发生了改变，reads也要因此进行调整。StringTie还提供了对每个转录本进行read-count，这个数据会传给BallGown

BallGown利用StringTie拼接的结果，计算得出差异表达转录本【接受StringTie传来的转录本和定量数据，根据设置的不同的实验条件进行分组】它包含了某些R包可以直接可视化

这三件套流程支持带有时间线的实验(比如研究某些病的发育历期)以及两个以上实验条件下的结果比较
流程很容易操作，但需要会使用命令行及R脚本运行

当然，这三个工具并不是缺一不可：如果你是要做无参转录组，可以用其他工具如trinity拼接好以后传给StringTie；也有一些项目只需要验证已知的某些基因在样本中的表达量；另外即使对于模式物种，基因注释文件也不完善，因此可以结合许多其他工具使用。虽然这三种工具能够检测差异基因表达，但对于外显子的差异不能探索，好在可以结合使用其他工具如：DEXseq 、rMATS 、MISO 。

这三件套不包括数据预处理（去接头、去污染、去低质量序列），但可以用第三方软件FASTX (http://hannonlab.cshl.edu/fastx_toolkit) 、FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) ；
适用于二代illumina测序，对于三代Pac Bio、Nanopore的比对环节需要用其他工具；
Ballgown目前支持三个组装软件的结果：StringTie、Cufflinks、RSEM

本文带了实例数据：
人类的RNA-seq数据，选取了其中基因比较丰富的X染色体数据151M，相当于全基因组的5%

技术路线如下：

技术路线

具体方法：

1. HISAT（Hierarchical Indexing for Spliced Alignment of Transcripts）比对
   比对的工具目前常用的有Bowtie2、BWA，他们运行时使用了BWT（Burrows-Wheeler transform）的数据格式，这是一种压缩格式，能将参考基因组高度压缩。然后再使用特殊的index机制——Ferragina–Manzini (FM) index，因此能快速在基因组上搜索匹配reads的区段，速度高达每小时比对几百万条reads。

   以上方法在DNA比对中最为常用，例如组装基因组需要比对拼接。但是对于RNA来讲，有一个DNA中不存在的问题需要注意：成熟的mRNA将introns切除，因此许多RNA测序的reads要跨越内含子进行比对，因此一条短序列可能会被拆开比对到相隔1万bp甚至更远的位置 【人类平均内含子平均长度大于6000bp，长的有超过1M的】人类RNA测序采用的得到的是100bp reads，会有超过35%的reads能跨越式比对到多个外显子。
   要将一条跨越多个外显子的reads比对到基因组上，其难度要大于 本就是在一个外显子中的reads。

   Hisat使用了叠加的比对算法：一是全局比对，整个基因组建立index；辅助成千上万个小的局部index；
   这两种算法也是基于BWT/FM体系建立起来的。因此这种特制的比对算法使用时比BWA、bowtie要快好几倍，但是内存用量仅是他们的两倍。相比较于他的前任Tophat2，速度提升了50倍。

   如果比对人类参考基因组（3G大小）的话，内存需要最少8G，比较人性化，一般电脑就能跑起来；
   8G内存，20个样本，每个样本100Mreads，一天能跑完比对环节。
   用户可以在比对的时候，可以添加基因注释--描述已知基因的位置（包括它们外显子、内含子边界），这样比对会增加结果准确性，可以更快找到新的剪切位点、转录起始和终止位点

2. StringTie拼接、融合
   首先将比对完的reads划分成不同的基因座组，然后把每个基因座拼接成isoform。运用了network flow algorithm ，从reads数最集中(也就是表达量最高)的转录本开始，然后移除已拼接完的reads，再次重复这个过程，最后拼接成众多的isoform，结束的标志是：全部的reads都被拼接完或者剩下的reads数目低于某个限制；
   Network flow algorithm这个算法能够同时计算转录本reads丰度和外显子-内含子结构，因此它重构各个转录本的速度比以前版本Cufflinks更快，消耗内存更少，但最少要求内存8G。

   如果用户提供了注释文件，stringtie会以它为指导进行拼接，另外还会将拼接的基因根据已注释好的直接命名。注释文件对低表达量的基因很有帮助，因为reads数量太少，不好拼接得到该基因。拼接完后，可以用gffcompare检测有多少拼接的转录本或全部或部分地map到了已知的基因，并且统计有多少是全新的基因。【此环节可以跳过】

   想要比较不同样本基因表达量差异，如果通过一个一个基因寻找进行比较是很难的，但是如果将每个样本中所有的基因融合到一起，再比较样本的差异就方便多。StringTie提供了merge 的功能，将各个相似的样本中的所有基因融合起来，再比较融合后的大序列不同。即使某一个样本中缺少了某一段外显子，通过和其他样本的比较参考，也能将这段缺少外显子序列找到和它相似的进行merge，**最根本的目的就是：相似的拼接alignments进行merge，再进行比较。

   merge

   这个图就解释了merge的用处：绿色的sample1-4是来自不同的四个样本的部分拼接片段；黑色的是基因组注释文件的部分；sample1和2都和参考基因组一致，所以把他们merge到一起形成了转录本A；sample3和4与参考基因组不一致，但彼此之间是一致的，所以把他们merge到一起形成转录本B

   merge完成以后，会重新估算一遍merged transcripts表达量。

   StringTie可以不依赖于注释文件进行新的基因、转录本预测，注释文件只是个佐证

3. Ballgown差异表达分析
   这算是上游Hisat、StringTie与下游R/bioconductor的中间桥梁，他将数据进行标准化
   基本上差异表达分析都会包括下面几个步骤：(i) data visualization and inspection, (ii) statistical tests for differential expression, (iii) multiple test correction and (iv) downstream inspection and summarization of results.

   • 数据可视化和检查
   • 差异表达的统计分析
   • 多重检验校正
   • 下游检查和数据summary

   在R里面使用Ballgown处理需要得到：

   • 表型数据：关于样本的信息
   • 表达数据：标准化和未标准化的外显子，转录本，基因的表达数量
   • 基因信息：有关外显子，转录本，基因的坐标以及注释信息

   使用Ballgown：

   • 读入数据：将上游Stringtie输出的转录本表达量数据与表型数据结合起来【注意保证二者的ID号一致】
   • 丰度预测：以FPKM（fragments per kilobase of transcript per million mapped reads）为单位的丰度预测根据library size进行标准化。
   • 差异表达分析：FPKM的数据解读需要取log转化一下，再建立线性模型
   • 可以在基因，转录本，外显子水平上进行差异分析
   • 结果以表格形式展现，样本量大的话还有p值和q值

第三部分：数据测试

软件准备：
HISAT2 software (http://ccb.jhu.edu/software/hisat2 or http://github.com/infphilo/hisat2, version 2.0.1 or later) 【hisat2支持--sra-acc从NCBI SRA数据库直接下载数据，然后自动转为fa/fq格式】
StringTie software (http://ccb.jhu.edu/software/stringtie or https://github.com/gpertea/stringtie , version 1.2.2 or later)
SAMtools (http://samtools.sourceforge.net, version 0.1.19 or later)
数据准备：
下载RNA-seq测试数据
其中包括人类基因组X染色体RNA-seq，基因注释文件 。【X染色体参考基因组、测试RNA-seq数据都在其中】

• 解压chrX_data.tar.gz会得到四个文件夹：samples、indexes、genome、genes
  samples中包含2个种群(GBR (British from England) and YRI (Yoruba from Ibadan, Nigeria) )各6个sample（男女各3个重复），一共12个sample。然后使用双端测序，结果共有24个测序fq文件
  

  测序文件

  indexes 中是hisat2构建好的chrX索引文件；
  genome 中是参考基因组chrX.fa (GRCh38 build 81 );
  genes 中是chrX.gtf , 是从RefSeq数据库中下载的GRCh38的基因组注释文件

• hisat2构建索引
  我们其实下载好了indexes文件，但是如果从头开始自己建立索引应该怎么做呢？
  新建一个文件夹index，用于存放一会生成的chrX_tran索引文件

  1. 第一步利用hisat2的两个python程序,将剪切位点、外显子区域找出来：
     extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
  2. 第二步使用hisat2-build
     hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
     【如果注释文件不存在，--ss、--exon可以省略】
     这一步仅构建chrX的索引就需要内存9G左右，如果针对全基因组构建需要160G内存以上。
     没有注释文件就不用这么多内存，另外对chrX构建索引，使用一个CPU核心差不多10min；
     构建全基因组的使用8个核心，大概需要2h
     我这里运行的时间是00:16:43，可能服务器有其他程序在跑

• 比对环节
  原来文章中是一个一个样本运行的，这样很麻烦，所以我写了一个小脚本，提高效率

  1. 将每个sample的reads比对到基因组上：
     基本的程序使用是这样的：

     hisat2 -p 24 --dta -x index/chrX_tran -1 chrX_data/samples/chrX_data/samples/ERR188273_chrX_1.fastq.gz -2
     chrX_data/samples/ERR188273_chrX_2.fastq.gz -S ERR188273_chrX.sam
     

     hisat2有几个选项设置：
     -p :线程数
     --dta：输出和拼接工具兼容的alignments(此处默认Stringtie)
     --dta-cufflinks:如果用cufflinks拼接，就要用这个选项
     --known-splicesite-infile ：如果之前没有建立剪切位点索引，直接用下载的gtf作为索引是不行的，可以用这个选项添加剪切位点文件【一般还是建议按上一步提前建立好索引文件】
     -x : 索引文件

     【这里就有一个问题：因为有12组24个文件，所以我不可能一个一个去运行，这里写一个循环脚本，让程序自动运行】
     【问题的关键是如何让输入的-1、-2双端测序文件匹配到12个sample数据上？】

     观察发现这些文件的命名都是有规律的，不同的只是前边数字部分，如果我们可以把1.fastq.gz之前的部分（ERRxxxxxx_chrX_）提取出来，那么运行程序的时候只需要引用名字，再加上对应的1.fastq.gz或者2.fastq.gz就可以啦！但是还要注意的是提取的时候只需要用1.fastq或者2.fastq其中一类就行，否则提取出来的会有重复

     文件命名

     【脚本如下】
     for i in *1.fastq.gz;do
     names=${i%_*} 
     #上面👆这一句的意思是：匹配提取出文件_前面的内容（或许你会问：为什么数字编号后面也有_，却不会匹配到前面的数字呢？这是因为存在一个贪婪匹配原则，匹配越多越好）
     hisat2 -p 24 --dta -x ../../index/chrX_tran -1 ${names}_1.fastq.gz -2 ${name    s}_2.fastq.gz -S ../../align/${names}.sam
     done
     

  2. 用samtools进行排序和转换：

     vi sam2bam.sh
     
     for i in *.sam;do
     names=${i%.*}
     samtools sort -@ 8 -o ${names}.bam ${names}.sam
     done
     
     time sh sam2bam.sh
     【运行时间: 4m33.061s】
     

  3. Stringtie拼接转录本：

     vi assemble.sh 
     
     for i in *.bam;do
     names=${i%.*}
     stringtie $i -p 8 -G ../chrX_data/genes/chrX.gtf -o ${names}.gtf 
     done
     
     【参数：-G：参考注释文件；-p：线程数；-o：输出的合并转录本GTF格式】
     【运行时间：1m53.666s】
     

  4. Stringtie转录本合并：
     标准格式：stringtie --merge [Options] { gtf_list | strg1.gtf ...}

     先构建merge_list: 其实就是把所有的输出gtf文件放到一个文本文档里，方便merge程序调取，
     当然如果merge的转录本量少的话，也可以一个一个输入

     vi mergelist.sh
     for i in *.gtf;do
     echo $i > mergelist.txt
     done
     sh mergelist.sh > mergelist.txt
     

     运行merge：
     stringtie --merge -p 8 -G ../chrX_data/genes/chrX.gtf -o stringtie_merged.gtf mergelist.txt
     【运行时间：0m2.597s 】
     【查看结果有31027行】
     

  5. 检查拼接的转录本有多少匹配到了参考基因组注释文件【可选】
     使用gffcompare
     安装：http://ccb.jhu.edu/software/stringtie/gff.shtml or http://github.com/gpertea/gffcompare

     标准使用：
     gffcompare –r chrX.gtf –G transcripts.gtf
     -r：参考注释文件
     -G：拼接的转录本，也是要被比较的对象
     -o：指定输出的结果前缀，否则默认gffcmp

     这里使用的程序是：
     gffcompare –r ../chrX_data/genes/chrX.gtf –G –o merged stringtie_merged.gtf
     匹配的结果存放在merged.annotated.gtf 中，其中在附加列中有一项内容"class code"
     可以帮助我们快速判断转录本与参考基因组的关系
     

     class code

     还有一个文件gffcmp.stats ,他计算了不同的gene features (e.g., exons, introns, transcripts and genes) 的Sensitivity和Precision分析数据
     【其中，Sensitivity是指：参考注释文件中的gene在拼接转录本中正确重构的比例；
     Precision：拼接转录本与参考注释重叠的比例】
     还计算了拼接转录本中新外显子、内含子、基因座的数量!

     其他信息

6. Stringtie估算转录本表达量并生成Ballgown能使用的统计表格

   首先新建一个文件夹ballgown，并在其子目录下新建sample名的文件夹，一会用来存放各个统计数据,像这样：

   Ballgown文件夹

【如果你要一个个mkdir敲进去肯定很慢，而且还容易错】
【批量新建文件夹简单的办法： 把文件夹名都放在一个文本文档中，这里我们先用mergelist.txt复制过来就行，但是mergelist.txt不能直接用，因为它里面只有前半部分是我们需要的】
【本着能用一行命令绝不写脚本的原则：
cat mergelist.txt | sed 's/_.*f//' |xargs -n1 mkdir】
直接按我们的需求新建好了文件夹

按需求新建文件夹

 > 接下来运行stringtie
 > vi estimate_abundance.sh
 > 
 > for i in *.bam;do
 > names=${i%_*}
 > stringtie -e -B -p 8 -G stringtie_merged.gtf -o ../ballgown/$names/${names}_chrX.gtf $i
 > done
 > 
 > -e:只估算给定的参考转录本（就是合并之后的）【与-G搭配】
 > -G：参考注释文件【这里用stringtie_merged.gtf】
 > -B：生成Ballgown格式的表格，和output文件放一起【与-o搭配】
 > 运行结果如下：

stringtie运行结果

7. Ballgown差异表达分析

   7.1首先下载并加载R包

   下载：“ballgown”
   source("https://bioconductor.org/biocLite.R") 
   biocLite("ballgown")
   
   “devtools”
   install.packages("devtools")
   
   "genefilter"
   source("https://bioconductor.org/biocLite.R") 
   biocLite("genefilter")
   
   "RSkittleBrewer"
   devtools::install_github('RSkittleBrewer', 'alyssafrazee')
   
   这里也能够看出来R语言安装包的三种途径：CRAN的install.packages、bioconductor的biocLite
   以及devtools::install_github
   

   library(devtools) #包安装器
   library(RSkittleBrewer) #设置颜色
   library(genefilter) #计算平均数和方差
   library(dplyr) #排序整理结果

   7.2 第二步加载表型数据

   pheno_data = read.csv("geuvadis_phenodata.csv")
   
   这里使用测试文件自带的csv文件。当然实际运行时要自己用excel创建，其中包含了关于测序样本的信息
   csv文件(comma-separated values)就是逗号隔开的；
   

   表型数据

   7.3 第三步加载表达量数据
   ballgown可以识别Stringtie、Cufflinks、RSEM的数据文件，读取时需要设置三个地方
   数据存放文件夹dataDir、样本识别号samplePAttern，表型数据pData
   bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
   将数据打包存放到bg_chrX这个对象中，它的核心就是一个矩阵，每一行是基因，每一列是样本

   7.4 第四步过滤低表达基因
   根据方差，过滤掉低于1的
   bg_chrX_filt = subset(bg_chrX, "rowVars(texpr(bg_chrX)) > 1", genomesubset=TRUE)

   texpr是提取bg_chrX的表达量(reads_count)，rowVars对行(基因)求方差，意思就是：过滤掉只在一个样本中表达，其他的样本中都不表达的基因

   7.5 第五步确定组间有差异的转录本/基因
   例如，表型数据有两个变量：性别与种群。要找不同性别(也就是这里的组)之间表达量差异的转录本，需要用种群的数据进行矫正。
   7.5.1 使用stattest 进行数据检验
   #官方解释：Test each transcript, gene, exon, or intron in a ballgown object for differential expression, using comparisons of linear models

   results_transcripts = stattest(bg_chrX_filt, feature="transcript",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")
   
   #featureL：在哪个水平进行计算，可选"gene", "transcript", "exon", or "intron"
   #covariate协变量：对因变量有影响的变量。它不是研究者研究的自变量,不为实验者所操纵,但仍影响实验结果。设置性别为协变量，因为本实验主要是分析不同国家人群之间的差异。
   #adjustvars：主变量，即不同的国家人群--> 字符串表示
   #getFC：得到不同性别组之间利用无关变量矫正后的差异倍数(FC=Exp male/female) --> 向量表示
   #meas：measurement采用哪种计算方式
   

   当然也可以确定组间有差异的基因
   results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")

   7.5.2 在差异转录本中添加基因名称与ID号：
   results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)

   7.5.3 将转录本/基因按P值从小到大排序：
   results_transcripts = arrange(results_transcripts, pval)
results_genes = arrange(results_genes, pval)
   转录本的结果如下：

   排序后的转录本

   7.5.4 将转录本/基因结果写入csv文件：
   write.csv(results_transcripts, "chrX_transcripts_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_genes_results.csv", row.names = FALSE)

   7.5.5 找q value小于0.05的基因和转录本：
   subset(results_transcripts, results_transcripts$qval < 0.05)
subset(results_genes, results_genes$qval < 0.05)

   差异表达转录本的结果如下：【匹配到了两个已知基因的isoform】

   表达差异转录本

   差异表达基因的结果如下：

   
   

   表达差异基因

8. 数据可视化

8.1 【可选】图层设置—使图片更美观

tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)

8.2 展示样本间基因丰度(FPKM值)的分布，根据颜色区分性别

fpkm = texpr(bg_chrX,meas="FPKM") #这里提取的是转录本的FPKM，当然也可以提取其中的exon、gene等
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
#取log的目的是为了让数据能存在正值、负值，否则数据都为正不容易看出差别。这样低表达小于0，高表达大于0，不表达等于0

8.3 展示单个转录本信息

比如要展示第12条：
ballgown::transcriptNames(bg_chrX)[12] --> 得到转录本名称
ballgown::geneNames(bg_chrX)[12] --> 得到包含此转录本的基因名称
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))

单个转录本信息

总结：最好的有参转录组分析软件组合应该是Hisat2+Stringtie+Deseq2
关于Deseq2以后会有介绍