查看: 5298|回复: 53

[转录组] 转录组组装新方法 - HISAT, StringTie and Ballgown

  [复制链接]

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

发表于 2016.9.27 11:20:24 | 显示全部楼层 |阅读模式
本帖最后由 johnlcd 于 2016.9.27 14:43 编辑

最近在做转录组数据分析,看到了nature protocol上8月11最新发的转录组分析新流程的文章,拿出来和大家分享,自己也准备好好研究研究,希望对自己以后的分析有所帮助。文章中主要用了HISAT/HISAT2, StringTie 和 Ballgown这几个工具和软件包,都是开源免费使用的,可方便下载安装。
首先了解一下几个工具作用;

HISAT/HISAT2:比对工具了,类似于tophat2;
Stringtie:组装与定量工具。
Ballgown:为差异表达计算工具



主要流程图





具体步骤


1、创建index
首先利用下面脚本提取剪接信息(有参考GFF前提下,没有忽略此步):
$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
然后构建HISAT2 index:
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
The --ss and --exon options(没有第一步可以不写)。indexing requires 9 GB of RAM for chromosome X and 160 GB for the whole human genome. The amount of memory is much smaller if one omits annotation information. Indexing chromosome X using one CPU core takes <10 min. It should take ~2 h to build an index for the whole human genome using eight CPU cores.

2、开始比对
各样本分别比对参考基因组
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188044_chrX_1.fastq.gz -2
chrX_data/samples/ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188104_chrX_1.fastq.gz -2
chrX_data/samples/ERR188104_chrX_2.fastq.gz -S ERR188104_chrX.sam
将SAM 转换为BAM:
$ samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
$ samtools sort -@ 8 -o ERR188104_chrX.bam ERR188104_chrX.sam

3、组装转录本
$ stringtie -p 8 -G chrX_data/genes/chrX.gtf -o
ERR188044_chrX.gtf –l ERR188044 ERR188044_chrX.bam
$ stringtie -p 8 -G chrX_data/genes/chrX.gtf -o
ERR188104_chrX.gtf –l ERR188104 ERR188104_chrX.bam

4、合并各个样本
$ stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o stringtie_merged.gtf
chrX_data/mergelist.txt
chrX_data/mergelist.txt:各个gtf路径放在里面。

5、估计表达丰度
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o
ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o
ballgown/ERR188104/ERR188104_chrX.gtf ERR188104_chrX.bam

6、加载 Ballgown  R包
$ R
R version 3.2.2 (2015-08-14) -- "Fire Safety"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)

7、 加载表型数据.
An example file called geuvadis_phenodata.csv is included with the data files for this protocol (ChrX_data). In general, you will have to create this file yourself. It contains information about your RNA-seq samples, formatted as illustrated in this csv (comma-separated values) file.
>pheno_data = read.csv("geuvadis_phenodata.csv")

8、加载表达丰度文件其来源于stingtie
To do this, we use the ballgown command with the following three parameters: the directory in which the data are stored (dataDir, which here is named simply ‘Ballgown’), a pattern that appears in
the sample names (samplePattern) and the phenotypic information that we loaded in the previous step (pData). Note that once a Ballgown object is created, any other Bioconductor32 package can be applied for data analysis
or data visualization.
>bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)

9、过滤低表达基因
>bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)

10、鉴定差异转录本
>results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")

11、鉴定差异基因
>results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")

12、添加基因名字和geneID。
>results_transcripts =
data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)

13、按照P值从小到大排序。
>results_transcripts = arrange(results_transcripts,pval)
>results_genes = arrange(results_genes,pval)

14、保存到文件。:
>write.csv(results_transcripts, "chrX_transcript_results.csv",
row.names=FALSE)
>write.csv(results_genes, "chrX_gene_results.csv",
row.names=FALSE)

15鉴定 q value <0.05的转录本:
>subset(results_transcripts,results_transcripts$qval<0.05)
>subset(results_genes,results_genes$qval<0.05)

16 作图。
>tropical= c('darkorange', 'dodgerblue',
'hotpink', 'limegreen', 'yellow')
>palette(tropical)

17、对于基因按照FPKM 值作图
>fpkm = texpr(bg_chrX,meas="FPKM")
>fpkm = log2(fpkm+1)
>boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')



18、对单个基因在不同样本中表达情况作图。. For example, here we show how to create a plot for the 12th transcript in the data set . The first two commands below show the name of the transcript (NM_012227)
and the name of the gene that contains it (GTP binding protein 6, GTPBP6):
>ballgown::transcriptNames(bg_chrX)[12]
## 12
## "NM_012227"
>ballgown::geneNames(bg_chrX)[12]
## 12
## "GTPBP6"
>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))



19、输出一个样本中一个基因座位的所有转录本的基因结构与表达丰度图
>plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))



20、我们也可以使用plotMeans 属于一个基因的所有转录本的平均表达值
>plotMeans('MSTRG.56', bg_chrX_filt,groupvar="sex",legend=FALSE)

就这些了,总之:This protocol does not require programming expertise, but it does assume familiarity with the Unix command-line interface and the ability to run basic R scripts. Users should be comfortable running programs from the command line and editing text files in the Unix environment.




