## ----set-options, echo=FALSE, cache=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- options(width = 400) ## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # source("https://bioconductor.org/biocLite.R") # biocLite("HiCcompare") # library(HiCcompare) ## ---- eval=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # # read in files # mat <- read.table("hic_1000000.matrix") # bed <- read.table("hic_1000000_abs.bed") # # convert to BEDPE # dat <- hicpro2bedpe(mat, bed) # # NOTE: hicpro2bedpe returns a list of lists. The first list, dat$cis, contains the intrachromosomal contact matrices # # NOTE: dat$trans contains the interchromosomal contact matrix which is not used in HiCcompare. ## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # cnv <- get_CNV(path2bam = 'path/to/bamfiles', out.file = 'path/to/bamfiles/outfile', # bin.size = 1000, genome = 'hg19', CNV.level = 2) ## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # data('hg19_blacklist') # # # combine cnv excluded regions with blacklist regions # exclude <- cbind(cnv, hg19_blacklist) ## ---- warning=FALSE, message=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- library(`HiCcompare`) # load the data data("HMEC.chr22") data("NHEK.chr22") head(HMEC.chr22) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # create the `hic.table` object chr22.table <- create.hic.table(HMEC.chr22, NHEK.chr22, chr = 'chr22') head(chr22.table) ## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # chr22.table <- create.hic.table(HMEC.chr22, NHEK.chr22, chr = 'chr22', exclude.regions = exclude, exclude.overlap = 0.2) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # create list of `hic.table` objects data("HMEC.chr10") data("NHEK.chr10") # create the `hic.table` object chr10.table <- create.hic.table(HMEC.chr10, NHEK.chr10, chr = 'chr10') hic.list <- list(chr10.table, chr22.table) head(hic.list) ## ---- echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- HMEC.chr22_BEDPE <- chr22.table[, 1:7, with=FALSE] NHEK.chr22_BEDPE <- chr22.table[, c(1:6, 8), with=FALSE] head(HMEC.chr22_BEDPE) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- bed.hic.tab <- create.hic.table(HMEC.chr22_BEDPE, NHEK.chr22_BEDPE) head(bed.hic.tab) ## ---- eval=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # # first dataset # mat1 <- read.table("hic1_1000000.matrix") # bed1 <- read.table("hic1_1000000_abs.bed") # dat1 <- hicpro2bedpe(mat, bed) # dat1 <- dat1$cis # extract intrachromosomal matrices # # second dataset # mat2 <- read.table("hic2_1000000.matrix") # bed2 <- read.table("hic2_1000000_abs.bed") # dat2 <- hicpro2bedpe(mat, bed) # dat2 <- dat2$cis # extract intrachromosomal matrices # # # for chr1 # hic.table <- create.hic.table(dat1[[1]], dat2[[1]]) # # # for all chromosomes # hic.list <- mapply(create.hic.table, dat1, dat2, SIMPLIFY = FALSE) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- data("hmec.IS") data("nhek.IS") head(hmec.IS) IS.hic.tab <- create.hic.table(hmec.IS, nhek.IS) ## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # hic.list <- total_sum(hic.list) ## ---- fig.width=7, fig.height=4------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Jointly normalize data for a single chromosome hic.table <- hic_loess(chr22.table, Plot = TRUE, Plot.smooth = FALSE) knitr::kable(head(hic.table)) ## ---- message=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Multiple hic.tables can be processed in parallel by entering a list of hic.tables hic.list <- hic_loess(hic.list, parallel = TRUE) ## ---- message=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- filter_params(hic.table) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- hic.table <- hic_compare(hic.table, A.min = 15, adjust.dist = TRUE, p.method = 'fdr', Plot = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- knitr::kable(head(hic.table)) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- IntSet <- make_InteractionSet(hic.table) ## ----fig.width=7, fig.height=4-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- number_of_unitdistances <- 100 # The dimensions of the square matrix to be simualted number_of_changes <- 250 # How many cells in the matrix will have changes i.range <- sample(1:number_of_unitdistances, number_of_changes, replace = TRUE) # Indexes of cells to have controlled changes j.range <- sample(1:number_of_unitdistances, number_of_changes, replace = TRUE) # Indexes of cells to have controlled changes sim_results <- hic_simulate(nrow = number_of_unitdistances, medianIF = 50000, sdIF = 14000, powerlaw.alpha = 1.8, fold.change = 4, i.range = i.range, j.range = j.range, Plot = TRUE, alpha = 0.1) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- names(sim_results) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- sims <- sim_matrix(nrow = number_of_unitdistances, medianIF = 50000, sdIF = 14000, powerlaw.alpha = 1.8, fold.change = 4, i.range = i.range, j.range = j.range) ## ---- fig.height=4, fig.width=7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- MD.plot1(M = hic.table$M, D = hic.table$D, mc = hic.table$mc, smooth = TRUE) ## ---- fig.height=4, fig.width=7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # no p-value coloring MD.plot2(M = hic.table$adj.M, D = hic.table$D, smooth = FALSE) # p-value coloring MD.plot2(M = hic.table$adj.M, D = hic.table$D, hic.table$p.value, smooth = FALSE) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- full.NHEK <- sparse2full(NHEK.chr22) full.NHEK[1:5, 1:5] sparse.NHEK <- full2sparse(full.NHEK) head(sparse.NHEK) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- KR.NHEK <- KRnorm(full.NHEK) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- SCN.NHEK <- SCN(full.NHEK) ## ---- message=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- result <- MA_norm(hic.table, Plot = TRUE)