RNA-sequencing数据处理流程

farui

RNA-sequencing数据处理流程
1 X7 o; n/ n; ]" m, Y& k4 \3 s- L  n, Q8 a
RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。
6 S# u, _5 H3 \9 q# O* ~  Q2 E+ _, i" w' A- j' G; I$ u
因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。# ]/ s+ ^# {0 m5 |
: K7 u0 \# \0 B: R2 c' o) R, B
RNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤
  q$ S$ ?3 O7 U& G第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：
, v% c0 A. X& C, k      library(SRAbd)) U! U+ D5 v9 U/ D
      library(DBI)
1 x( O. f/ C3 w, d% X0 n      if(!file.exists('SRAmetadb.sqlite'))
. q, l  T/ |5 p( ~! e; A* ]        srafile <- getSRAdbFile()7 @' D" Z4 G4 d5 d4 F9 N! I) u
      else/ @7 h& [6 C% Y, \
        srafile <- 'SRAmetadb.sqlite'. r% C2 y7 Z# ^: V8 @  }' a$ |
      con = dbConnect(RSQLite::SQLite(), srafile)! h6 |9 \/ f# h0 o3 ^& u  N- l
      df_srr <- listSRAfile(sra_accn, con)                #sra_accn为已知SRA Accession号
( d, z/ i# p5 J4 W      l_run <- unlist(df_srr["run"])& R( j: O0 O# D2 \' ?5 m8 j% [
      getSRAfile(l_run, con, fileType = 'sra')
3 S$ j1 ]8 E* d2 O/ ^如果不知道SRA Accession号，可以通过NCBI网站搜索。
+ Q& F  y1 X$ d! H( ^7 h第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。
: ]6 d9 E' F, w1 `1 J, G+ }第三步：对应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 I# `$ b  G2 v/ K
      STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径2 C# z- b( v+ |6 I* b9 R
产生完参考基因组之后，就是对应RNA片段到参考基因组
0 I- W. |1 p# D9 Z0 H      STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
( ~  L/ p6 M% P# b/ s: ]这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下/ T4 i4 K& {! f2 Y  f+ H
      samtools view -bS xxxxx.sam -o xxxx.bam
5 [& _. T: {6 ]. M3 m1 B* U* o- q第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：
8 G& A9 S# {, S2 X. ^+ _/ ~0 t/ ^        1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；
5 g: [) b  y, H8 I& g        2） 产生一个GRangesList对象；0 P" a' z4 v7 S& H
        3） 把产生的若干.bam文件放入BamFileList对象中；3 v% G* a9 x2 f9 |0 i: _, i. _5 O
        4） 使用summarizeOverlaps方法产生count matrix。
; o5 F, R4 R, R7 `9 ^4 N6 ~) }0 C具体操作如下：/ d+ i" Q; A/ I# j' O! L: [
      library("TxDb.Hsapiens.UCSC.hg19.knownGene")6 G2 ~# l( R' J+ T3 s: U
      library("Rsamtools")  W$ I* r" a# H; r" l* x1 j
      library(GenomicAlignments)) Z) ]0 t9 W1 u6 d5 V
      bamfiles <- BamFileList(filenames, yieldSize=2000000)8 I( n( ^1 F: `8 u  W0 r: c, F5 y0 O* f
      exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene"); V) |( h5 G. s* _
      flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)3 T& s7 u8 K  S- r( c/ {
      param <- ScanBamParam(flag=flag)
5 l1 m$ `, W8 ?7 _; }  m( q( |      CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)
: h# K: [6 P1 m! N8 U% m, r! d最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。
! y+ a4 h+ r3 Y  U2 |

farui

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

misser217

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

yinchong42

同期待下文

lchanlon

回复 farui 的帖子
* \% B) b6 B& I  {2 X/ |. z7 g% p
# @5 D; c) @' j. Q# h; N% X4 Z; J, X* |Welcome! Do you have experience on how to evaluate quality of RNA-seq reads?

farui

回复 lchanlon 的帖子% z* a4 u3 o- H5 W
! O2 [$ i4 S4 l6 k
你提醒得很好，RNA-seq的Quality control也是一个重要步骤。我目前正在学习和测试两个工具，一个是RSeQC(python)，另一个是QoRTs(java和R)，等我整理以后在发上来。