软件链接
HISAT2:http://ccb.jhu.edu/software/hisat2/index.shtml
StringTie:http://ccb.jhu.edu/software/stringtie/index.shtml
Ballgown:
1. 运行R;
2. source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x

评分

参与人数 2奥币 +40 贡献 +11 收起 理由
zhouyulu + 20 + 5 非常详细的,给力
基迪奥-周煌凯 + 20 + 6 这个分享很给力!

查看全部评分

回复

使用道具 举报

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

 楼主| 发表于 2016.9.27 11:28:03 | 显示全部楼层
附上软件链接:
HISAT2:http://ccb.jhu.edu/software/hisat2/index.shtml
StringTie:http://ccb.jhu.edu/software/stringtie/index.shtml
Ballgown:
1. 运行R;
2. source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")
回复 支持 1 反对 0

使用道具 举报

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

 楼主| 发表于 2016.9.27 12:41:07 | 显示全部楼层
本帖最后由 johnlcd 于 2016.9.27 13:09 编辑

看到了一篇帖子,对三个软件的算法功能进行了总结,也和其他转录组分析软件进行了比较,比如tophat,cufflinks,cuffdiff等,可供参考理解。
PS:重点在文末


##############################################
HISAT:
1:index算法基于BWT和以BWT为基础的FM index: FM 算法是以BWT为基础,在计算的过程中加了两个参数,一个是OCC,
Occ[c,r]表示在BWT(T)中第r行之前出现字符c的个数,因为如果把整个基因组存进去,每次都要重头数一遍,消耗量太大,所以是以几百行为一组为一个check point,这样内存就小了,也方便检索;另外一个是SA, 记录第r行在参考基因组中是什么位置,这个过程是在BWT操作中实现的
2:建立index的方式:全基因组FM index和局部index(特色,新的建索引库策略)核心算法会用到tophat3上
HISAT在基因组范围内有48000个局部FM index
每个长度64k,通过测试数据和真实数据的比较,发现这是一个最快的软件
HISAT的基本设计原理:
HISAT利用bowtie2来建立很多低水平结果的FM index,这些index包含两种类型:1:全基因组index;2:大量的小FMindex,每个代表64k,也因为建了很多小的index,所以总的内存使用情况也较低
3 :100bp的reads至少跨越两个exon的占34.5%,分为三类:1:每个外显子至少16bp(25.1);每个外显子8-15bp(5.1),每个外显子1-7bp(4.2);
另外跨越三个的3.1%。针对每种reads都有相应的比对策略
基本的比对策略的思想是先在众多index中找,如果比对上多个地方,就往两边延伸;
4:基本的性能对比
速度:110,193条reads/s, tophat 1,954
利用模拟数据正确比对率(99.2%),tophat2(97.4)
在跨越2个exon,只有1-7bp的reads来说,
  uniq比对率HISAT(94.4)tophat2(77.8)
另外敏感度和准确性
HISAT(97.3,94.8) Tophat2(90.6,82.6)
100bp,20M的reads
用时26.7min,tophat2 1170分钟

#############################
1:StringTie和Cufflinks算法对比
cufflinks parsimony算法  (简约算法):生成最少的亚型
说这种算法没有考虑转录丰度,在isoforms方面算的不准。其在算表达量的时候,按照图上的说法是用了最大似然冗余算法。
stringTie先将reads分为不同的类,然后再针对每个类的reads生成一个拼接图来确定转录本,之后每个转录本产生一个流神经网络的最大流算法来评估表达水平
http://www.cnblogs.com/ZJUT-jiangnan/p/3632525.html
这个算法的意思对应过来就是在一个基因处的若干个转录本,如何分配reads的数目才能让每个转录本的数目都处在最多的状态,发的链接上是用了一个工厂的例子。
这个算法是求解最优化的,在这篇文章上也是第一次用。
2:在RNA组装方面优劣势对比
在组装方面StringTie具有一些优势,在低表达的部分,阈值过滤5%的StringTie比阈值过滤10%的准确度和敏感度还要高(这里的准确度和敏感度是把原始数据随机抽取出来一部分数据,看看这两组随机抽出来的数据的重合度如何和ROC是一个事情)
关于组装效果,StringTie要好于cufflinks,同时他们又远远好于其他软件,就组装错误而言,文章上列举了一个例子,即附表四,cufflinks预测出来的FPKM值可以到293509.68,而StringTie最大是935,这个问题,我们在项目处理中也遇到过
3:性能对比
时间上来说:StringTie 30min      cufflink 81min 其他更多
真实数据测试35-76分钟,比其他快3倍
内存比cufflink少一半
找出来的基因中,cufflink找出来的70%在StringTie中有重合,相比于cufflink,StringTie在基因重构方面对三种类型的基因更有效,分别是:低冗余,高exon数目,和多重转录本。
按照作者的说法,StringTie之所以能取得好的效果,是因为模型中有个求最有的过程;其可以通过
转录组的每个组装平衡覆盖度,算法自己可以通过合并覆盖到的深度来限制组装。
作者倒是希望通过StringTie来取代cufflink

########################################

