文件类

类型	示例使用	的名字	包
请全部	范围注释	Bedfile（）	rtracklayer
.wig	报道	WigFile (),BigWigFile ()	rtracklayer
.gtf	记录模型	GTFFile ()	rtracklayer
		makeTxDbFromGFF ()	GenomicFeatures
.2bit	基因组序列	TwoBitFile ()	rtracklayer
.fastq	读取和品质	FastqFile ()	ShortRead
本演讲	对齐的读取	BamFile ()	Rsamtools
.tbx	索引一样	tabixfile（）	Rsamtools
.vcf	变体呼叫	VcfFile ()	VariantAnnotation.

## rtracklayer menagerie库(rtracklayer)名称(getClass("RTLFile")@subclasses)

##[1]“UCSCFile”“GFFFile”“BEDFile”“##[4]”“WIGFile”“BigWigFile”“ChainFile”“##[7]”“TwoBitFile”“FastaFile”“TabSeparatedFile”“##[10]”“CompressedFile”“GFF1File”“GFF2File”“##[13]”“GFF3File”“BEDGraphFile”“BED15File”“##[16]”“GVFFile”“GZFile”“BGZFile”“##[22]”BZ2File”“XZFile”

笔记

• 不一致的接口
• open (),close (),进口（）/产量()/读* ()
• 一些:通过索引选择性导入(例如，.baibam指数);选择(“列”);限制(“行”)

管理文件集合

* filelist（）类

• reduceByYield ()-遍历单个大文件
• bplapply ()(BiocParallel） - 并行对多个文件执行独立操作

GenomicFiles ()

• “行”作为基因组范围限制，'列'作为文件
• 每个行x列是一个地图从文件数据到有用的表示R
• reduceByRange (),reduceByFile ()：折叠地图入摘要表示
• 看到genomicfiles小插图图1

高效的代码

将following定义为一个函数。

f0 < - 函数（n）{##效率低下！ANS < -  numeric（）for（i在seq_len（n））ANS < -  C（ANS，EXP（I））ANS}

使用system.time ()为了探究执行as需要多长时间n从100增加到10000。使用相同的()和微基准测试比较替代方案f1 (),f2 (),f3 ()这三个不同函数的正确性和性能。这些函数使用什么策略?

f <- function(n) {ans <- numeric(n) for (i in seq_len(n)) ans[[i] <- exp(i) ans} f2 <- function(n) sapply(seq_len(n)， exp) f3 <- function(n) exp(seq_len(n))

连续睡眠和平行睡眠

去睡觉1秒钟，然后返回我。这需要8秒钟。

## parallel f1 <- function(i) {Sys.sleep(1) i} ## serial f0 <- function(n) lapply(seq_len(n)， fun) ## parallel f1 <- function(n) bplapply(seq_len(n)， fun)

读取重叠的窗口

此练习使用以下软件包：

图书馆（基因组）库（基因组）库（BiocomicFallels）库（Rsamtools）库（Genomeinfodb）

这是一个基本的重新实现csaw被计数算法。假设芯片片段长度为110nt，我们将覆盖范围置于宽度为50的窗口中。我们关注chr1。

碎片弹。len <- 110间距<- 50 chr <- "chr1"

这里我们指向bam文件，表示我们将以大小为1,000,000的块处理这些文件。

FLS < -  dir（“〜/ usebiocumon-data / chipseq / nfya /”，“.bam $”，full = true）名称（fls）< -  sub（“fastq。*”，“，baseName（vs））BFL < -  Bamfilerist（FLS，ExitingSize = 1000000）

我们将使用tileGenome ()，聚焦于“标准”染色体。

[chr] tiles <- tileGenome(len, tilewidth=spacing, cut.last.tile.in.chrom=TRUE)

我们将使用reduceByYield ()来遍历单个文件。我们读来告诉这个函数收益率一大块文件，我们怎么地图从它的输入表示到每个窗口的数据块计数，最后我们将如何减少连续的块变成最终的表示。

收益率应该是一个函数，它接受一个参数，即输入源，并返回一组记录

yield <- function(x，…)

地图必须采取产量的输出，也许额外的参数，并返回数量的矢量。我们将调整描述对齐的基因组范围大小，使得它们具有等于片段长度的宽度

map <-函数(x, tiles, frag。{gr <- keepstandardchromosome (granges(x)) countOverlaps(tiles, resize(gr, fragl .len))}

减少从地图（在我们的情况下，计数的载体）并将它们结合成一个结果。我们只是添加了我们的向量（+实际上是一个功能！）

减少< -“+”

为了处理一个文件，我们使用reduceByYield ()，传递我们想要处理的文件、yield、map和reduce函数。我们的“包装器”函数将任何附加参数传递给reduceByYield ()使用…:

count1file <- function(bf，…)reduceByYield(bf, yield, map, reduce，…)

使用yieldSize和reduceByYield ()意味着我们在处理每个文件时不会消耗太多的内存，这样我们就可以并行地处理文件BiocParallel。这simplify2array ()函数将一个向量列表转换为一个矩阵。

count <- bplapply(bfl, count1file, tiles=tiles, fragl .len= fragl .len) count <- simplify2array(counts) dim(counts) colsum (counts)