在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文件到本地。

image

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

多组学分析一样

也有教程：https://www.bioconductor.org/help/course-materials/2015/EMBO/A01_RBiocForSequenceAnalysis.html

值得花10个小时学习的。