|
 
- 积分
- 326
- 威望
- 326
- 包包
- 946
|
RNA-sequencing数据处理流程/ q3 U5 k J3 f U/ k$ W' b
# \" `! n5 [ L3 n/ ~4 T6 w, A4 ?
RNA-sequencing技术已经逐渐取代microarray技术，成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析，也就是最近才开始作RNA-sequencing数据分析，现在把我了解的RNA-sequencing数据分析的流程写下来，一方面是给我自己增加印象，而且如果有错误的话希望得到大家的指正，另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。
" w3 W! l4 O A' P5 ]& p
1 k5 Z+ c, b7 @因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法，这里我着重叙述我使用的方法，对于其他的方法我只能提及，但无法一一详述。另一方面，我这里处理流程从原始sra文件开始，到形成表达数据计数矩阵（count matrix）为止，所以严格说，这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析，聚类分析，GO term分析，binding site分析以及pathway分析，在此我也不作详述，如果有朋友感兴趣，我可以在未来一一介绍。4 I: d( i# B1 k
9 I6 k0 h8 v: T; l0 p
RNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大，而在后续的处理差异较小（仅differential expression分析有差别）。RNA-sequencing“预”处理主要分四个步骤
8 {- h9 v- a8 w8 f* A' I- w第一步：下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号，具体的操作如下：" @, ^# h( ?7 ~# C; \
library(SRAbd)
, y! b% c' x) I% X( ~: X library(DBI)
+ ^7 b/ y7 J3 C if(!file.exists('SRAmetadb.sqlite')) 9 S( Z+ h6 X/ X; J# W: B
srafile <- getSRAdbFile()
3 _# C9 M& N4 L. }& h3 C1 j* ` else
6 f9 P! S; u9 {% e8 u- `7 n1 ~ srafile <- 'SRAmetadb.sqlite'
3 l# L9 @ X6 X con = dbConnect(RSQLite::SQLite(), srafile)
4 O- R& I3 Q+ R% [7 w; L. E df_srr <- listSRAfile(sra_accn, con) #sra_accn为已知SRA Accession号
# n* \ S) |) U! D* k+ T l_run <- unlist(df_srr["run"])+ e: o. L4 z/ r8 a6 h( m3 \8 _
getSRAfile(l_run, con, fileType = 'sra')
5 t9 ~% n4 w Y0 j4 m4 Z# g如果不知道SRA Accession号，可以通过NCBI网站搜索。5 t* m: X$ i( D+ P
第二歩：转化.sra文件为.fastq文件。我们需要下载SRA Toolkit，使用其中的fastq-dump命令来完成转换。
+ ]! A6 `& _' i0 ]1 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片段到参考基因组。具体的操作如下（命令行）：. B9 q& h; h" p- {8 n
STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径; a/ W4 s7 P0 k% R& D$ i% _2 @
产生完参考基因组之后，就是对应RNA片段到参考基因组: m3 d- L ^! M
STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径
$ m) S f0 _; }: G8 K" I这样就可以产生.sam文件，我们再需要通过samtools把.sam转化成.bam文件，操作如下
9 ~7 I3 @$ |( ~; Z1 \8 |2 y samtools view -bS xxxxx.sam -o xxxx.bam0 T! }# K" Y, f v0 a
第四步，我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现，比如python的HTSeq包，Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务，也需要几个步骤：) y4 _: ], n9 O F D0 d( ^& _7 F* y
1） 我们需要定义基因模型（gene models）。简单的方法是直接使用Bioconductor现有的，比如TxDb.Hsapiens.UCSC.hg19.knowngene包；
. H4 X1 [1 y9 A! R1 m) k 2） 产生一个GRangesList对象；
& `( z/ A" d+ t' { 3） 把产生的若干.bam文件放入BamFileList对象中；' w3 O* x ^3 h/ h/ ^% f/ A
4） 使用summarizeOverlaps方法产生count matrix。9 l' u+ T I2 [4 g. B$ y" @8 {
具体操作如下：
# s( ?" M v, Z1 `& B' | library("TxDb.Hsapiens.UCSC.hg19.knownGene")- p% i# F0 x$ e% S" v8 k' Q
library("Rsamtools")
" F$ `$ K1 ~# p8 J# G library(GenomicAlignments)
7 _$ W1 v2 \8 |% H% U" _3 v bamfiles <- BamFileList(filenames, yieldSize=2000000)
" r2 ]# T8 x7 f6 { ^ R exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")% U* z2 ?$ _, Q" g) O
flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
" Q" u# y8 U3 o ` param <- ScanBamParam(flag=flag)5 o# S( [0 K p' u, N1 b& S# E
CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param); y4 a- V0 J* `4 W; ]: b
最终，我们就获得了CntMat，可以开始一系列后续的分析处理了。* \$ ^# x% b# F! K2 h
|
-
总评分: 威望 + 20
包包 + 50
查看全部评分
|
发表于 2015-8-14 18:14
发表于 2015-8-17 23:40
