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

最新推荐文章于 2024-03-06 23:00:16 发布

wangyunpeng_bio

最新推荐文章于 2024-03-06 23:00:16 发布

阅读量8.6k

点赞数 10

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

本文链接：https://blog.csdn.net/qq_29300341/article/details/77094998

版权

coursera 同时被 2 个专栏收录

27 篇文章 3 订阅

订阅专栏

6 篇文章 0 订阅

订阅专栏

博主自学了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