RNA-sequencing数据处理流程

farui

RNA-sequencing数据处理流程2 `* |" F' t+ T8 T
  p; M! Z! x4 v8 h- l
RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。0 \' L% s; q# M; V: J0 D& X
" A+ A2 L: ~" W) ^5 b
因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。# ?9 ~2 c, b7 x4 O! A8 b3 P2 m* S

8 G0 s5 l2 J& d, `! HRNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤+ O7 x, T# x' Z+ w0 W- u8 R. L
第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：
3 P8 N1 Z9 ~: n9 X5 ~      library(SRAbd)/ D3 ~& o& T$ T  e0 G
      library(DBI)
: q( k9 O, ~7 `4 O' E  v0 e      if(!file.exists('SRAmetadb.sqlite'))
0 ^' G/ F* w+ _1 ?3 F& S+ r        srafile <- getSRAdbFile()
, f( \. D2 p5 H/ E# h      else' x1 l, d7 f" v1 w, i0 ?- C% U6 k" |
        srafile <- 'SRAmetadb.sqlite'
" w/ P" T' D% K8 a! n! t0 C      con = dbConnect(RSQLite::SQLite(), srafile)0 X6 n6 N( A* z  H
      df_srr <- listSRAfile(sra_accn, con)                #sra_accn为已知SRA Accession号  h% b- n( ]. Q8 K: W6 {
      l_run <- unlist(df_srr["run"])
% [$ w  I3 t5 {      getSRAfile(l_run, con, fileType = 'sra')  X8 I4 ]8 J2 L9 R6 c
如果不知道SRA Accession号，可以通过NCBI网站搜索。
6 Z, k* P0 j3 T, o' C7 D" i6 ]第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。/ R' `! J  m' H9 h/ k9 ^/ R
第三步：对应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片段到参考基因组。具体的操作如下（命令行）：
! I5 R: `: m( P/ q$ L8 y; x      STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
! M4 Q+ w3 I$ e2 g6 l) I产生完参考基因组之后，就是对应RNA片段到参考基因组
8 U6 {2 h9 V. J( R' |      STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
" k( p$ P6 J* I" R% J/ g这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下
4 C5 J" H, [, ]9 o9 Z% |  O6 c      samtools view -bS xxxxx.sam -o xxxx.bam6 n$ n( `1 q( e* \" y
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：
( }& e- y, f( @  T% b5 \        1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；. Q7 \3 Q0 I8 y
        2） 产生一个GRangesList对象；
5 b7 z1 {' J  n# R, |% M/ a1 m" L  n- T        3） 把产生的若干.bam文件放入BamFileList对象中；; t8 |* s) j, x8 B( R2 B. C* O
        4） 使用summarizeOverlaps方法产生count matrix。) U/ r) i+ i+ L: Y7 K% |: Q
具体操作如下：! n* e  L$ b& k& _/ T$ m
      library("TxDb.Hsapiens.UCSC.hg19.knownGene")! A0 |, c* w! k/ L. x% |
      library("Rsamtools")- R9 L, I! P2 n; u
      library(GenomicAlignments)4 G5 L, |4 s, J6 I$ h
      bamfiles <- BamFileList(filenames, yieldSize=2000000)9 I- x) A; C: M: q. s3 F! B! r
      exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene"). N- q6 j. e6 f9 l
      flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
) s" E0 M: w  e. p/ l      param <- ScanBamParam(flag=flag)
; d2 E7 u  g( ]; ]      CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)
/ V6 e* ~# l& z/ f) V7 p最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。" f6 X' Z' O( I6 @( ^/ _

farui

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

misser217

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

yinchong42

同期待下文

lchanlon

回复 farui 的帖子
7 M) F5 f7 g# N3 R  H4 w1 `  \" f, `2 }
Welcome! Do you have experience on how to evaluate quality of RNA-seq reads?

farui

回复 lchanlon 的帖子9 j1 Y2 c/ N" N( C# Z& _  z: K8 G

3 p; v( C$ a& k2 q5 Y6 M你提醒得很好，RNA-seq的Quality control也是一个重要步骤。我目前正在学习和测试两个工具，一个是RSeQC(python)，另一个是QoRTs(java和R)，等我整理以后在发上来。