关于tablemaker和ballgown
tablemaker的设计是为了处理cufflinks和cuffmerge的结果,一遍后面用ballgown处理,不过按照文章上说的ballgown也可以直接输入,可以是针对StringTie的结果
ballgown是一个R脚本,用来分析差异结果的,差异分布的结果是基于F-test来检测的,和我们以前用的有些不同;
就测试内存而言,cuffdiff用148G6 9个小时,而ballgown用18秒8G,

就差异结果而言,cuffdiff更保守一些,文章中用了一个例子,就是男性Y染色体的表达量数据和女性落到Y染色体的表达量数据对比,ballgown可以起到很好的效果,而cuffdiff找到的差异基因少。

所有这些软件的开发目的就是为了能在笔记本上跑开生物程序



转自:http://blog.csdn.net/hill_night/article/details/44829965

评分

参与人数 1奥币 +20 贡献 +10 收起 理由
zhouqian2617 + 20 + 10 Good job!

查看全部评分

回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
52
奥币
5549
积分
1231
注册时间
2016.1.8
在线时间
264 小时

突出贡献优秀版主论坛元老


发表于 2016.9.27 14:45:54 | 显示全部楼层
给力
回复

使用道具 举报

帝王蝶

Rank: 4

主题
7
奥币
758
积分
243
注册时间
2016.9.5
在线时间
64 小时

发表于 2016.9.27 17:39:00 | 显示全部楼层
好巧,昨天刚开始看这篇。学习了。
回复 支持 反对

使用道具 举报

版主

Rank: 10Rank: 10Rank: 10

主题
19
奥币
4029
积分
1459
注册时间
2015.12.29
在线时间
356 小时

突出贡献优秀版主热心会员


发表于 2016.9.28 10:58:14 | 显示全部楼层
楼主能不能把geuvadis_phenodata.csv下载下来,放到帖子里。
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

 楼主| 发表于 2016.9.28 11:53:37 | 显示全部楼层
“An example file called geuvadis_phenodata.csv is included with the data files for this protocol (ChrX_data). In general, you will have to create this file yourself. It contains information about your RNA-seq samples, formatted as illustrated in this csv (comma-separated values) file. ”
那不是个示例文件吗,自己的数据自己建表型文件啊,有的数据就包括有phenotype文件
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
1
奥币
627
积分
266
注册时间
2016.7.20
在线时间
52 小时

发表于 2016.9.28 23:19:05 | 显示全部楼层
这么给力,必须顶
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
5
奥币
114
积分
247
注册时间
2016.1.6
在线时间
85 小时

发表于 2016.9.30 11:11:06 | 显示全部楼层
写的好详细~
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
1
奥币
1013
积分
159
注册时间
2016.7.13
在线时间
69 小时

发表于 2016.9.30 11:12:01 | 显示全部楼层
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
奥币
254
积分
156
注册时间
2016.9.6
在线时间
88 小时

发表于 2016.10.1 20:06:12 | 显示全部楼层
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
奥币
2571
积分
986
注册时间
2016.1.15
在线时间
217 小时

发表于 2016.10.2 08:49:13 | 显示全部楼层
给力,收藏!
回复

使用道具 举报

草履虫

Rank: 2

主题
0
奥币
248
积分
2
注册时间
2016.10.2
在线时间
0 小时

发表于 2016.10.2 15:38:00 | 显示全部楼层
楼主,我也不知道哪个phenotypedata长什么样,nature文章没买看不了,求指教
回复 支持 反对

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
1
奥币
15768
积分
3330
注册时间
2016.4.25
在线时间
662 小时

发表于 2016.10.3 22:28:45 | 显示全部楼层
刚看完。收藏哟。
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
奥币
804
积分
124
注册时间
2016.6.7
在线时间
26 小时

发表于 2016.10.8 20:27:00 | 显示全部楼层
这个文章看到了,可惜生物信息小白一个,看不懂。
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

 楼主| 发表于 2016.10.9 09:47:28 | 显示全部楼层
andywang6137 发表于 2016.10.8 20:27
这个文章看到了,可惜生物信息小白一个,看不懂。

是你
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
5
奥币
1555
积分
408
注册时间
2016.6.21
在线时间
118 小时

发表于 2016.10.16 10:06:32 | 显示全部楼层
学习了
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
7
奥币
1786
积分
794
注册时间
2016.4.8
在线时间
196 小时

发表于 2016.10.24 18:57:26 | 显示全部楼层
回复

使用道具 举报

帝王蝶

Rank: 4

主题
5
奥币
193
积分
258
注册时间
2016.3.5
在线时间
60 小时

发表于 2016.10.27 15:21:16 | 显示全部楼层
这个ballgown是不是只能处理有重复的情况?
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
3
奥币
1130
积分
206
注册时间
2016.3.16
在线时间
39 小时

 楼主| 发表于 2016.11.2 22:25:10 | 显示全部楼层
zyr258 发表于 2016.10.27 15:21
这个ballgown是不是只能处理有重复的情况?

文章中的默认参数对样本量有要求,样本少的话ballgown是可选的,用DEseq或者edgeR也可以做差异表达
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表