RNA-seq差异表达分析工作流程-2 11

参考文献:
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")

11 thoughts on “RNA-seq差异表达分析工作流程-2

  1. Reply naturehuner 3月 8,2017 1:09 上午

    这个新的差异表达分析流程是否可以支持没有重复样本的分析?比如phenoData.tab文件为:
    sampleID group
    file1 WT
    file2 KD1
    file3 KD2
    file4 KD3
    这里面的统计模型是否适用?

  2. Reply Guojian 3月 17,2017 3:08 下午

    请教下博主,Ballgown分析后会有提供normalized的表达值吗?比如说有多个时间点。经过stringTie计数后可否直接用于DESeq2来分析?

  3. Reply 生信新手啊 7月 18,2017 11:18 下午

    我也想问个问题。我从文章上下的数据,用hisat2 mapping,然后用featureCounts默认参数算reads数,然后DESeq2出个差异分析,和文章里面的基本对不上,然后用edgeR之类再来一遍,和DESeq2的结果一致,和文章差异很大。然后将featureCounts的参数标准降低,比如多处mapping的都计数等等,并且参考了featureCounts的指引和这个软件的文章,感觉参数还是比较合理的,然后再来个差异分析和原文差异依旧很大。
    这个测序数据是人类的,测序质量好,mapping率也高。问题只可能出在reads数上。
    第一个问题是遇到这种,和原文对不上或者分析结果和实验结果对不上的时候只有换软件吗?您有没有遇到过,怎么处理的。还有就是对于这种注释的比较好的物种,你们一般用什么软件去差异分析,就是从mapping得到的东西到差异分析这中间一般用什么软件。
    第二个问题是,featureCounts调整参数以后的结果和之前结果相差的都是挺大的,包括原文的结果。不知道哪个是正确的。我的意思是拿到差异分析的结果,在没有实验的情况下,能不能判断哪个是合理正确的。
    谢谢。

    • Reply admin 7月 19,2017 10:37 上午

      我试用过hisat2-stringtie-ballgown.感觉它的p-value有问题。很多人抱怨这个问题。所以我暂时放弃使用hisat2.仔细阅读文档才发现,ballgown并不适合小样本量的比较。原文我找不到了。
      计数的影响是挺大的,但不应该是决定因素。使用不同的软件,最后出来的结果应该差不多才是。

      • Reply admin 7月 20,2017 10:31 上午

        找到原文了,是这样写的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页.

  4. Reply super 11月 29,2017 9:44 上午

    欧老师你好,感谢更新!
    您的意思stringtie-ballgown只适合大样本比较么?那您觉得hisat2-featureCount-edgeR/DESeq2这个pipeline怎么样?如果只是小样本量比如 3-4 per group。

    • Reply admin 11月 30,2017 2:11 下午

      我记得最近有一篇文献是比较不同的pipeline的,你可以读一下: Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis

  5. Reply super 11月 29,2017 9:46 上午

    欧老师你好,请问关于small rna-seq或者clip rna-seq您有经验么?
    是不是不能再使用tophat或者hisat这种splice aware的aligner了?得用bwa或者bowtie了是么?
    在预处理和下游处理中有没有需要注意的?谢谢!

Leave a Reply

  

  

  

%d 博主赞过: