⼀步⼀步教你做转录组分析（HISAT,StringTieandBallgown）该分析流程主要根据2016年发表在Nature Protocols上的⼀篇名为Transcript-level expression
analysis of RNA-q experiments with HISAT, StringTie and Ballgown的⽂章撰写的，主要⽤
到以下三个软件：
HISAT (ccb.jhu.edu/software/hisat/index.shtml)利⽤⼤量FM索引，以覆盖整个基因组，能
够将RNA-Seq的读取与基因组进⾏快速⽐对，相较于STAR、Tophat，该软件⽐对速度快，占
⽤内存少。
StringTie(ccb.jhu.edu/software/stringtie/)能够应⽤流神经⽹络算法和可选的de novo组装进
⾏转录本组装并预计表达⽔平。与Cufflinks等程序相⽐，StringTie实现了更完整、更准确的基
因重建，并更好地预测了表达⽔平。
Ballgown (/alyssafrazee/ballgown)是R语⾔中基因差异表达分析的⼯具，能利
⽤RNA-Seq实验的数据(StringTie, RSEM, Cufflinks)的结果预测基因、转录本的差异表达。然⽽
Ballgown并没有不能很好地检测差异外显⼦，⽽ DEXq、rMATS和MISO可以很好解决该问
题。
⼀、数据下载
Linux系统下常⽤的下载⼯具是wget，但该⼯具是单线程下载，当使⽤它下载较⼤数据时⽐较
慢，所以选择axel，终端中输⼊安装命令：
$sudo yum install axel
然后提⽰输⼊密码获得root权限后即可⾃动安装，安装完成后，输⼊命令axel，终端会显⽰如下
请说话内容，表⽰安装成功。
Axel⼯具常⽤参数有：
希腊语axel ［选项］［下载⽬录］［下载地址］
-s ：指定每秒下载最⼤⽐特数
-
n：指定同时打开的线程数
-o：指定本地输出⽂件
-S：搜索镜像并从X rvers服务器下载
-N：不使⽤代理服务器
奥林匹克圣歌-v：打印更多状态信息
-a：打印进度信息
-h：该版本命令帮助
-V：查看版本信息号
#Axel安装成功后在终端中输⼊命令：
$axel ftp://ftpb.jhu.edu/pub/RNAq_protocol/chrX_
此时在终端中会显⽰如下图信息，如果不想该信息刷屏，添加参数q，采⽤静默模式即可。
#数据下载后，进⾏解压：
$tar–zxvfchrX_
解压后利⽤tree命令查看数据结构，它会以树状图的形式列出⽬录的内容。整个数据的结构如下图所⽰：
chrX_gtf是X号染⾊体的注释⽂件
chrX.fa是X号染⾊体的序列⽂件
indexes⽂件夹中是HISAT对于X号染⾊体的index⽂件，该⽂件是根据序列⽂件chrX.fa利⽤hisat2-build构建的，samples⽂件夹中的12个fastq⽂件是英格兰岛和约鲁巴住民的X号染⾊体的数据。
⼆、软件安装
⾸先安装bioconda，它是⼀个⾃动化管理⽣物信息软件的⼯具，安装简单，且各个软件依赖的环境⼀同打包且相互隔离，⾮常适合在服务器中搭建⽣信分析环境。ratatouille
#下载和安装miniconda
$ inuum.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 三个软件，安装的命令为：
$ 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从上⼀步获得所有转录本及其丰度，根据实验条件进⾏分类统计。
ong四、实战
⾸先使⽤hisat2进⾏⽐对，具体⽤法：
hisat2 [options]* -x {-1 -2 | -U | –sra-acc } [-S ]
主要参数：
-x ：参考基因组索引⽂件的前缀。
-1 ：双端测序结果的第⼀个⽂件。若有多组数据，使⽤逗号将⽂件分隔。Reads的长度可以不⼀致。
-2 ：双端测序结果的第⼆个⽂件。若有多组数据，使⽤逗号将⽂件分隔，并且⽂件顺序要和-1参数对应。Reads的长度可以不⼀致。
-S ：指定输出的SAM⽂件。
由于该样本采⽤双端测序，⽂件数稍多，利⽤脚本⼀次性执⾏
$ for i in ;do
hisat2 -p 4 -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR$_chrX_ -2 chrX_data/samples/ERR$_chrX_ -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 ./f -o ERR$_f -l ERR$ ERR$_chrX.bam
Done
具体结果如下图：
然后利⽤软件stringtie将12个含有转录本信息的gtf⽂件合并成⼀个gtf，此时需要预先将12个GTF ⽂件的⽂件名录⼊到⽂件中，下载的数据中已经给出该⽂件，执⾏完会多出⼀个GTF⽂件，即f：
$stringtie--merge -p 8 -G ./f -o f./
参数--merge为转录本合并模式。在合并模式下，stringtie将所有样品的GTF/GFF⽂件列表作为输⼊，并将这些转录本合并/组装成⾮冗余的转录本集合。这种模式被⽤于新的差异分析流程中，⽤以⽣成⼀个跨多个RNA-Seq样品的全局的、统⼀的转录本。
接下来，重新组装转录本并估算基因表达丰度，并为ballgown创建读⼊⽂件。利⽤组装好的⾮冗余的转录本⽂件即f和12个bam⽂件，执⾏下⾯的脚本
$ for i in ;do
stringtie -e -B -p 8 -G f -o ballgown/ERR$/ERR$_fERR$_chrX.bam doneshmily是什么意思
输出⽂件在ballgown⽂件夹中，每个输出结果包含4个⽂件，如下图
接下来要⽤到R语⾔分析，选择在Windows中的Rstudio软件中进⾏分析，前提是系统中已经正确安装R语⾔，才能使⽤Rstudio
#安装需要的R包
>source('bioconductor/biocLite.R')
voodoo
>biocLite('ballgown')
>source('bioconductor/biocLite.R')
>biocLite('genefilter')
bass是什么意思>source('bioconductor/biocLite.R')
>biocLite('devtools')
>source('bioconductor/biocLite.R')
>biocLite('RSkittleBrewer')
>install.packages('dplyr')
#加载要⽤到的语⾔包
> library(RSkittleBrewer)
>library(ballgown)
>library(genefilter)
>library(dplyr)
>library(devtools)
#设置R语⾔的⼯作路径
>twd('F:/data/R')
#读取表型数据如下图所⽰：
>read.csv('geuvadis_phenodata.csv')
>pheno_data<-read.csv('f: ata/r/geuvadis_phenodata.csv="">
#dataDir告知数据路径， samplePattern则依据样本的名字来，pheno_data则指明了样本数据的关系，这个⾥⾯第⼀列样本名需要和ballgown下⾯的⽂件夹的样本名⼀样，不然会报错
>bg_chrX= ballgown(dataDir = “F:/data/R/ballgown',samplePattern = 'ERR',
pData=pheno_data)
#滤掉低丰度的基因，这⾥选择过滤掉样本间差异少于⼀个转录本的数据tast是什么意思
>bg_chrX_filt=subt(bg_chrX,'rowVars(texpr(bg_chrX))>1',genomesubt=TRUE)
#确认组间有差异的转录本，在这⾥我们⽐较male和famle之间的基因差异，指定的分析参数
为“transcripts”，主变量是“x”，修正变量是“population”，getFC可以指定输出结果显⽰组间表达量的foldchange。
>result_transcripts=stattest(bg_chrX_filt,feature = 'transcript',covariate = 'x',adjustvars =
c('population'), getFC=TRUE,meas='FPKM')
#查看有差异转录本的输出结果，如下图
> result_transcripts
#确定各组间有差异的基因，如下图
