版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、一步一步教你做转录组分析(HISAT, StringTie and Ballgown)该分析流程主要根据2016年发表在Nature Protocols 上的一篇名为 Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown 的 文章撰写的,主要用到以下三个软件:HISAT ()利用大量FM 索引,以覆盖整个基因组,能够将RNA-Seq的读取与基因 组进行快速比对,相较于STAR、Tophat,该软件比对速度 快,占用内存少。StringTie()能够应用流神经网
2、络算法和可选 的de nov。组装进行转录本组装并预计表达水平。与 Cufflinks等程序相比,StringTie实现了更完整、更准确的基 因重建,并更好地预测了表达水平。Ballgown ()是R语言中 基因差异表达分析的工具,能利用RNA-Seq实验的数据 (StringTie, RSEM, Cufflinks)的结果预测基因、转录本的差异表 达。然而Ballgown并没有不能很好地检测差异外显子,而 DEXseq、rMATS和MISO可以很好解决该问题。一、数据下载Linux系统下常用的下载工具是wget,但该工 具是单线程下载,当使用它下载较大数据时比较慢,所以选 择axel,终端中
3、输入安装命令:$sudo yum install axel然后 提示输入密码获得root权限后即可自动安装,安装完成后, 输入命令axel,终端会显示如下内容,表示安装成功。 AxelZ具常用参数有:axel 选项下载目录下载地址 -s :指定每秒下载最大比特数-n :指定同时打开的线程数 -o :指定本地输出文件-S :搜索镜像并从X servers服务器下 载-N :不使用代理服务器-v :打印更多状态信息-a :打印进 度信息-h :该版本命令帮助-V :查看版本信息号#Axel安装 成功后在终端中输入命令:$axel此时在终端中会显示如下 图信息,如果不想该信息刷屏,添加参数q,采用静
4、默模式 即可。#数据下载后,进行解压:$tar-zxvfchrX_data.tar.gz解压后 利用tree命令查看数据结构,它会以树状图的形式列出目录 的内容。整个数据的结构如下图所示:chrX_gtf是X号染色体的注释文件chrX.fa是X号染色体的序 列文件indexes文件夹中是HISAT对于X号染色体的index 文件,该文件是根据序列文件chrX.fa利用hisat2-build构建 的,samples文件夹中的12个fastq文件是英格兰岛和约鲁 巴住民的X号染色体的数据。二、软件安装首先安装bioconda,它是一个自动化管理生物 信息软件的工具,安装简单,且各个软件依赖的环境
5、一同打 包且相互隔离,非常适合在服务器中搭建生信分析环境。# 下载和安装miniconda$ wget下载完成后在终端中安装 $bash Miniconda-latest-Linux-x86_64.sh 按照提示安装,完 成后$source/.bashrc #使以上的安装立即生效#输入以下 命令检验miniconda是否安装成功$ conda list显示如下图信 息说明安装成功然后利用conda install软件名+版本号安装软件即可,我们 需要安装hisat2、stringtie、samtools三个软件,安装的命令 为:$ condainstall hisat2$ condainsta
6、ll stringtie$ condainstall samtools三、分析流程1、使用HISAT将读段匹配到参考基因组上, 使用者可以提供注释文件,但HISAT依旧会检测注释文件没 有列出来的剪切位点。2、比对上的reads将会被呈递给 StringTie进行转录本组装,StringTie单独的对每个样本进行 组装,在组装的过程中顺带估算每个基因及isoform的表达 水平。3、所有的转录本都被呈递给StringTie的merge函数 进行merge,这一步是必须的,因为有些样本的转录本可能 仅仅被部分reads覆盖,无法被第二步的StringTie组装出来。 merge步骤可以创建出所有
7、样本里面都有的转录本,方便下 一步的对比。4、merge的数据再一次被呈递给StringTie, StringTie可以利用merge的数据重新估算转录本的丰度,还 能额外的提供转录本reads数量的数据给下一步的 ballgown。5、Ballgown从上一步获得所有转录本及其丰度, 根据实验条件进行分类统计。四、实战首先使用hisat2进行比对,具体用法:hisat2 options* -x (-1 -2 | -U | -sra-acc -S 主要参数:-x :参考基因组索引文件的前缀。-1 :双端测序结果的第一个文件。若有多组数据,使用逗 号将文件分隔。Reads的长度可以不一致。-2
8、:双端测序结果的第二个文件。若有多组数据,使用逗 号将文件分隔,并且文件顺序要和-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 -2 chrX_data/samples/ERR$_chrX_2.fastq.gz -SERR$_chrX.samdone将该脚本保存为1.sh,在终端中运行即 可,即:sh/脚本/所处/位置
9、/1.sh脚本执行完即可得到右图 中 12 个 sam 文件。SAM(Sequence Alignment/Map)格式是 一种通用的比对格式,用来存储reads到参考序列的比对信 息。下图是输出的比对结果,可以看到在比对样本ERR188044 时,共有1321477条reads,其中8.53%一次也未比对上, 89.68%比对上了一次,1.79%不止一次比对上,将其中8.53% 一次也未比对上的不按照顺序进行比对,发现有4.08%比对 上了一次,再将剩余的108188条reads进行单端比对,发 现50.47%未比对上,48.33%比对上了一次,1.20%比对上不止 一次,最后结果是,总共比
10、对上了 95.87%。其他的比对结果 就不一一解释了。最终我们获得了 12个sam文件: 然后通过samtools将sam文件转换为bam文件,作为 stringtie 的输入文件,具体脚本为:$ for i in ;dosamtools sort - 4 -o ERR$_chrX.bamERR$_chrX.samdone 此处 sort 默认 输出的bam文件是按其基因组位置排序,而tophat的输出 的bam文件即是按此顺序排序的。sort -n则是按reads的 ID排序。bam文件为二进制文件,占用的磁盘空间比sam 文本文件小;利用bam二进制文件的运算速度快将脚本保存 为3.sh,
11、直接在终端中执行脚本sh/脚本/所在/位置/3.sh, 最终得到的结果如下图。接下来利用stringtie对转录组进行 组装,会针对每个bam文件生成一个gtf文件,它主要记录 了转录本的组装信息,同样用一个小脚本执行该步操作: $ for i in ;dostringtie -p 8 -G ./genes/chrX.gtf -o ERR$_chrX.gtf -l ERR$ ERR$_chrX.bamDone 具体结果如下 图:然后利用软件stringtie将12个含有转录本信息的gtf文件合 并成一个gtf,此时需要预先将12个GTF文件的文件名录入 到mergelist.txt文件中,下载
12、的数据中已经给出该文件,执 行完会多出一个GTF文件,即tringtie_merged.gtf : $stringtie-merge -p 8 -G ./genes/chrX.gtf -o stringtie_merged.gtf./mergelist.txt参数-merge为转录本合并模式。在合并模式下,stringtie 将所有样品的GTF/GFF文件列表作为输入,并将这些转录本 合并/组装成非冗余的转录本集合。这种模式被用于新的差异 分析流程中,用以生成一个跨多个RNA-Seq样品的全局的、 统一的转录本。接下来,重新组装转录本并估算基因表达丰度,并为 ballgown创建读入文件。利用
13、组装好的非冗余的转录本文件 即stringtie_merged.gif和12个bam文件,执行下面的脚本 $ for i in ;dostringtie -e -B -p 8 -G stringtie_merged.gif -o ballgown/ERR$/ERR$_chrX.gtfERR$_chrX.bamdone 输出文件 在ballgown文件夹中,每个输出结果包含4个文件,如下图 接下来要用到R语言分析,选择在Windows中的Rstudio软 件中进行分析,前提是系统中已经正确安装R语言,才能使 用 Rstudio#安装需要的Rsource()biocLite(ballgown)so
14、urce()biocLite(gene filter)source()biocLite(devtools)source()biocLite(RS kittleBrewer)install.packages(dplyr)#加载要用到的语言包library(RSkittleBrewer)library(ballgown)library(genefilter) library(dplyr)library(devtools)#设置 R 语言的工作路 径setwd(F:/data/R)#读取表型数据如下图所示: read.csv(geuvadis_phenodata.csv)pheno_data#dat
15、aDir 告知数据路径,samplePattern则依据样本的名字来, pheno_data则指明了样本数据的关系,这个里面第一列样本 名需要和ballgown下面的文件夹的样本名一样,不然会报 bg_chrX= ballgown(dataDir = F:/data/R/ballgown,samplePattern = ERR, pData=pheno_data)#滤掉低丰度的基因,这里选择过滤掉样 本间差异少于一个转录本的数据 bg_chrX_filt=subset(bg_chrX,rowVars(texpr(bg_chrX)1 ,genomesubset=TRUE)#确认组间有差异的转录本
16、,在这里 我们比较male和famle之间的基因差异,指定的分析参数为 “transcripts”,主变量是“sex”,修正变量是“population,getFC 可以指定输出结果显示组间表达量的foldchange。result_transcripts=stattest(bg_chrX_filt,feature =transcript,covariate = sex,adjustvars = c(population), getFC=TRUE,meas=FPKM)#查看有差异转录本的输出结果, 如下图result_transcripts#确定各组间有差异的基因,如下图result_gene
17、s=stattest(bg_chrX_filt,feature= gene,covariate = sex,adjustvars = c(population), getFC=TRUE,meas=FPKM)#为组间有差异的转录本添加基 因名,如下图:result_transcripts=data.frame(geneNames=ballgown:geneN ames(bg_chrX_filt), geneIDs = ballgown:geneIDs(bg_chrX_filt), result_transcripts)# 按照 p-value从小到大排 result_transcripts=ar
18、range(result_transcripts,pval)result _genes=arrange(result_genes,pval)#把两个结果写入到 csv 文件中 write.csv(result_transcripts, chrX_transcript_results.csv,s=FALSE)write.csv(res ult_genes, chrX_gene_results.csv,s=FALSE)# 从以 上的输出中筛选出q值小于0.05的transcripts和genes,即 性别间表达有差异的基因 subset(result_transcr
19、ipts,result_transcripts$qval subset(result_genes,result_genes$qval#输出结果如下图: subset(result_genes,result_genes$qval#输出结果如下图:# 赋予调色板五个指定颜色tropical= c(darkorange, dodgerblue,hotpink, limegreen, yellow) 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=l og2(FPKM + 1)#查看单个转录本在样品中的分布,如图,并绘制箱线图 ballgown:transcriptNames(bg
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 2024年放射性污染防治项目合作计划书
- 2024年光纤用GECL4项目合作计划书
- 2024年气血循环机项目发展计划
- 2024年橡胶加工专用设备项目发展计划
- 吉林省吉林市第一中学2023-2024学年高二下学期7月期末考试语文试题(解析版)
- 2024复合材料和增强纤维 碳纤维增强塑料(CFRP)和金属组件十字拉伸强度的测定
- 第5课+工业革命与工厂制度+导学案 高二历史统编版(2019)选择性必修2经济与社会生活
- 北师大版小学信息技术五年级上册期中试卷含参考答案
- 有限空间作业卡
- 蜀河国土资源创建文明单位实施方案
- 人工智能在医疗领域的应用课件
- 2023年中级经济师《中级运输经济》历年真题汇编(共288题)
- 广西中小学幼儿园教师职业道德负面清单和失范行为处理办法
- VI设计教学大纲
- 【教学能力比赛】教案-水平角测量
- 智能化工程施工重难点分析
- 《阅读理解六层次》课件
- 2024届高三化学二轮复习 溶液中离子浓度的分步系数图像问题的解题思路
- 西藏安琪珠峰生物科技有限公司年产200吨高原益生菌绿色制造项目环评报告书
- 厨余垃圾餐厨垃圾堆肥系统设计方案
- 静设备吊装技术方案
评论
0/150
提交评论