在R里面使用Rsubread完成组学分析全套流程
你是否害怕linux的黑白命令行操作，是否对去可视化畏畏缩缩，那么你会爱上它：Rsubread
这里演示一下传统的RNA-seq数据的表达量分析全流程, 安装Rsubread包后会有自带的测序数据如下：
-rw-r--r-- 1 jmzeng admin 25K Nov 9 18:04 longreads.txt.gz
-rw-r--r-- 1 jmzeng admin 80K Nov 9 18:04 reads.txt.gz
-rw-r--r-- 1 jmzeng admin 80K Nov 9 18:04 reads1.txt.gz
-rw-r--r-- 1 jmzeng admin 80K Nov 9 18:04 reads2.txt.gz
-rw-r--r-- 1 jmzeng admin 89K Nov 9 18:04 reference.fa下面的分析流程也以此为例子，不过要切记，一旦切换到人类真实数据，下面的步骤都会耗时很可观，要有心理准备哈！
step1：构建索引
需要有参考基因组文件，这里使用Rsubread自带的数据作为例子：
library(Rsubread)
ref <- system.file("extdata","reference.fa",package="Rsubread")
buildindex(basename="reference_index",reference=ref)step2：比对
需要有fastq格式的测序数据，还是使用Rsubread自带的数据作为例子：
## 首先是单端数据
reads <- system.file("extdata","reads.txt.gz",package="Rsubread")
align(index="reference_index",readfile1=reads,output_file="alignResults.BAM",phredOffset=64)
## 下面是双端
reads1 <- system.file("extdata","reads1.txt.gz",package="Rsubread")
reads2 <- system.file("extdata","reads2.txt.gz",package="Rsubread")
align(index="reference_index",readfile1=reads1,readfile2=reads2,
output_file="alignResultsPE.BAM",phredOffset=64)测试数据比对很迅速，也会同步输出bam文件到本地。
step3：定量
需要有基因组特征描述文件，通常是gtf格式的基因，转录本，外显子的染色体，起始终止坐标，这里还是使用测试数据，自己制作 特征描述文件如下：
ann <- data.frame(
GeneID=c("gene1","gene1","gene2","gene2"),
Chr="chr_dummy",
Start=c(100,1000,3000,5000),
End=c(500,1800,4000,5500),
Strand=c("+","+","-","-"),
stringsAsFactors=FALSE)
ann
fc_SE <- featureCounts("alignResults.BAM",annot.ext=ann)
fc_SE
fc_PE <- featureCounts("alignResultsPE.BAM",annot.ext=ann,isPairedEnd=TRUE)
fc_PE是不是很激动，这么简单就完成了NGS组学数据分析一条龙分析啊！！！
还有一些小细节
x <- qualityScores(filename=reads,offset=64,nreads=1000)
x[1:10,1:10]
propmapped("alignResults.BAM")值得注意的是，你只是看了看这个包的用法而已，要想用得好，请听下回分解哦！
推荐看原版50页的PDF说明书哈：https://bioconductor.org/packages/release/bioc/vignettes/Rsubread/inst/doc/SubreadUsersGuide.pdf
其它例子：http://combine-australia.github.io/RNAseq-R/07-rnaseq-day2.html
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2018-12-06，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录