一步一步教你做转录组分析HISATStringTieandBallgownWord格式.docx
《一步一步教你做转录组分析HISATStringTieandBallgownWord格式.docx》由会员分享，可在线阅读，更多相关《一步一步教你做转录组分析HISATStringTieandBallgownWord格式.docx（6页珍藏版）》请在冰豆网上搜索。
打印更多状态信息-a：
打印进度信息-h：
该版本命令帮助-V：
查看版本信息号#Axel安装成功后在终端中输入命令：
$axel此时在终端中会显示如下图信息，如果不想该信息刷屏，添加参数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下载完成后在终端中安装$bashMiniconda-latest-Linux-x86_64.sh按照提示安装，完成后$source~/.bashrc#使以上的安装立即生效#输入以下命令检验miniconda是否安装成功$condalist显示如下图信息说明安装成功
然后利用condainstall软件名+版本号安装软件即可，我们需要安装hisat2、stringtie、samtools三个软件，安装的命令为：
$condainstallhisat2$condainstallstringtie$condainstallsamtools
三、分析流程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{-1-2|-U|–sra-acc}[-S]主要参数：
-x：
参考基因组索引文件的前缀。
-1：
双端测序结果的第一个文件。
若有多组数据，使用逗号将文件分隔。
Reads的长度可以不一致。
-2：
双端测序结果的第二个文件。
若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。
-S：
指定输出的SAM文件。
由于该样本采用双端测序，文件数稍多，利用脚本一次性执行$foriin;
dohisat2-p4-xchrX_data/indexes/chrX_tran-1chrX_data/samples/ERR$_chrX_1.fastq.gz-2chrX_data/samples/ERR$_chrX_2.fastq.gz-SERR$_chrX.samdone将该脚本保存为1.sh，在终端中运行即可，即：
sh~/脚本/所处/位置/1.sh脚本执行完即可得到右图中12个sam文件。
SAM（SequenceAlignment/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的输入文件，具体脚本为：
$foriin;
dosamtoolssort-@4-oERR$_chrX.bamERR$_chrX.samdone此处sort默认输出的bam文件是按其基因组位置排序，而tophat的输出的bam文件即是按此顺序排序的。
sort-n则是按reads的ID排序。
bam文件为二进制文件，占用的磁盘空间比sam文本文件小；
利用bam二进制文件的运算速度快将脚本保存为3.sh，直接在终端中执行脚本sh~/脚本/所在/位置/3.sh，最终得到的结果如下图。
接下来利用stringtie对转录组进行组装，会针对每个bam文件生成一个gtf文件，它主要记录了转录本的组装信息，同样用一个小脚本执行该步操作：
dostringtie-p8-G./genes/chrX.gtf-oERR$_chrX.gtf-lERR$ERR$_chrX.bamDone具体结果如下图：
然后利用软件stringtie将12个含有转录本信息的gtf文件合并成一个gtf，此时需要预先将12个GTF文件的文件名录入到mergelist.txt文件中，下载的数据中已经给出该文件，执行完会多出一个GTF文件，即tringtie_merged.gtf：
$stringtie--merge-p8-G./genes/chrX.gtf-ostringtie_merged.gtf./mergelist.txt
参数--merge为转录本合并模式。
在合并模式下，stringtie将所有样品的GTF/GFF文件列表作为输入，并将这些转录本合并/组装成非冗余的转录本集合。
这种模式被用于新的差异分析流程中，用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。
接下来，重新组装转录本并估算基因表达丰度，并为ballgown创建读入文件。
利用组装好的非冗余的转录本文件即stringtie_merged.gtf和12个bam文件，执行下面的脚本$foriin;
dostringtie-e-B-p8-Gstringtie_merged.gtf-oballgown/ERR$/ERR$_chrX.gtfERR$_chrX.bamdone输出文件在ballgown文件夹中，每个输出结果包含4个文件，如下图
接下来要用到R语言分析，选择在Windows中的Rstudio软件中进行分析，前提是系统中已经正确安装R语言，才能使用Rstudio
#安装需要的R包>
source（'
'
）>
biocLite（'
ballgown'
genefilter'
devtools'
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#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'
）#为组间有差异的转录本添加基因名，如下图：
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#输出结果如下图：
#赋予调色板五个指定颜色>
tropical=c（'
darkorange'
'
dodgerblue'
'
hotpink'
limegreen'
yellow'
palette（tropical）#以FPKM为参考值作图，以性别作为区分条件>
fpkm=texpr（bg_chrX,meas='
）#方便作图将其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]>
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='
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（'
GeneMSTRG.575insampleERR188234'
）,sample=c（'
ERR188234'
））#这里以id=575的基因为例（对应上一步作图）>
plotMeans（'
MSTRG.575'
bg_chrX_filt,groupvar='
legend=FALSE）
广告：
转录组分析视频教程，欢迎登陆网易云课堂观看：
扫码获取如有任何问题欢迎大家加入WeGAP讨论社区扫描二维码即可加入
升级会员 