干细胞之家 - 中国干细胞行业门户第一站

标题: RNA-sequencing数据处理流程 [打印本页]

作者: farui 时间: 2015-8-14 18:14 标题: RNA-sequencing数据处理流程

RNA-sequencing数据处理流程
5 y+ F) Z# D" L% Z, E! O
RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。

因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。

RNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤
第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：' y6 W* z" w5 }" U  ~
   library(SRAbd)
   library(DBI)
   if(!file.exists('SRAmetadb.sqlite')) 7 w# \# K9 y- h2 o! D7 e2 D1 X
      srafile <- getSRAdbFile()
   else& ~! s, @1 U1 e# Y0 y' W
      srafile <- 'SRAmetadb.sqlite'' s5 H) t" ~8 E/ N- D! B* s
   con = dbConnect(RSQLite::SQLite(), srafile)' R+ j  ]0 |, h* B4 s' s# s  H/ [
   df_srr <- listSRAfile(sra_accn, con) #sra_accn为已知SRA Accession号4 A8 `2 [: }# {0 n, t  {4 k
   l_run <- unlist(df_srr["run"])2 X9 D. m; _4 p( n8 n: U
   getSRAfile(l_run, con, fileType = 'sra')
如果不知道SRA Accession号，可以通过NCBI网站搜索。& I  Q1 `' A# F$ W
第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。, @4 g+ s- T: J  y; d  H( ^( v
第三步：对应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片段到参考基因组。具体的操作如下（命令行）：+ s. j/ z* g3 f1 ]: H+ K  t5 C; j
   STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
产生完参考基因组之后，就是对应RNA片段到参考基因组7 U, q( Q# C; }, Y0 D
   STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下6 x0 P5 _! x9 @- c3 f  j- P* e/ q
   samtools view -bS xxxxx.sam -o xxxx.bam
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：
1）我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；
2）产生一个GRangesList对象；
3）把产生的若干.bam文件放入BamFileList对象中；+ e6 x& [8 P, f$ X6 M5 Q# J$ P' d
4）使用summarizeOverlaps方法产生count matrix。
具体操作如下：
   library("TxDb.Hsapiens.UCSC.hg19.knownGene")
   library("Rsamtools")
   library(GenomicAlignments)5 d3 m; V# r! Z' P+ d8 B
   bamfiles <- BamFileList(filenames, yieldSize=2000000)3 Q& d$ _" S7 F5 {
   exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")0 X; n5 T( D+ l" S* ^2 K, z+ p
   flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
   param <- ScanBamParam(flag=flag)
   CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)
最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。% h1 n. c# M! T% w4 w  _

作者: farui 时间: 2015-8-16 00:40

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

作者: misser217 时间: 2015-8-16 08:15

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

作者: yinchong42 时间: 2015-8-17 14:52

同期待下文

作者: lchanlon 时间: 2015-8-17 23:40

回复 farui 的帖子0 ?+ X2 G9 n" f6 \$ y9 ]; T
# R( J4 L- Y, l2 u- o; N6 `, g
Welcome! Do you have experience on how to evaluate quality of RNA-seq reads?

作者: farui 时间: 2015-8-18 04:18

回复 lchanlon 的帖子8 n; ^. [2 c; [7 `: n. C, W
1 ^5 a0 O9 C5 A2 @5 o" B
你提醒得很好，RNA-seq的Quality control也是一个重要步骤。我目前正在学习和测试两个工具，一个是RSeQC(python)，另一个是QoRTs(java和R)，等我整理以后在发上来。

欢迎光临干细胞之家 - 中国干细胞行业门户第一站 (http://www.stemcell8.cn/)