## ----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)

## ----showProduceTable------------------------------------------------------
mg$produce_table()

## ----produceTableDesign----------------------------------------------------
mg$produce_table(design = design)

## ----produceDataFrame------------------------------------------------------
mg$produce_data_frame()

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

## ----getTable--------------------------------------------------------------
mg <- get_demo_metagene()
mg$produce_table()
mg$get_table()

## ----getMatrices-----------------------------------------------------------
mg <- get_demo_metagene()
mg$produce_table()
m <- mg$get_matrices()
# m$list1$ctrl$input to access to region 'list1' and 'ctrl' design

## ----getDataFrame----------------------------------------------------------
mg <- get_demo_metagene()
mg$produce_table()
mg$produce_data_frame()
mg$get_data_frame()

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

## ----getDesign-------------------------------------------------------------
mg$produce_table(design = get_demo_design())
## Alternatively, it is also possible to add a design without producing the
## table:
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_table(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")

## ----extract_subtables-----------------------------------------------------
tab <- mg$get_table()
tab0 <- tab[which(tab$region == "list1"),]
tab1 <- tab0[which(tab0$design == "align1"),]
tab2 <- tab0[which(tab0$design == "align2"),]

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

## ----calculateRNI----------------------------------------------------------
ratio_normalized_intersect <- 
 perm_fun(tab1[, .(moy=mean(value)), by=bin]$moy, 
        tab2[, .(moy=mean(value)), by=bin]$moy)
ratio_normalized_intersect

## ----permTest--------------------------------------------------------------
permutation_results <- permutation_test(tab1, tab2, sample_size = 50,
                                        sample_count = 1000, FUN = perm_fun)

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