- 积分
- 326
- 威望
- 326
- 包包
- 946
|
RNA-sequencing数据处理流程9 L2 v3 p# z2 Q0 |, Y
P _4 l/ V9 D0 q* Y2 d; ARNA-sequencing技术已经逐渐取代microarray技术,成为转录组研究中的重要的高通量技术。我以前作了若干年microarray基因表达数据分析,也就是最近才开始作RNA-sequencing数据分析,现在把我了解的RNA-sequencing数据分析的流程写下来,一方面是给我自己增加印象,而且如果有错误的话希望得到大家的指正,另一方面是希望能对也使用RNA-sequencing技术的朋友有一些启发作用。
1 O% x8 @) b. p7 n6 v: @) O, U9 G
; v6 h6 R1 r) M2 X0 ~" e$ l4 z9 j$ w因为目前RNA-sequencing处理流程的每个环节上都有不止一种方法,这里我着重叙述我使用的方法,对于其他的方法我只能提及,但无法一一详述。另一方面,我这里处理流程从原始sra文件开始,到形成表达数据计数矩阵(count matrix)为止,所以严格说,这篇应该叫RNA-sequencing数据“预”处理流程。之后的differential expression分析,聚类分析,GO term分析,binding site分析以及pathway分析,在此我也不作详述,如果有朋友感兴趣,我可以在未来一一介绍。8 B& j2 I* m Y+ l
: W, l- U2 G* T, F) k2 v
RNA-sequencing处理流程和microarray处理流程在“预”处理阶段差异较大,而在后续的处理差异较小(仅differential expression分析有差别)。RNA-sequencing“预”处理主要分四个步骤
7 t* {2 m' L3 D% J第一步:下载.sra文件。在Bioconductor的SRAdb包里有方法可以下载.sra文件。假设有SRA Accession号,具体的操作如下:. r) R- K; i- v% A- P' F$ I
library(SRAbd), P( Z2 e% V8 b7 r
library(DBI)
' ~/ K2 X! p! h: W/ ^ if(!file.exists('SRAmetadb.sqlite')) ) t9 k1 c( P( f6 R/ U
srafile <- getSRAdbFile()1 J" ` D3 P$ l5 O$ f
else
1 A" E* ^. M& u srafile <- 'SRAmetadb.sqlite'
8 B1 Z0 E3 N; U& _ con = dbConnect(RSQLite::SQLite(), srafile)
" s8 o0 X$ t. J! b- _* [ df_srr <- listSRAfile(sra_accn, con) #sra_accn为已知SRA Accession号3 L. G& k- `9 `
l_run <- unlist(df_srr["run"])
7 f- R# w$ [: u8 l" t0 a getSRAfile(l_run, con, fileType = 'sra')5 _* D, q _3 S/ A
如果不知道SRA Accession号,可以通过NCBI网站搜索。
% u! e: G) T, _, s第二歩:转化.sra文件为.fastq文件。我们需要下载SRA Toolkit,使用其中的fastq-dump命令来完成转换。- N3 C: z9 M4 L9 Q+ I+ ^
第三步:对应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片段到参考基因组。具体的操作如下(命令行):, i: R3 Y+ o3 p6 X
STAR --runThreadN 这里是个数字表示多少线程 --runMode 在产生参考基因组阶段这里必须填genomeGenerate选项 --genomeDir 这里填你希望产生的参考基因组存放的路径 --genomeFastaFiles 这里填下载的FASTA文件路径 --sjdbGTFfile 这里填下载的GTF文件路径
% S: h/ ?4 a$ w) [% V产生完参考基因组之后,就是对应RNA片段到参考基因组! P( M9 L7 c8 p
STAR --runThreadN 这里是个数字表示多少线程 --genomeDir 这里填产生好的参考基因组存放的路径 --readFilesIn 这里填已经准备好的.fasta文件路径: d; _* ^% Q1 Z
这样就可以产生.sam文件,我们再需要通过samtools把.sam转化成.bam文件,操作如下
; U- J/ G+ V; n3 Z+ | samtools view -bS xxxxx.sam -o xxxx.bam
1 O5 d" s. [# E' a第四步,我们要利用产生的.bam文件获得count matrix。这也有几个工具可以实现,比如python的HTSeq包,Bioconductor的Rsubread包和GenomicAlignments包。我使用GenomicAlignments包。完成这个任务,也需要几个步骤:$ {$ w0 ~& c" A7 n h* X
1) 我们需要定义基因模型(gene models)。简单的方法是直接使用Bioconductor现有的,比如TxDb.Hsapiens.UCSC.hg19.knowngene包;
2 Z* i/ Z- e. Y 2) 产生一个GRangesList对象;9 X- `% o( E6 O6 ]% x1 T( g
3) 把产生的若干.bam文件放入BamFileList对象中;
4 e. k% {+ }$ P/ }8 q1 n3 q 4) 使用summarizeOverlaps方法产生count matrix。3 v9 Z J; i, P6 d
具体操作如下:/ A" T# e5 H8 l' y
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
; {( z+ l; p, C9 H4 h$ r library("Rsamtools")9 P' A( I M2 l
library(GenomicAlignments)
4 f! A5 t$ @1 {" C$ v bamfiles <- BamFileList(filenames, yieldSize=2000000)
" {9 [& p( r) k. I" @1 ` exByGn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")
& }; c! W3 m+ y$ W flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=TRUE)
* o" i- {& {) m# B. f* b param <- ScanBamParam(flag=flag)
3 \8 b2 p y3 e CntMat <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=FALSE, param=param)
$ E. z# M, r/ \, i* f最终,我们就获得了CntMat,可以开始一系列后续的分析处理了。6 R" j9 e; U C( s/ S
|
-
总评分: 威望 + 20
包包 + 50
查看全部评分
|