该分析流程主要根据2016年发表在Nature Protocols上的一篇名为Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown的文章撰写的,主要用到以下三个软件: HISAT (http://ccb./software/hisat/index.shtml)利用大量FM索引,以覆盖整个基因组,能够将RNA-Seq的读取与基因组进行快速比对,相较于STAR、Tophat,该软件比对速度快,占用内存少。 StringTie(http://ccb./software/stringtie/)能够应用流神经网络算法和可选的de novo组装进行转录本组装并预计表达水平。与Cufflinks等程序相比,StringTie实现了更完整、更准确的基因重建,并更好地预测了表达水平。 Ballgown (https://github.com/alyssafrazee/ballgown)是R语言中基因差异表达分析的工具,能利用RNA-Seq实验的数据(StringTie, RSEM, Cufflinks)的结果预测基因、转录本的差异表达。然而Ballgown并没有不能很好地检测差异外显子,而 DEXseq、rMATS和MISO可以很好解决该问题。 一、数据下载 Linux系统下常用的下载工具是wget,但该工具是单线程下载,当使用它下载较大数据时比较慢,所以选择axel,终端中输入安装命令: $sudo yum install axel 然后提示输入密码获得root权限后即可自动安装,安装完成后,输入命令axel,终端会显示如下内容,表示安装成功。 Axel工具常用参数有: axel [选项][下载目录][下载地址] -s :指定每秒下载最大比特数 -n:指定同时打开的线程数 -o:指定本地输出文件 -S:搜索镜像并从X servers服务器下载 -N:不使用代理服务器 -v:打印更多状态信息 -a:打印进度信息 -h:该版本命令帮助 -V:查看版本信息号 #Axel安装成功后在终端中输入命令: $axel ftp://ftp.ccb./pub/RNAseq_protocol/chrX_data.tar.gz 此时在终端中会显示如下图信息,如果不想该信息刷屏,添加参数q,采用静默模式即可。 #数据下载后,进行解压: $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,它是一个自动化管理生物信息软件的工具,安装简单,且各个软件依赖的环境一同打包且相互隔离,非常适合在服务器中搭建生信分析环境。 #下载和安装miniconda $ wget https://repo./miniconda/Miniconda3-latest-Linux-x86_64.sh #下载完成后在终端中安装 $bash Miniconda-latest-Linux-x86_64.sh 按照提示安装,完成后 $source ~/.bashrc #使以上的安装立即生效 #输入以下命令检验miniconda是否安装成功 $ conda list 显示如下图信息说明安装成功 然后利用conda install 软件名+版本号安装软件即可,我们需要安装hisat2、stringtie、samtools三个软件,安装的命令为: $ condainstall hisat2 $ condainstall stringtie $ condainstall samtools 三、分析流程 1、使用HISAT将读段匹配到参考基因组上,使用者可以提供注释文件,但HISAT依旧会检测注释文件没有列出来的剪切位点。 2、比对上的reads将会被呈递给StringTie进行转录本组装,StringTie单独的对每个样本进行组装,在组装的过程中顺带估算每个基因及isoform的表达水平。 3、所有的转录本都被呈递给StringTie的merge函数进行merge,这一步是必须的,因为有些样本的转录本可能仅仅被部分reads覆盖,无法被第二步的StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本,方便下一步的对比。 4、merge的数据再一次被呈递给StringTie,StringTie可以利用merge的数据重新估算转录本的丰度,还能额外的提供转录本reads数量的数据给下一步的ballgown。 5、Ballgown从上一步获得所有转录本及其丰度,根据实验条件进行分类统计。 四、实战 首先使用hisat2进行比对,具体用法: hisat2 [options]* -x 主要参数: -x -1 -2 -S 由于该样本采用双端测序,文件数稍多,利用脚本一次性执行 $ for i in ;do hisat2 -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 -S ERR$_chrX.sam done 将该脚本保存为1.sh,在终端中运行即可,即:sh ~/脚本/所处/位置/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%比对上不止一次,最后结果是,总共比对上了95.87%。其他的比对结果就不一一解释了。 最终我们获得了12个sam文件: 然后通过samtools将sam文件转换为bam文件,作为stringtie的输入文件,具体脚本为: $ for i in ;do samtools sort -@ 4 -o ERR$_chrX.bamERR$_chrX.sam done 此处sort默认输出的bam文件是按其基因组位置排序,而tophat的输出的bam文件即是按此顺序排序的。 sort -n 则是按reads的ID排序 。 bam文件为二进制文件,占用的磁盘空间比sam文本文件小;利用bam二进制文件的运算速度快 将脚本保存为3.sh,直接在终端中执行脚本sh ~/脚本/所在/位置/3.sh,最终得到的结果如下图。 接下来利用stringtie对转录组进行组装,会针对每个bam文件生成一个gtf文件,它主要记录了转录本的组装信息,同样用一个小脚本执行该步操作: $ for i in ;do stringtie -p 8 -G ./genes/chrX.gtf -o ERR$_chrX.gtf -l ERR$ ERR$_chrX.bam Done 具体结果如下图: 然后利用软件stringtie将12个含有转录本信息的gtf文件合并成一个gtf,此时需要预先将12个GTF文件的文件名录入到mergelist.txt文件中,下载的数据中已经给出该文件,执行完会多出一个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创建读入文件。利用组装好的非冗余的转录本文件即stringtie_merged.gtf和12个bam文件,执行下面的脚本 $ for i in ;do stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR$/ERR$_chrX.gtfERR$_chrX.bam done 输出文件在ballgown文件夹中,每个输出结果包含4个文件,如下图 接下来要用到R语言分析,选择在Windows中的Rstudio软件中进行分析,前提是系统中已经正确安装R语言,才能使用Rstudio #安装需要的R包 >source('https:///biocLite.R') >biocLite('ballgown') >source('https:///biocLite.R') >biocLite('genefilter') >source('https:///biocLite.R') >biocLite('devtools') >source('https:///biocLite.R') >biocLite('RSkittleBrewer') >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<-read.csv('f: ata/r/geuvadis_phenodata.csv="">-read.csv('f:> #dataDir告知数据路径, 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) #确认组间有差异的转录本,在这里我们比较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_genes=stattest(bg_chrX_filt,feature= 'gene',covariate = 'sex',adjustvars = c('population'), getFC=TRUE,meas='FPKM') #为组间有差异的转录本添加基因名,如下图: > result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs = ballgown::geneIDs(bg_chrX_filt), result_transcripts) #按照p-value从小到大排序 >result_transcripts=arrange(result_transcripts,pval) >result_genes=arrange(result_genes,pval) #把两个结果写入到csv文件中 >write.csv(result_transcripts, 'chrX_transcript_results.csv',row.names=FALSE) >write.csv(result_genes, 'chrX_gene_results.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', '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='log2(FPKM+1)') #查看单个转录本在样品中的分布,如图,并绘制箱线图 >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)) #plotTranscripts函数可以根据指定基因的id画出在特定区段的转录本,可以通过sample函数指定看在某个样本中的表达情况,这里选用id=1750, sample=ERR188234 >plotTranscripts(ballgown::geneIDs(bg_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讨论社区 扫描二维码即可加入
|
|