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

博主自学了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]#读取序列IDas(quality(reads),"matrix")[1:2,1:10]#转化序列质量library(Rsamtools)bamPath <- system.file("extdata","ex1.bam",package = "Rsamtools")bamFile <- BamFile(bamPath)bamFileseqinfo(bamFile)aln <- scanBam(bamFile)length(aln)class(aln)names(aln)aln <- aln[[1]]names(aln)lapply(aln,function(xx)xx[1])#取出每个列表里面的第一个元素yieldSize(bamFile) <- 1bamFile#此刻yieldSize: 1 每次只读取一行open(bamFile)scanBam(bamFile)[[1]]$seqscanBam(bamFile)[[1]]$seqscanBam(bamFile)[[1]]$seqgr <- GRanges(seqnames = "seq2",ranges = IRanges(start = c(100,1000),end =c(1500,2000)))grparams<- ScanBamParam(which = gr,what = scanBamWhat())scanBamWhat()aln <- scanBam(bamFile,param = params)names(aln)head(aln[[1]]$pos)#有些reads很长,与设置的100边界重叠bamView <- BamViews(bamPath)#读取多个bam文件bamViewaln <- scanBam(bamView)#读入bam文件names(aln[[1]][[1]])bamRanges(bamView) <- gr#对BamViews设置rangesaln<- scanBam(bamView)names(aln[[1]])quickBamFlagSummary(bamFile)#快速读取bam文件,给出summary