RNA-sequencing数据处理流程

farui

RNA-sequencing数据处理流程
& I5 [9 H+ l/ ^1 R2 l- U4 x( }, k
0 e6 Q5 I0 ^' Z! M% |RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。6 B8 s+ _+ a% X9 Z
, O6 L/ a# x5 @
因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。
( Z$ n  Z% R& r1 d
: \- X. A; ^8 u9 ^8 xRNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤% M* }+ Y/ ^: g2 M- g, ]+ E
第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：
( V4 w' Y3 n" e$ V8 g      library(SRAbd)
3 z  p5 T) f" y4 L" T9 Q2 s, `* k      library(DBI): O: @: _* a9 W) ]+ e$ R
      if(!file.exists('SRAmetadb.sqlite')) 5 X1 `8 h( _/ Z
        srafile <- getSRAdbFile()
/ \3 Q: b, \  n+ N3 Q      else( u% t) H7 F9 Y4 x
        srafile <- 'SRAmetadb.sqlite'; p. I. Q$ C" ]1 c5 `1 m5 I1 Z9 s+ \( T
      con = dbConnect(RSQLite::SQLite(), srafile), x7 ]7 W* z( [9 D8 a
      df_srr <- listSRAfile(sra_accn, con)                #sra_accn为已知SRA Accession号
/ I8 v! F2 `$ |) C0 y# y      l_run <- unlist(df_srr["run"])( ]7 N8 d/ F" C/ K
      getSRAfile(l_run, con, fileType = 'sra')
; n8 d4 e  s( o. E$ w( l+ g如果不知道SRA Accession号，可以通过NCBI网站搜索。
5 N1 l. ]0 n% k# C% l第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。6 c; {) R- g$ z- D8 A+ D
第三步：对应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片段到参考基因组。具体的操作如下（命令行）：
; x: I$ F3 K; l/ R( U' a. Y- d      STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
: s) b: G- h! P产生完参考基因组之后，就是对应RNA片段到参考基因组* q" |% J* w* y5 U
      STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
: E# p4 O) c' R& J这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下
5 e8 U5 |; k. p2 w- N; b' j      samtools view -bS xxxxx.sam -o xxxx.bam) A2 B. M+ Z) N' n* ~1 a
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：
0 b$ B0 _# m0 B; k" H( c+ z        1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；
3 u3 f( A" m: A% T2 q  T3 v4 Y        2） 产生一个GRangesList对象；
+ m- k' |: B# o9 K3 i) t, H% y# G" D        3） 把产生的若干.bam文件放入BamFileList对象中；
  d) _& ^4 Q7 f% {" d/ j& s. [        4） 使用summarizeOverlaps方法产生count matrix。
8 K4 g; ]: q' T" B0 I  P$ y$ ?具体操作如下：2 R1 _8 }7 j- Z5 f
      library("TxDb.Hsapiens.UCSC.hg19.knownGene")
# z* c7 J' Q6 J& Q      library("Rsamtools")
7 w8 O9 w$ D1 K8 v      library(GenomicAlignments)
1 z( F+ \& v/ e7 g      bamfiles <- BamFileList(filenames, yieldSize=2000000)
+ P. {- N' b0 I& Y  h8 m      exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
- n/ ]  y8 s1 }$ T, |, ^      flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)+ P5 r/ F; e* ]2 @" Q9 W
      param <- ScanBamParam(flag=flag)
( |- }5 T+ O$ @      CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)% I/ S3 e% ?( G- Q8 U5 o
最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。- Q7 W5 _6 `& q* q

farui

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

misser217

感谢分享，期待下一个帖子哈

yinchong42

同期待下文

lchanlon

回复 farui 的帖子% }* E" J& ]8 n$ I6 H/ q3 R
# A! @  E$ M% L% H: Q; r
Welcome! Do you have experience on how to evaluate quality of RNA-seq reads?

farui

回复 lchanlon 的帖子+ S$ L% E  L5 ^
) U# ~5 O$ k8 j! v+ x
你提醒得很好，RNA-seq的Quality control也是一个重要步骤。我目前正在学习和测试两个工具，一个是RSeQC(python)，另一个是QoRTs(java和R)，等我整理以后在发上来。