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

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