|
 
- 积分
- 326
- 威望
- 326
- 包包
- 946
|
RNA-sequencing数据处理流程4 e9 F+ M4 h. U7 K( G8 {/ N! L
# M* _# u* q) S. m3 H
RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。
( ]. P1 a O) S3 @+ D: R' N2 A* ^& |
因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。
, M3 ~0 m$ _; I4 C: b N1 |$ j3 o. F" e% p! k0 X
RNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤; f8 E5 }$ o3 M9 d9 ~
第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：
+ M; q/ s7 A t! q; h library(SRAbd)' q' u* C; ~; b0 s: H3 T+ o
library(DBI)5 D& v. J$ C8 W5 `
if(!file.exists('SRAmetadb.sqlite')) ! E( ~3 L: G I1 \: \9 b! k
srafile <- getSRAdbFile()
1 ~/ |* Z- i& K! u# ^+ z else! u) u* z& l/ W9 R8 ~7 F) F$ y
srafile <- 'SRAmetadb.sqlite'
. N6 X4 W1 S* s" w8 V. k4 `. [ con = dbConnect(RSQLite::SQLite(), srafile)
3 R" J6 A1 U, h df_srr <- listSRAfile(sra_accn, con) #sra_accn为已知SRA Accession号; y8 r/ N- ^& c( c- X" `
l_run <- unlist(df_srr["run"])
0 \9 o0 K3 A" @- X" v* Q getSRAfile(l_run, con, fileType = 'sra')3 h4 V+ |4 t9 y( O) }) T6 t
如果不知道SRA Accession号，可以通过NCBI网站搜索。
) e/ q# w! C$ O! p2 H _5 g7 E第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。
! \, l, b7 c$ H% F8 i; h第三步：对应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片段到参考基因组。具体的操作如下（命令行）：
$ \" a2 K/ z$ q( G" q$ } STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
( C P2 s# J, l4 W产生完参考基因组之后，就是对应RNA片段到参考基因组; t& W; K4 l5 F( Y$ x9 ^
STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
! k j8 i% J1 l' V这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下
6 d6 b, \! `3 s samtools view -bS xxxxx.sam -o xxxx.bam, B; p8 m. U! \" }1 @8 _
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：" ?7 g. q4 v8 c1 w
1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；
4 q. D* }# [. }! K% i 2） 产生一个GRangesList对象；: p8 a; L4 w4 B( O" H
3） 把产生的若干.bam文件放入BamFileList对象中；
* Q3 R2 T( h& c: Z1 w 4） 使用summarizeOverlaps方法产生count matrix。7 d* _: q: M0 r2 W
具体操作如下：
5 Q3 S! X( C" n$ h library("TxDb.Hsapiens.UCSC.hg19.knownGene")
0 E/ |! K" v* I, w0 g8 _ library("Rsamtools")- O6 i1 o, k6 |* s0 q! p, r
library(GenomicAlignments)
N$ v" ]5 a; E bamfiles <- BamFileList(filenames, yieldSize=2000000)/ P/ l7 r+ m5 ^" H
exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
& X1 C8 X. B3 u4 _ N( o/ u$ F' ` flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
, y0 R5 B1 u4 E' o9 K3 z param <- ScanBamParam(flag=flag)+ H0 P( \4 {/ z' R0 w
CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)
7 I: u% k3 a% [% Q6 |0 |最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。
+ v% h" a) i3 U* C# y |
-
总评分: 威望 + 20
包包 + 50
查看全部评分
|
发表于 2015-8-14 18:14
发表于 2015-8-17 23:40
