一、软件介绍
该分析流程主要根据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索引，以覆盖整个基因组，能够将RNA-Seq的读取与基因组进行快速比对，相较于STAR、Tophat，该软件比对速度快，占用内存少。
StringTie(http://ccb.jhu.edu/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
Axel工具常用参数有：
axel ［选项］［下载目录］［下载地址］
-s ：指定每秒下载最大比特数
-n：指定同时打开的线程数
-o：指定本地输出文件
-S：搜索镜像并从X servers服务器下载
-N：不使用代理服务器
-v：打印更多状态信息
-a：打印进度信息
-h：该版本命令帮助
-V：查看版本信息号
#Axel安装成功后在终端中输入命令：
$ axel ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz
此时在终端中会显示如上图信息，如果不想该信息刷屏，添加参数q，采用静默模式即可。
#数据下载后，进行解压：
$ tar –zxvf chrX_data.tar.gz
chrX_gtf是X号染色体的注释文件
chrX.fa是X号染色体的序列文件
indexes文件夹中是HISAT对于X号染色体的index文件，该文件是根据序列文件chrX.fa利用hisat2-build构建的，samples文件夹中的12个fastq文件是英格兰岛和约鲁巴住民的X号染色体的数据。
三、软件安装：
首先安装bioconda，它是一个自动化管理生物信息软件的工具，安装简单，且各个软件依赖的环境一同打包且相互隔离，非常适合在服务器中搭建生信分析环境。
#下载和安装miniconda
$ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
#下载完成后在终端中安装
$ bash Miniconda-latest-Linux-x86_64.sh
#按照提示安装，完成后
$ source ~/.bashrc #使以上的安装立即生效
#输入以下命令检验miniconda是否安装成功
$ conda list
然后利用conda install 软件名+版本号安装软件即可，我们需要安装hisat2、stringtie、samtools三个软件，安装的命令为：
$ conda install hisat2
$ conda install stringtie
$ conda install 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 <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]
主要参数：
-x <hisat2-idx>：参考基因组索引文件的前缀。
-1 <m1>：双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 <m2>：双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-S <hit>：指定输出的SAM文件。
由于该样本采用双端测序，文件数稍多，利用脚本一次性执行
$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
hisat2 -p 4 -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR${i}_chrX_1.fastq.gz -2 chrX_data/samples/ERR${i}_chrX_2.fastq.gz -S ERR${i}_chrX.sam
done
将该脚本保存为1.sh，在终端中运行即可，即：
sh ~/脚本/所处/位置/1.sh
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 {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
samtools sort -@ 4 -o ERR${i}_chrX.bam ERR${i}_chrX.sam
done
此处sort默认输出的bam文件是按其基因组位置排序，而tophat的输出的bam文件即是按此顺序排序的
sort -n 则是按reads的ID排序 。
接下来利用stringtie对转录组进行组装，会针对每个bam文件生成一个gtf文件，它主要记录了转录本的组装信息，同样用一个小脚本执行该步操作：
$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
stringtie -p 8 -G ./genes/chrX.gtf -o ERR${i}_chrX.gtf -l ERR${i} ERR${i}_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 {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR${i}/ERR${i}_chrX.gtf ERR${i}_chrX.bam
done
接下来要用到R语言分析，选择在Windows中的Rstudio软件中进行分析，前提是系统中已经正确安装R语言，才能使用Rstudio
#安装需要的R包
> source("https://bioconductor.org/biocLite.R")
> biocLite("ballgown")
> source("https://bioconductor.org/biocLite.R")
> biocLite("genefilter")
> source("https://bioconductor.org/biocLite.R")
> biocLite("devtools")
> source("https://bioconductor.org/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:/data/R/geuvadis_phenodata.csv ")
#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<0.05)
> subset(result_genes,result_genes$qval<0.05)
#赋予调色板五个指定颜色
> 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)
请关注我的公众号----生信栈，不定期分享实用的生物信息干货！！！
网友评论