一步一步教你做转录组分析HISATStringTieandBallgown.docx

1、一步一步教你做转录组分析HISATStringTieandBallgown步一步教你做转录组分析( HISAT, StringTie andBallgown )该分析流程主要根据 2016 年发表在 Nature Protocols 上的一篇名为 Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown 的文章撰写的，主要用到以下三个软件： HISAT (http:/ccb.jhu.edu/software/hisat/index.shtml) 利用大量 FM 索

2、引，以覆盖整个基因组，能够将 RNA-Seq 的读取与基因 组进行快速比对，相较于 STAR 、Tophat ，该软件比对速度 快，占用内存少。 StringTie(http:/ccb.jhu.edu/software/stringtie/) 能够应用流 神经网络算法和可选的 de novo 组装进行转录本组装并预计 表达水平。与 Cufflinks 等程序相比， StringTie 实现了更完 整、更准确的基因重建， 并更好地预测了表达水平。 Ballgown ( 是 R 语言中基因 差异表达分析的工具，能利用 RNA-Seq 实验的数据 (StringTie, RSEM, Cufflink

3、s) 的结果预测基因、转录本的差 异表达。然而 Ballgown 并没有不能很好地检测差异外显子， 而 DEXseq 、 rMATS 和 MISO 可以很好解决该问题。一、 数据下载 Linux 系统下常用的下载工具是 wget ，但该工 具是单线程下载，当使用它下载较大数据时比较慢，所以选 择 axel ，终端中输入安装命令： $sudo yum install axel 然后 提示输入密码获得 root 权限后即可自动安装，安装完成后， 输入命令 axel ，终端会显示如下内容，表示安装成功。Axel 工具常用参数有： axel 选项下载目录 下载地 址-s :指定每秒下载最大比特数-n

4、:指定同时打开的线程 数-o :指定本地输出文件-S :搜索镜像并从 X servers服务 器下载-N :不使用代理服务器-v :打印更多状态信息-a :打 印进度信息-h :该版本命令帮助-V :查看版本信息号#Axel 安装成功后在终端中输入命令: $axelftp:/ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.g z 此时在终端中会显示如下图信息，如果不想该信息刷屏， 添加参数q，采用静默模式即可。#数据下载后，进行解压：$tar之xvfchrX_data.tar.gz 解压后 利用 tree 命令查看数据结构， 它会以树状图的形

5、式列出目录 的内容。整个数据的结构如下图所示:chrX_gtf是X号染色体的注释文件 chrX.fa是X号染色体的 序列文件 indexes 文件夹中是 HISAT 对于 X 号染色体的 index 文件，该文件是根据序列文件 chrX.fa 利用 hisat2-build 构建的， samples 文件夹中的 12 个 fastq 文件是英格兰岛和 约鲁巴住民的 X 号染色体的数据。二、软件安装首先安装 bioconda ，它是一个自动化管理生物 信息软件的工具，安装简单，且各个软件依赖的环境一同打 包且相互隔离，非常适合在服务器中搭建生信分析环境。 下载和安装 miniconda$ wge

6、t https:/repo.continuum.io/miniconda/Miniconda3-latest-Lin ux-x86_64.sh# 下载完成后在终端中安装 $bash Miniconda-latest-Linux-x86_64.sh 按照提示安装，完成后 $source /.bashrc # 使以上的安装立即生效 #输入以下命令 检验 miniconda 是否安装成功 $ conda list 显示如下图信息说 明安装成功 然后利用 conda install 软件名 + 版本号安装软件即可，我们 需要安装 hisat2 、 stringtie 、 samtools 三个软件，安

7、装的命 令为： $ condainstall hisat2$ condainstall stringtie$ condainstall samtools三、分析流程 1、使用 HISAT 将读段匹配到参考基因组上 , 使用者可以提供注释文件， 但 HISAT 依旧会检测注释文件没 有列出来的剪切位点。 2 、比对上的 reads 将会被呈递给 StringTie 进行转录本组装， StringTie 单独的对每个样本进 行组装，在组装的过程中顺带估算每个基因及 isoform 的表 达水平。 3 、所有的转录本都被呈递给 StringTie 的 merge 函数进行 merge ，这一步是必须

8、的，因为有些样本的转录本 可能仅仅被部分 reads 覆盖，无法被第二步的 StringTie 组 装出来。 merge 步骤可以创建出所有样本里面都有的转录本， 方便下一步的对比。 4、 merge 的数据再一次被呈递给StringTie ， StringTie 可以利用 merge 的数据重新估算转录 本的丰度，还能额外的提供转录本 reads 数量的数据给下一 步的 ballgown 。 5、Ballgown 从上一步获得所有转录本及其 丰度，根据实验条件进行分类统计。四、实战首先使用 hisat2 进行比对，具体用法： hisat2 options* -x-1 -2 卜U | sa-a

9、cc -S 主要参数：-x :参考基因组 索引文件的前缀。-1 ：双端测序结果的第一个文件。 若有多组数据， 使用逗号 将文件分隔。 Reads 的长度可以不一致。-2 ：双端测序结果的第二个文件。 若有多组数据， 使用逗号 将文件分隔， 并且文件顺序要和 -1 参数对应。 Reads 的长度 可以不一致。-S ：指定输出的 SAM 文件。由于该样本采用双端测序，文 件数稍多，利用脚本一次性执行 $ for i in ;dohisat2 -p 4 -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR$_chrX_1.fastq.gz -

10、2 chrX_data/samples/ERR$_chrX_2.fastq.gz -S ERR$_chrX.samdone 将该脚本保存为 1.sh ，在终端中运行 即可，即： sh /脚本/所处/位置/1 .sh 脚本执行完即可得到右 图中 12 个 sam 文件。SAM(Sequence Alignment/Map) 格式 是一种通用的比对格式，用来存储 reads 到参考序列的比对 信息。下图是输出的比对结果，可以看到在比对样本ERR188044 时，共有 1321477 条 reads ，其中 8.53% 一次 也未比对上， 89.68% 比对上了一次， 1.79% 不止一次比对上，

11、将其中 8.53% 一次也未比对上的不按照顺序进行比对， 发现 有 4.08% 比对上了一次，再将剩余的 108188 条 reads 进行 单端比对，发现 50.47% 未比对上， 48.33% 比对上了一次， 1.20% 比对上不止一次， 最后结果是， 总共比对上了 95.87% 。 其他的比对结果就不一一解释了。 最终我们获得了 12 个 sam 文件：然后通过 samtools 将 sam 文件转换为 bam 文件，作为 stringtie 的输入文件， 具体脚本为： $ for i in ;dosamtools sort - 4 -o ERR$_chrX.bamERR$_chrX.s

12、amdone 此处 sort 默 认输出的 bam 文件是按其基因组位置排序，而 tophat 的输 出的 bam 文件即是按此顺序排序的。 sort -n 则是按 reads 的ID排序bam文件为二进制文件，占用的磁盘空间比sam 文本文件小；利用 bam 二进制文件的运算速度快将脚本保存 为 3.sh ，直接在终端中执行脚本 sh / 脚本 /所在 /位置 /3.sh ，最终得到的结果如下图。接下来利用 stringtie 对转录组进行 组装，会针对每个 bam 文件生成一个 gtf 文件，它主要记录 了转录本的组装信息，同样用一个小脚本执行该步操作： $ for i in ;dostr

13、ingtie -p 8 -G ./genes/chrX.gtf -o ERR$_chrX.gtf -l ERR$ ERR$_chrX.bamDone 具体结果如 下图： 然后利用软件 stringtie 将 12 个含有转录本信息的 gtf 文件合 并成一个 gtf ，此时需要预先将 12 个 GTF 文件的文件名录 入到 mergelist.txt 文件中，下载的数据中已经给出该文件， 执行完会多出一个 GTF 文件，即 tringtie_merged.gtf ： $stringtie-merge -p 8 -G ./genes/chrX.gtf -o stringtie_merged.gt

14、f./mergelist.txt 参数 -merge 为转录本合并模式。 在合并模式下， stringtie 将所有样品的 GTF/GFF 文件列表作为输入，并将这些转录 本合并 / 组装成非冗余的转录本集合。 这种模式被用于新的差 异分析流程中，用以生成一个跨多个 RNA-Seq 样品的全局 的、统一的转录本。 接下来，重新组装转录本并估算基因表达丰度，并为 ballgown 创建读入文件。 利用组装好的非冗余的转录本文件 即 stringtie_merged.gtf 和 12 个 bam 文件，执行下面的脚本 $ for i in ;dostringtie -e -B -p 8 -G st

15、ringtie_merged.gtf -o ballgown/ERR$/ERR$_chrX.gtfERR$_chrX.bamdone 输出 文件在 ballgown 文件夹中，每个输出结果包含 4 个文件， 如下图 接下来要用到 R 语言分析，选择在 Windows 中的 Rstudio 软件中进行分析， 前提是系统中已经正确安装 R 语言， 才能 使用 Rstudio#安装需要的 R包 source(https:/bioconductor.org/biocLite.R)biocLite(b allgown)source(https:/bioconductor.org/biocLite.R)b

16、io cLite(genefilter)source(https:/bioconductor.org/biocLite.R)biocLite(devtools)source(https:/bioconductor.org/bi ocLite.R)biocLite(RSkittleBrewer)install.packages(dply r)#加载要用到的语言包 library(RSkittleBrewer)library(ballgown)library(genefilter)library(dplyr)library(devtools)# 设置 R 语言的工作路 径setwd(F:/data

17、/R)# 读取表型数据如下图所示： read.csv(geuvadis_phenodata.csv)pheno_data#dataDi r 告知数据路径， samplePattern 则依据样本的名字来， pheno_data 则指明了样本数据的关系，这个里面第一列样 本名需要和 ballgown 下面的文件夹的样本名一样， 不然会报 错 bg_chrX= ballgown(dataDir =“F:/data/R/ballgown,samplePattern = ERR, pData=pheno_data)# 滤掉低丰度的基因， 这里选择过滤掉样 本间差异少于一个转录本的数据 bg_chrX_

