最近在做转录组数据分析，看到了nature protocol上8月11最新发的转录组分析新流程的文章，拿出来和大家分享，自己也准备好好研究研究，希望对自己以后的分析有所帮助。文章中主要用了HISAT/HISAT2, StringTie 和 Ballgown这几个工具和软件包，都是开源免费使用的，可方便下载安装。
首先了解一下几个工具作用;
HISAT/HISAT2:比对工具了，类似于tophat2；
Stringtie：组装与定量工具。
Ballgown：为差异表达计算工具
主要流程图
具体步骤
1、创建index
首先利用下面脚本提取剪接信息（有参考GFF前提下，没有忽略此步）：
extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
然后构建HISAT2 index:
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
The --ss and --exon options（没有第一步可以不写）。indexing requires 9 GB of RAM for chromosome X and 160 GB for the whole human genome. The amount of memory is much smaller if one omits annotation information. Indexing chromosome X using one CPU core takes <10 min. It should take ~2 h to build an index for the whole human genome using eight CPU cores.
2、开始比对
各样本分别比对参考基因组
hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188104_chrX_1.fastq.gz -2
chrX_data/samples/ERR188104_chrX_2.fastq.gz -S ERR188104_chrX.sam
将SAM 转换为BAM:
samtools sort -@ 8 -o ERR188104_chrX.bam ERR188104_chrX.sam
3****、组装转录本
stringtie -p 8 -G chrX_data/genes/chrX.gtf -o
ERR188104_chrX.gtf –l ERR188104 ERR188104_chrX.bam
4、合并各个样本
$ stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o stringtie_merged.gtf
chrX_data/mergelist.txt
chrX_data/mergelist.txt：各个gtf路径放在里面。
5、估计表达丰度
stringtie –e –B -p 8 -G stringtie_merged.gtf -o
ballgown/ERR188104/ERR188104_chrX.gtf ERR188104_chrX.bam
6、加载 Ballgown R包
$ R
R version 3.2.2 (2015-08-14) -- "Fire Safety"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)
7、 加载表型数据.
An example file called geuvadis_phenodata.csv is included with the data files for this protocol (ChrX_data). In general, you will have to create this file yourself. It contains information about your RNA-seq samples, formatted as illustrated in this csv (comma-separated values) file.
pheno_data = read.csv("geuvadis_phenodata.csv")
8、加载表达丰度文件其来源于stingtie
To do this, we use the ballgown command with the following three parameters: the directory in which the data are stored (dataDir, which here is named simply ‘Ballgown’), a pattern that appears in
the sample names (samplePattern) and the phenotypic information that we loaded in the previous step (pData). Note that once a Ballgown object is created, any other Bioconductor32 package can be applied for data analysis
or data visualization.
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
9、过滤低表达基因。
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
10、鉴定差异转录本
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
11、鉴定差异基因
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
12、添加基因名字和geneID。
results_transcripts =
data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
13、按照P值从小到大排序。
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
14、保存到文件。:
write.csv(results_transcripts, "chrX_transcript_results.csv",
row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv",
row.names=FALSE)
15鉴定 q value <0.05的转录本:
subset(results_transcripts,results_transcripts
qval<0.05)
16 作图。
tropical= c('darkorange', 'dodgerblue',
'hotpink', 'limegreen', 'yellow')
palette(tropical)
17、对于基因按照FPKM 值作图 。
fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
18、对单个基因在不同样本中表达情况作图
For example, here we show how to create a plot for the 12th transcript in the data set . The first two commands below show the name of the transcript (NM_012227)
and the name of the gene that contains it (GTP binding protein 6, GTPBP6):
ballgown::transcriptNames(bg_chrX)[12]
12 "NM_012227"
ballgown::geneNames(bg_chrX)[12]
12 "GTPBP6"
>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))
</ignore_js_op>
**19、输出一个样本中一个基因座位的所有转录本的基因结构与表达丰度图 **
plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))
20、我们也可以使用plotMeans 属于一个基因的所有转录本的平均表达值。
plotMeans('MSTRG.56', bg_chrX_filt,groupvar="sex",legend=FALSE)
就这些了，总之：This protocol does not require programming expertise, but it does assume familiarity with the Unix command-line interface and the ability to run basic R scripts. Users should be comfortable running programs from the command line and editing text files in the Unix environment.
附：
软件链接
HISAT2：http://ccb.jhu.edu/software/hisat2/index.shtml
StringTie：http://ccb.jhu.edu/software/stringtie/index.shtml
Ballgown：
1. 运行R；
2. source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")
转载自基迪奥生物
