最近在做轉(zhuǎn)錄組數(shù)據(jù)分析，看到了nature protocol上8月11最新發(fā)的轉(zhuǎn)錄組分析新流程的文章，拿出來和大家分享，自己也準(zhǔn)備好好研究研究，希望對自己以后的分析有所幫助。文章中主要用了HISAT/HISAT2, StringTie 和 Ballgown這幾個工具和軟件包，都是開源免費(fèi)使用的，可方便下載安裝。
首先了解一下幾個工具作用;
HISAT/HISAT2:比對工具了，類似于tophat2；
Stringtie：組裝與定量工具。
Ballgown：為差異表達(dá)計算工具
主要流程圖
具體步驟
1、創(chuàng)建index
首先利用下面腳本提取剪接信息（有參考GFF前提下，沒有忽略此步）：
extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
然后構(gòu)建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 轉(zhuǎn)換為BAM:
samtools sort -@ 8 -o ERR188104_chrX.bam ERR188104_chrX.sam
3****、組裝轉(zhuǎn)錄本
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、估計表達(dá)豐度
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、 加載表型數(shù)據(jù).
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、加載表達(dá)豐度文件其來源于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、過濾低表達(dá)基因。
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
10、鑒定差異轉(zhuǎn)錄本
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的轉(zhuǎn)錄本:
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、對單個基因在不同樣本中表達(dá)情況作圖
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、輸出一個樣本中一個基因座位的所有轉(zhuǎn)錄本的基因結(jié)構(gòu)與表達(dá)豐度圖 **
plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))
20、我們也可以使用plotMeans 屬于一個基因的所有轉(zhuǎn)錄本的平均表達(dá)值。
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. 運(yùn)行R；
2. source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")
轉(zhuǎn)載自基迪奧生物
