R与bioconductor--Short Read(读取fastq) Rsamtools

最新推荐文章于 2024-10-12 22:00:13 发布

wangyunpeng_bio

最新推荐文章于 2024-10-12 22:00:13 发布

阅读量9.3k

点赞数 11

CC 4.0 BY-SA版权

分类专栏： R coursera 文章标签： coursera r语言 bioconductor ShortRead Rsamtools

本文链接：https://blog.youkuaiyun.com/qq_29300341/article/details/77094998

coursera 同时被 2 个专栏收录

27 篇文章

订阅专栏

6 篇文章

订阅专栏

本教程通过Coursera上的课程分享了如何使用Bioconductor进行基因组科学数据分析，包括使用ShortRead包读取fastq文件及利用Rsamtools包操作bam文件的具体步骤。

博主自学了coursera上来自约翰霍普金斯大学<使用Bioconductor分析基因组科学数据>，很不错，推荐给大家

Short Read(读取fastq)

  library(ShortRead)
 
  fl <- system.file(package="ShortRead", "extdata", "E-MTAB-1147",
 
                    "ERR127302_1_subset.fastq.gz")
 
  reads <- readFastq(fl)
 
  fqFile<- FastqFile(fl)
 
  reads <- readFastq(fl)
 
  sread(reads)[1:2]#读取序列
 
  quality(reads)[1:2]#读取序列质量
 
  id(reads)[1:2]#读取序列ID
 
  as(quality(reads),"matrix")[1:2,1:10]#转化序列质量

Rsamtools

  library(Rsamtools)
 
  bamPath <- system.file("extdata","ex1.bam",package = "Rsamtools")
 
  bamFile <- BamFile(bamPath)
 
  bamFile
 
  seqinfo(bamFile)
 
  aln <- scanBam(bamFile)
 
  length(aln)
 
  class(aln)
 
  names(aln)
 
  aln <- aln[[1]]
 
  names(aln)
 
  lapply(aln,function(xx)xx[1])#取出每个列表里面的第一个元素
 
  yieldSize(bamFile) <- 1
 
  bamFile#此刻yieldSize: 1 每次只读取一行
 
  open(bamFile)
 
  scanBam(bamFile)[[1]]$seq
 
  scanBam(bamFile)[[1]]$seq
 
  scanBam(bamFile)[[1]]$seq
 
  gr <- GRanges(seqnames = "seq2",
 
  ranges = IRanges(start = c(100,1000),end =c(1500,2000)))
 
  gr
 
  params<- ScanBamParam(which = gr,what = scanBamWhat())
 
  scanBamWhat()
 
  aln <- scanBam(bamFile,param = params)
 
  names(aln)
 
  head(aln[[1]]$pos)#有些reads很长,与设置的100边界重叠
 
  bamView <- BamViews(bamPath)#读取多个bam文件
 
  bamView
 
  aln <- scanBam(bamView)#读入bam文件
 
  names(aln[[1]][[1]])
 
  bamRanges(bamView) <- gr#对BamViews设置ranges
 
  aln<- scanBam(bamView)
 
  names(aln[[1]])
 
  quickBamFlagSummary(bamFile)#快速读取bam文件,给出summary
 
  最后是完整代码片段

library(ShortRead)
fl <- system.file(package="ShortRead", "extdata", "E-MTAB-1147", 
                  "ERR127302_1_subset.fastq.gz")
reads <- readFastq(fl)
fqFile<- FastqFile(fl)
reads <- readFastq(fl)
sread(reads)[1:2]#读取序列
quality(reads)[1:2]#读取序列质量
id(reads)[1:2]#读取序列ID
as(quality(reads),"matrix")[1:2,1:10]#转化序列质量

library(Rsamtools)
bamPath <- system.file("extdata","ex1.bam",package = "Rsamtools")
bamFile <- BamFile(bamPath)
bamFile
seqinfo(bamFile)
aln <- scanBam(bamFile)
length(aln)
class(aln)
names(aln)
aln <- aln[[1]]
names(aln)
lapply(aln,function(xx)xx[1])#取出每个列表里面的第一个元素
yieldSize(bamFile) <- 1
bamFile#此刻yieldSize: 1 每次只读取一行
open(bamFile)
scanBam(bamFile)[[1]]$seq
scanBam(bamFile)[[1]]$seq
scanBam(bamFile)[[1]]$seq
gr <- GRanges(seqnames = "seq2",
ranges = IRanges(start = c(100,1000),end =c(1500,2000)))
gr
params<- ScanBamParam(which = gr,what = scanBamWhat())
scanBamWhat()
aln <- scanBam(bamFile,param = params)
names(aln)
head(aln[[1]]$pos)#有些reads很长,与设置的100边界重叠
bamView <- BamViews(bamPath)#读取多个bam文件
bamView
aln <- scanBam(bamView)#读入bam文件
names(aln[[1]][[1]])
bamRanges(bamView) <- gr#对BamViews设置ranges
aln<- scanBam(bamView)
names(aln[[1]])
quickBamFlagSummary(bamFile)#快速读取bam文件,给出summary