18、filt=subset(bg_chrX,rowVars(texpr(bg_chrX) 1,genomesubset=TRUE)# 确认组间有差异的转录本， 在这 里我们比较 male 和 famle 之间的基因差异，指定的分析参 数为“ transcripts 主变”量，是“ sex ”修，正变量是“ population getFC 可以指定输出结果显示组间表达量的 foldchange 。 result_transcripts=stattest(bg_chrX_filt,featu re = transcript,covariate = sex,adjustvars = c(populat

19、ion), getFC=TRUE,meas=FPKM)# 查看有差异转录本的输出结 果，如下图 result_transcripts# 确定各组间有差异的基因，如下图 result_genes=stattest(bg_chrX_filt,feature= gene,covariate = sex,adjustvars = c(population), getFC=TRUE,meas=FPKM)# 为组间有差异的转录本添加 基因名，如下图： result_transcripts=data.frame(geneNames=ballgown:gen eNames(bg_chrX_filt), gen

20、eIDs = ballgown:geneIDs(bg_chrX_filt), result_transcripts)# 按照 p-value 从小到大排序 result_transcripts=arrange(result_transcripts,pval)res ult_genes=arrange(result_genes,pval)# 把两个结果写入到 csv 文件中 write.csv(result_transcripts,chrX_transcript_results.csv,row.names=FALSE)write.cs v(result_genes,chrX_gene_resul

