参考文献:
1. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT & Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads Nature Biotechnology 2015, doi:10.1038/nbt.3122
2. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, Nature Protocols 11, 1650-1667 (2016), doi:10.1038/nprot.2016.095
之前写过基于tuxedo的RNA-seq差异表达分析工作流程。时过境迁,RNA-seq差异表达分析的pipeline也在不断进化过程中。tophat2的缺点在于效率低,无法面对新一代mapping工具的挑战。cufflinks对于复杂实验设计无能为力,还被批评作者的数学是体育老师教的(我完全不能同意这一观点,因为我就没有真正看懂过cufflinks的统计方法,所以我认为原作者非常牛)。
现在,steven小组在去年强力推出了新一代的RNA-seq差异表达的分析流程。经过我的实践,运行效率的确提高了很多。至于验证率,还有待时间的考验。不过,因为DESeq和edgeR等工具的敏感度非常高,(当然假阳性也非常高,)非常有利于编写完整的故事。所以我认为基于HISAT-StringTie和Ballgown的流程也必然会逐渐成为主流。
本文的目的不在于解释软件的原理,只是想把这个新的流程写出来,希望能对大家有所帮助。当然这个流程不包括前处理步骤,比如切除barcodes以及fastqc等。
整个流程与之前的tuxedo非常类似,HISAT2对应tophat2, StringTie对应cufflink和cuffmerge,Ballgown对应cuffdiff。在参考文献2中,作者提供了非常详细的代码,比之直接阅读软件的帮助文件,更容易让人理解。
软件的下载与安装这里就不多言了,可以去http://www.ccb.jhu.edu/software.shtml找。
第一步,使用HISAT2来mapping。
使用HISAT2,和bowtie一样,需要事前准备indexes文件。对于大多数模式生物而言,可以直接去http://www.ccb.jhu.edu/software/hisat2/index.shtml右侧提供的链接下载。如果想自己生成一个,可以参考:http://www.ccb.jhu.edu/software/hisat2/manual.shtml#the-hisat2-build-indexer的文档。
在参考文献2中,作者提供了示例代码:
$ extract_splice_sites.py gene.gtf > gene.ss $ extract_exons.py gene.gtf > gene.exon $ ## build index $ hisat2-build --ss gene.ss --exon gene.exon genome.fa prefix_tran
需要注意的是,这一步需要大量的内存,对于人类基因组而言,需要160GB的内存,约2个小时(单核)来完成计算。
下面是比对的代码:
qsub, single end
#$ -cwd #$ -o log_$JOB_ID.txt -j y #$ -S /bin/bash #$ -pe openmpi 8 #$ -t 1-4:1 #$ -M qiuworld@qiubio.com #$ -m beas #$ -r no i=$(($SGE_TASK_ID - 1)) fastq=(file1 file2 file3 file4) s=grcm38 INDEX="/path/to/genome/hsat2/ensembl_grcm38/genome_snp_tran" ChromSizeInfo="/path/to/genome/hsat2/ensembl_grcm38/ChromInfo.txt" mkdir -p fastqc fastqc -o fastqc ${fastq[$i]} sample=${fastq[$i]%%.*} mkdir -p hisat2 mkdir -p bam hisat2 -p 8 --dta \ -x $INDEX \ -U ${fastq[$i]} | samtools view -Sbho bam/hisat2.$s.$sample.bam - samtools sort bam/hisat2.$s.$sample.bam bam/hisat2.$s.$sample.srt samtools index bam/hisat2.$s.$sample.srt.bam rm bam/hisat2.$s.$sample.bam mkdir -p bedgraph bedtools genomecov -bg -split -ibam bam/hisat2.$s.$sample.srt.bam \ -g $ChromSizeInfo \ > bedgraph/hisat2.$s.$sample.bedgraph mkdir -p bigwig bedGraphToBigWig bedgraph/hisat2.$s.$sample.bedgraph \ $ChromSizeInfo \ bigwig/hisat2.$s.$sample.bw
这里有几点需要注意的,如果是SE的话,使用-U来传递fastq文件。如果是PE的话,使用-1, -2来传递。其中–dta参数会为mapping结果提供一个XS tag。这有利于下游StringTie的运算。HISAT的输出是一个SAM文件,在这里直接pipe到了一个bam文件。如果不在意文件大小的话,使用在samtools view后面加上-u参数来关闭压缩以节约时间,反正最后这个文件也是要删除的。
第二步,使用StringTie计数。
在这里,我大量抄写了文献2中的代码(应该不会有版权问题吧?)。
fastq=(file1.fq file2.fq file3.fq file4.fq) s=grcm38 GTFFILE="/path/to/genome/hsat2/ensembl_GRCm38/Mus_musculus.GRCm38.84.gtf" NUMCPUS=8 mkdir -p stringtie for ((i=0; i<=${#fastq[@]}-1; i++)); do echo [`date +"%Y-%m-%d %H:%M:%S"`] " * Assemble transcripts (StringTie)" sample=${fastq[$i]%%.*} stringtie bam/hisat2.$s.$sample.srt.bam -p $NUMCPUS \ -G ${GTFFILE} \ -o stringtie/$sample.gtf -l $sample done ## step 2 merge mkdir -p ballgown echo [`date +"%Y-%m-%d %H:%M:%S"`] "#> Merge all transcripts (StringTie)" ls -1 stringtie/*.gtf > stringtie/mergelist.txt stringtie --merge -p $NUMCPUS -G ${GTFFILE} -o \ ballgown/stringtie_merged.gtf stringtie/mergelist.txt ## step 3 estimate transcript abundance echo [`date +"%Y-%m-%d %H:%M:%S"`] "#> Estimate abundance for each sample (StringTie)" for ((i=0; i<=${#fastq[@]}-1; i++ )); do sample=${fastq[$i]%%.*} mkdir -p ballgown/$sample stringtie -e -B -p $NUMCPUS -G ballgown/stringtie_merged.gtf \ -o ballgown/$sample/$sample.gtf bam/hisat2.$s.$sample.srt.bam done echo [`date +"%Y-%m-%d %H:%M:%S"`] "#> DONE."
在step 3中,-o 参数跟着两层目录,这一点要注意,因为每个样品都会输出一些ctab文件,它们的文件名都一样,如果都放在一个目录内,会被后来生成的文件覆盖。这一点在说明文档中没有体现。
在ballgown目录中的文件会是这样
extdata/ sample01/ e2t.ctab e_data.ctab i2t.ctab i_data.ctab t_data.ctab sample02/ e2t.ctab e_data.ctab i2t.ctab i_data.ctab t_data.ctab ... sampleNN/ e2t.ctab e_data.ctab i2t.ctab i_data.ctab t_data.ctab
第三步,使用Ballgown来做差异表达分析
首先我们需要一个实验设计的表格。比如文件phenoData.tab:
sampleID group file1 WT file2 WT file3 KD file4 KD
然后中在R中运行代码:
library(ballgown) pheno_data <- read.delim("phenoData.tab") bg <- ballgown(dataDir=".", samplePattern="^file") ##必须要有这个samplePattern, 可以依照自己的文件名更改 pData(bg) = pheno_data[match(sampleNames(bg), pheno_data$sampleID), ] save(bg, file='bg.rda') ## 因为ballgown的构千时间比较长,保存一下,以后还可以再用。 library(genefilter) bg.fil <- subset(bg, "rowVars(texpr(bg)) > 1", genomesubset=TRUE) ## transcript level res_tx <- stattest(bg.fil, feature="transcript", covariate="group", getFC=TRUE, meas="FPKM") ## gene level res_genes <- stattest(bg.fil, feature="gene", covariate="group", getFC=TRUE, meas="FPKM") ## annotation anno <- data.frame(symbol=geneNames(bg.fil), geneIDs=geneIDs(bg.fil)) res_tx <- cbind(anno, res_tx) res_genes <- cbind(anno[match(res_genes$id, anno$geneIDs), ], res_genes) res_tx <- res_tx[, -c(3, 4)] res_genes <- res_genes[, -c(3, 4)] ## sort by p-value library(dplyr) res_tx <- arrange(res_tx, pval) res_genes <- arrange(res_genes, pval) ## output library(WriteXLS) WriteXLS(res_tx, "transcript_results.xls") WriteXLS(res_genes, "gene_results.xls") ## subset res_tx_sig <- subset(res_tx, res_tx$qval < 0.05) res_genes_sig <- subset(res_genes, res_genes$qval < 0.05) WriteXLS(res_tx_sig, "transcript_results_significant.xls") WriteXLS(res_genes_sig, "gene_results_significant.xls") |
这个新的差异表达分析流程是否可以支持没有重复样本的分析?比如phenoData.tab文件为:
sampleID group
file1 WT
file2 KD1
file3 KD2
file4 KD3
这里面的统计模型是否适用?
因为它背后是基于limma的,所以我认为不合适。
请教下博主,Ballgown分析后会有提供normalized的表达值吗?比如说有多个时间点。经过stringTie计数后可否直接用于DESeq2来分析?
可以的。你需要注意的是,stringTie专门有一个程序来输出用于DESeq2分析的计数.
我也想问个问题。我从文章上下的数据,用hisat2 mapping,然后用featureCounts默认参数算reads数,然后DESeq2出个差异分析,和文章里面的基本对不上,然后用edgeR之类再来一遍,和DESeq2的结果一致,和文章差异很大。然后将featureCounts的参数标准降低,比如多处mapping的都计数等等,并且参考了featureCounts的指引和这个软件的文章,感觉参数还是比较合理的,然后再来个差异分析和原文差异依旧很大。
这个测序数据是人类的,测序质量好,mapping率也高。问题只可能出在reads数上。
第一个问题是遇到这种,和原文对不上或者分析结果和实验结果对不上的时候只有换软件吗?您有没有遇到过,怎么处理的。还有就是对于这种注释的比较好的物种,你们一般用什么软件去差异分析,就是从mapping得到的东西到差异分析这中间一般用什么软件。
第二个问题是,featureCounts调整参数以后的结果和之前结果相差的都是挺大的,包括原文的结果。不知道哪个是正确的。我的意思是拿到差异分析的结果,在没有实验的情况下,能不能判断哪个是合理正确的。
谢谢。
我试用过hisat2-stringtie-ballgown.感觉它的p-value有问题。很多人抱怨这个问题。所以我暂时放弃使用hisat2.仔细阅读文档才发现,ballgown并不适合小样本量的比较。原文我找不到了。
计数的影响是挺大的,但不应该是决定因素。使用不同的软件,最后出来的结果应该差不多才是。
找到原文了,是这样写的For small sample sizes (n < 4 per group), it is often better to perform regularization. This can be done using the limma package in Bioconductor.在Nature protocols vol.11 No.9 2016, 1661页.
欧老师你好,感谢更新!
您的意思stringtie-ballgown只适合大样本比较么?那您觉得hisat2-featureCount-edgeR/DESeq2这个pipeline怎么样?如果只是小样本量比如 3-4 per group。
我记得最近有一篇文献是比较不同的pipeline的,你可以读一下: Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
欧老师你好,请问关于small rna-seq或者clip rna-seq您有经验么?
是不是不能再使用tophat或者hisat这种splice aware的aligner了?得用bwa或者bowtie了是么?
在预处理和下游处理中有没有需要注意的?谢谢!
你可以搜一下ule lab的东西。我个人的观点是看你的目的对象是pre-RNA还是mature RNA了。保险的办法就是使用bwa.