# 使用R按照指定宽度计算覆盖率

[code lang="R"]
library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg19) ## human
tiles <- tileGenome(seqinfo(Hsapiens), tilewidth=200, cut.last.tile.in.chrom=TRUE) ##
tail(tiles[seqnames(tiles)=="chr1"]) ##
library(Rsamtools)
library(GenomicAlignments)
binnedAverage <- function(bins, numvar, mcolname)
{
stopifnot(is(bins, "GRanges"))
stopifnot(is(numvar, "RleList"))
stopifnot(identical(seqlevels(bins), names(numvar)))
bins_per_chrom <- split(ranges(bins), seqnames(bins))
means_list <- lapply(names(numvar),
function(seqname) {
views <- Views(numvar[[seqname]],
bins_per_chrom[[seqname]])
viewMeans(views)
})
new_mcol <- unsplit(means_list, as.factor(seqnames(bins)))
mcols(bins)[[mcolname]] <- new_mcol
bins
}

library(rtracklayer)
bams <- dir(pattern="*.bam\$")
for(bam in bams){
dat <- granges(dat)
score <- coverage(dat)
tiles.s <- keepSeqlevels(tiles, names(score))
bins <- binnedAverage(tiles.s, score, "score")
export(bins, gsub(".bam", ".bedgraph", bam), "BedGraph")
}
[/code]

2.将会跨窗口的按每个窗口分配一次，并使其的宽度设定为1.

4.计算coverage。

5.计算binnedSums。

[code lang="R"]
library(GenomicRanges)
library(BSgenome.Mmusculus.UCSC.mm10) ## mouse
tiles <- tileGenome(seqinfo(Mmusculus), tilewidth=200, cut.last.tile.in.chrom=TRUE) ##
edges <- tiles
width(edges) <- 1
tail(tiles[seqnames(tiles)=="chr1"]) ##
library(Rsamtools)
library(GenomicAlignments)
binnedSums <- function(bins, numvar, mcolname)
{
stopifnot(is(bins, "GRanges"))
stopifnot(is(numvar, "RleList"))
stopifnot(identical(seqlevels(bins), names(numvar)))
bins_per_chrom <- split(ranges(bins), seqnames(bins))
means_list <- lapply(names(numvar),
function(seqname) {
views <- Views(numvar[[seqname]],
bins_per_chrom[[seqname]])
viewSums(views)
})
new_mcol <- unsplit(means_list, as.factor(seqnames(bins)))
mcols(bins)[[mcolname]] <- new_mcol
bins
}

library(rtracklayer)
bams <- dir(pattern="*.bam\$")
for(bam in bams){
dat <- granges(dat)
ol <- findOverlaps(dat, edges)
dat.ol <- dat[queryHits(ol)]
dat.nol <- dat[-queryHits(ol)]
dat.ols <- dat.ol ## note the bin size must larger than reads length
width(dat.ols) <- 1
dat.ole <- dat.ol
start(dat.ole) <- end(dat.ole)
width(dat.nol) <- 1
dat <- c(dat.nol, dat.ols, dat.ole)
dat <- sort(dat)
score <- coverage(dat)
tiles.s <- keepSeqlevels(tiles, names(score))
bins <- binnedSums(tiles.s, score, "score")
bins <- bins[bins\$score>0]
export(bins, gsub(".bam", ".bedgraph", bam), "BedGraph")
}
[/code]