avoid_gaps_update: Is is a function designed to remove values <= to...

In CharlesJB/metagene: A package to produce metagene plots

Description

Is is a function designed to remove values <= to 'gaps_threshold'. Nucleotides local and global positions, bins, size of regions/genes and exons will be recalculated. To use on metagene's table during RNA-seq analysis. Not made for ChIP-Seq analysis or to apply on matagene's data_frame. A similar function is implemented in produce_data_frame() with same arguments. The unique goal of this function is to allow permutation_test which match the plot created using avoid_gaps, bam_name and gaps_threshold arguments in the produce_data_frame function.

Usage

avoid_gaps_update(table, bam_name, gaps_threshold = 0)

Arguments

table	A data.table from produce_table(...) function of metagene.
bam_name	A reference bam_name to allow the same removal (position in bam) of values for other bam file.
gaps_threshold	A threshold under which values will be removed.

Value

A data.table with values <= to 'gaps_threshold' removed

Examples

## Not run: 
bam_files <- c(
system.file("extdata/c_al4_945MLM_demo_sorted.bam", package="metagene"),
system.file("extdata/c_al3_362PYX_demo_sorted.bam", package="metagene"),
system.file("extdata/n_al4_310HII_demo_sorted.bam", package="metagene"),
system.file("extdata/n_al3_588WMR_demo_sorted.bam", package="metagene"))
region <- c(
system.file("extdata/ENCFF355RXX_DPM1less.bed", package="metagene"),
system.file("extdata/ENCFF355RXX_NDUFAB1less.bed", package="metagene"),
system.file("extdata/ENCFF355RXX_SLC25A5less.bed", package="metagene"))
mydesign <- matrix(c(1,1,0,0,0,0,1,1),ncol=2, byrow=FALSE)
mydesign <- cbind(c("c_al4_945MLM_demo_sorted.bam",
                   "c_al3_362PYX_demo_sorted.bam",
                   "n_al4_310HII_demo_sorted.bam",
                   "n_al3_588WMR_demo_sorted.bam"), mydesign)
colnames(mydesign) <- c('Samples', 'cyto', 'nucleo')
mydesign <- data.frame(mydesign)
mydesign[,2] <- as.numeric(mydesign[,2])-1
mydesign[,3] <- as.numeric(mydesign[,3])-1

mg <- metagene$new(regions = region, bam_files = bam_files, 
                                                           assay = 'rnaseq')
mg$produce_table(flip_regions = FALSE, bin_count = 100, 
                               design = mydesign, normalization = 'RPM')
mg$produce_data_frame(avoid_gaps = TRUE, 
                       bam_name = "c_al4_945MLM_demo_sorted", 
                       gaps_threshold = 10)
mg$plot()
tab <- mg$get_table()
tab <- avoid_gaps_update(tab, 
       bam_name = 'c_al4_945MLM_demo_sorted', gaps_threshold = 10)
tab1 <- tab[which(tab0$design == "cyto"),]
tab2 <- tab[which(tab0$design == "nucleo"),]

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

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

permutation_results <- permutation_test(tab1, tab2, sample_size = 2,
                               sample_count = 1000, FUN = perm_fun)
hist(permutation_results, 
       main="ratio_normalized_intersect (1=total overlapping area)")
abline(v=ratio_normalized_intersect, col = 'red')
sum(ratio_normalized_intersect >= permutation_results) / 
       length(permutation_results)

## End(Not run)

CharlesJB/metagene documentation built on July 11, 2021, 11:48 a.m.