R与bioconductor--Short Read(读取fastq) Rsamtools
2017-08-11 11:39
387 查看
博主自学了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
最后是完整代码片段
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
相关文章推荐
- Calling SNPs/INDELs with SAMtools/BCFtools
- html5 使用FileReader对象的readAsDataURL方法来读取图像文件
- READ TABLE 读取数据(转载自ignativshoo)
- 如何使用read命令读取文件的每一行
- ABAP:DYNP_VALUES_READ读取屏幕字段值
- samtools faidx 命令处理fasta序列
- 安装虚拟机出现“正在进行简易安装时,无法手动启动vmware tools"以及出现解压vmware tools时,Read-only file system
- java 读取配置(read config from file)
- 【R语言】R读取含中文excel文件,read.xlsx乱码问题
- Android-Can't read [D:\tools\android\SDK\platforms\android-23\optional\org.apache.http.legacy.jar
- 有关fstream::read()读取错误问题
- 用DataSet.ReadXml读取无Scheme的XML提速方法
- java实现附件预览(openoffice+swftools+flexpaper)(解决jsp读取全盘文件问题)
- OSG osgDB::readNodeFile 读取ply文件
- [转]SAP ABAP中使用Read_Text函数读取项目文本的方法
- 黑马程序员之C#学习笔记:使用Stream.BeginRead方法读取FileStream的流内容
- [Linux文件]使用read函数从文件读取数据的实例
- insmod: short read错误
- samtools error
- java 为什么InputStream.read()读取一个byte却返回一个int呢?