RNA-sequencing数据处理流程

farui

RNA-sequencing数据处理流程
% d; p( _8 b+ ]: Q
" n, ]5 l0 ]  s. J1 E5 hRNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。
8 k" y, P, I; A3 j/ ]- {* h3 g- m" D5 z5 j6 `: @
因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。
/ |, D' @. \4 Q* K/ w
$ M! u+ v# e0 ]+ Q) \7 @/ J/ k# bRNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤
4 V% r/ E8 H4 c- Q" m第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：: I, M1 o; F: i( E
      library(SRAbd)0 h- G8 p2 f6 P
      library(DBI)
8 V) F1 r/ M" c. v" I3 R4 w) t      if(!file.exists('SRAmetadb.sqlite')) 7 ?+ e: g" i! u" q, L! i' z
        srafile <- getSRAdbFile()% {" T6 C6 k! p; x
      else
4 L7 \( n. R4 y) J        srafile <- 'SRAmetadb.sqlite'
7 K5 j5 P' F. q0 z; Z5 H      con = dbConnect(RSQLite::SQLite(), srafile)4 V8 k6 n) B5 v, r0 }: ~
      df_srr <- listSRAfile(sra_accn, con)                #sra_accn为已知SRA Accession号& T- N- t9 j0 z# u% `6 k
      l_run <- unlist(df_srr["run"])2 G/ i8 a& L, D. _: U
      getSRAfile(l_run, con, fileType = 'sra')3 v/ n; c/ U$ x! f7 j; V+ B
如果不知道SRA Accession号，可以通过NCBI网站搜索。
0 m* T2 K7 E$ E' L- U4 M: Y+ O第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。
% ]6 w, {& [- d4 j0 l第三步：对应RNA片段到参考基因组。这一步，有很多工具可以实现，比如BWA，Bowtie/Bowtie2，GSNAP，TopHat2，或是STAR。这里我介绍STAR工具。STAR工具包可以在https://github.com/alexdobin/STAR/releases下载。同时，我们还需要下载genome fasta sequence文件和annotation GTF文件来产生参考基因组。以人类homo sapiens为例，我们从ENSEMBL（ftp://ftp.ensembl.org/pub/release-81/......）下载FASTA和GTF文件。使用STAR命令两次但是不同的参数来分别产生参考基因组和对应RNA片段到参考基因组。具体的操作如下（命令行）：
2 A$ d' x& {7 S5 V" ^9 |" q      STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
& G; S$ j0 H4 b: a产生完参考基因组之后，就是对应RNA片段到参考基因组
+ J8 e5 Y) E9 l- Z- O8 y      STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
. G5 e+ f) R0 V: c+ i' e. C* p这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下$ c* x4 X; h+ k+ w2 l9 Q
      samtools view -bS xxxxx.sam -o xxxx.bam- @8 i1 f0 o+ i2 o2 ~2 x
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：
% D# G9 y2 Y: o; [5 x/ S        1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；$ Y, L9 w; T' _( |
        2） 产生一个GRangesList对象；
( s9 n4 |! h" Z/ v3 C        3） 把产生的若干.bam文件放入BamFileList对象中；8 F" n6 d) }4 Z
        4） 使用summarizeOverlaps方法产生count matrix。: v% f3 _4 I2 `% T
具体操作如下：9 F9 H0 c4 L# K( o
      library("TxDb.Hsapiens.UCSC.hg19.knownGene")8 Z3 {5 P+ Z+ d, X3 ^5 k
      library("Rsamtools")
% G. {, {% U7 q& W6 A      library(GenomicAlignments)
0 ~0 L- I: ]2 y5 r) w! o. |# b% ~      bamfiles <- BamFileList(filenames, yieldSize=2000000)7 f) o* T8 y( U' s9 `# S
      exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
4 G+ d" _. ]0 e      flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
4 I& Z$ G; @6 [' o      param <- ScanBamParam(flag=flag)' d' j* o/ i- c- T  g. r
      CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)- h- l0 J( _# G+ f2 Y) R
最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。. k$ k' l" ?  d4 @$ E. B( f: e7 V

farui

如果有朋友感兴趣differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，请跟帖，我可以分享一些我的经验。

farui

回复 lchanlon 的帖子
+ C0 h, ]" a4 ?  v; V
0 w* Q- D- J9 _, V' |你提醒得很好，RNA-seq的Quality control也是一个重要步骤。我目前正在学习和测试两个工具，一个是RSeQC(python)，另一个是QoRTs(java和R)，等我整理以后在发上来。