1 生物信息学科研 博士ＧＷＡＳ

豆豆文献阅读第二天

2018-07-09  本文已影响195人  刘小泽

文章题目：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数据都在其中】

  1. 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运行结果

  1. 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以后会有介绍

上一篇 下一篇

猜你喜欢

热点阅读