## ----style, echo = FALSE, results = 'asis', message = FALSE--------------
BiocStyle::markdown()
library(knitr)

## ----libraryLoad, message = FALSE----------------------------------------
library(metagene)

## ----bamFiles------------------------------------------------------------
bam_files <- get_demo_bam_files()
bam_files

## ----namedBamFiles-------------------------------------------------------
named_bam_files <- bam_files
names(named_bam_files) <- letters[seq_along(bam_files)]
named_bam_files

## ----regionsArgument-----------------------------------------------------
regions <- get_demo_regions()
regions

## ----showDatasets--------------------------------------------------------
data(promoters_hg19)
promoters_hg19

## ----designFile----------------------------------------------------------
fileDesign <- system.file("extdata/design.txt", package="metagene")
design <- read.table(fileDesign, header=TRUE, stringsAsFactors=FALSE)
design$Samples <- paste(system.file("extdata", package="metagene"),
                        design$Samples, sep="/")
kable(design)

## ----alternateDesign-----------------------------------------------------
design <- data.frame(Samples = c("align1_rep1.bam", "align1_rep2.bam",
                     "align2_rep1.bam", "align2_rep2.bam", "ctrl.bam"),
                     align1 = c(1,1,0,0,2), align2 = c(0,0,1,1,2))
design$Samples <- paste0(system.file("extdata", package="metagene"), "/",
                         design$Samples)
kable(design)

## ----minimalAnalysis-----------------------------------------------------
regions <- get_demo_regions()
bam_files <- get_demo_bam_files()
# Initialization
mg <- metagene$new(regions = regions, bam_files = bam_files)
# Plotting
mg$plot(title = "Demo metagene plot")

## ----initialization------------------------------------------------------
regions <- get_demo_regions()
bam_files <- get_demo_bam_files()
mg <- metagene$new(regions = regions, bam_files = bam_files)

## ----showProduceMatrices-------------------------------------------------
mg$produce_matrices()

## ----produceMatricesDesign-----------------------------------------------
mg$produce_matrices(design = design)

## ----produceDataFrame----------------------------------------------------
mg$produce_data_frame(stat = "basic")

## ----showPlot------------------------------------------------------------
mg$plot(region_names = "list1", title = "Demo plot subset")

## ----getParams-----------------------------------------------------------
mg <- get_demo_metagene()
mg$get_params()

## ----getDesign-----------------------------------------------------------
mg$produce_matrices(design = get_demo_design())
## Alternatively, it is also possible to add a design without producing the
## matrices:
#mg$add_design(get_demo_design())
mg$get_design()

## ----getBamCount---------------------------------------------------------
mg$get_bam_count(bam_files[1])

## ----getRegions----------------------------------------------------------
mg$get_regions()

## ----getRegionsSubset----------------------------------------------------
mg$get_regions(region_names = c(regions[1]))

## ----getRawCoverages-----------------------------------------------------
coverages <- mg$get_raw_coverages()
coverages[[1]]
length(coverages)

## ----getRawCoveragesSubset-----------------------------------------------
coverages <- mg$get_raw_coverages(filenames = bam_files[1:2])
length(coverages)

## ----showChain-----------------------------------------------------------
rg <- get_demo_regions()
bam <- get_demo_bam_files()
d <- get_demo_design()
title <- "Show chain"
mg <- metagene$new(rg, bam)$produce_matrices(design = d)$plot(title = title)

## ----copyMetagene--------------------------------------------------------
mg_copy <- mg$clone()

## ----memory, collapse=TRUE-----------------------------------------------
mg1 <- metagene$new(bam_files = bam_files, regions = regions[1])
mg1$produce_data_frame()
mg2 <- metagene$new(bam_files = bam_files, regions = regions[2])
mg2$produce_data_frame()

## ----extractDF-----------------------------------------------------------
df1 <- mg1$get_data_frame()
df2 <- mg2$get_data_frame()
df <- rbind(df1, df2)

## ----plotMetagene--------------------------------------------------------
p <- plot_metagene(df)
p + ggplot2::ggtitle("Managing large datasets")

## ----extractMatrices-----------------------------------------------------
matrices <- mg$get_matrices()
m1 <- matrices[["list1"]][["align1"]][["input"]]
m2 <- matrices[["list1"]][["align2"]][["input"]]

## ----similaRpeak---------------------------------------------------------
library(similaRpeak)
perm_fun <- function(profile1, profile2) {
    similarity(profile1, profile2)[["metrics"]][["RATIO_NORMALIZED_INTERSECT"]]
}

## ----calculateRNI--------------------------------------------------------
ratio_normalized_intersect <- perm_fun(colMeans(m1), colMeans(m2))
ratio_normalized_intersect

## ----permTest------------------------------------------------------------
permutation_results <-  permutation_test(m1, m2, sample_size = nrow(m1),
                                         sample_count = 1000, FUN = perm_fun)

## ----perm_pval-----------------------------------------------------------
sum(ratio_normalized_intersect >= permutation_results) / length(permutation_results)