if(interactive()) { # see example in seekRIP
# Retrieve system files
extdata.dir <- system.file("extdata", package="RIPSeeker")
bamFiles <- list.files(extdata.dir, ".bam$", recursive=TRUE, full.names=TRUE)
bamFiles <- grep("PRC2", bamFiles, value=TRUE)
# Parameters setting
binSize <- 1e5 # use a large fixed bin size for demo only
multicore <- FALSE # use multicore
strandType <- "-" # set strand type to minus strand
################ run main function for HMM inference on all chromosomes ################
mainSeekOutputRIP <- mainSeek(bamFiles=
grep(pattern="SRR039214", bamFiles, value=TRUE, invert=TRUE),
binSize=binSize, strandType=strandType,
reverseComplement=TRUE, genomeBuild="mm9",
uniqueHit = TRUE, assignMultihits = TRUE,
rerunWithDisambiguatedMultihits = TRUE,
multicore=multicore, silentMain=FALSE, verbose=TRUE)
nbhGRRIP <- mainSeekOutputRIP$nbhGRList$chrX
logOddScore <- computeLogOdd(nbhGRRIP)
values(nbhGRRIP) <- cbind(as.data.frame(values(nbhGRRIP)), logOddScore)
enrichIdx <- which(values(nbhGRRIP)$viterbi_state == 2)
unmergedRIP <- nbhGRRIP[enrichIdx]
mergedRIP <- reduce(unmergedRIP, min.gapwidth = median(width(unmergedRIP) ))
overlapIdx <- findOverlaps(mergedRIP, unmergedRIP)
# a list with query hits as names and subject hits as items
findOverlapsHits <- split(overlapIdx, queryHits(overlapIdx))
# get the score for the first merged region
x <- scoreMergedBins(findOverlapsHits[[1]], unmergedRIP, mergedRIP)
# get scores for all of the merged regions
mergedRIPList <- lapply(split(overlapIdx, queryHits(overlapIdx)),
scoreMergedBins, unmergedRIP, mergedRIP)
names(mergedRIPList) <- NULL
mergedRIP <- do.call(c, mergedRIPList)
# logOddScore is the averaged logOddScore across merged bins
mergedRIP
}
Run the code above in your browser using DataLab