21、ts.csv,row.names=FALSE)# 从以上的输出中筛选出 q 值小于 0.05 的 transcripts 和 genes ，即性别 间表达有差异的基因 subset(result_transcripts,result_transcripts$qval subset(result_genes,result_genes$qval# 输出结果如下图： subset(result_genes,result_genes$qval# 输出结果如下图： #赋予调色板五个指定颜色 tropical= c(darkorange, dodgerblue,hotpink, limegreen, ye

22、llow) palette(tropical)# 以 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) #查看单个转录本在样品中的分布，如图，并绘制箱线图 ballgown:transcriptNames(bg_chrX)12ballgown:gene Names(bg_chrX)12pl

23、ot(fpkm12, 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(fpkm12, jitter(as.numeric(pheno_data$sex), col=as.numeric(pheno_data$sex)#plotTranscripts 函数可以根据指定基因的 id 画出在特定区 段的转录本，可以通过 sample 函数指定

24、看在某个样本中的 表达情况，这里选用 id=1750, sample=ERR188234plotTranscripts(ballgown:geneIDs(b g_chrX)1750,bg_chrX, main=c(Gene MSTRG.575 in sample ERR188234), sample=c(ERR188234)# 这里以id=575 的基因为例(对应上一步作 图) plotMeans(MSTRG.575,bg_chrX_filt,groupvar=sex, legend=FALSE) 广告：转录组分析视频教程，欢迎登陆网易云课堂观看：扫 码获取如有任何问题欢迎大家加入 WeGAP 讨论社区扫描二 维码即可加入