## ----setup, echo=FALSE-------------------------------------------------------- knitr::opts_chunk$set( message = FALSE, warning = FALSE, fig.height = 8, fig.width = 10 ) ## ----install, eval = FALSE---------------------------------------------------- # if (!"BiocManager" %in% rownames(installed.packages())) # install.packages("BiocManager") # pkg <- c("tidyverse", "Rsamtools", "csaw", "BiocParallel", "rtracklayer") # BiocManager::install(pkg) # BiocManager::install("extraChIPs") ## ----wincounts---------------------------------------------------------------- library(tidyverse) library(Rsamtools) library(csaw) library(BiocParallel) library(rtracklayer) bfl <- system.file( "extdata", "bam", c("ex1.bam", "ex2.bam", "input.bam"), package = "extraChIPs" ) %>% BamFileList() names(bfl) <- c("ex1", "ex2", "input") rp <- readParam( pe = "none", dedup = TRUE, restrict = "chr10" ) wincounts <- windowCounts( bam.files = bfl, spacing = 60, width = 180, ext = 200, filter = 1, param = rp ) ## ----add-totals--------------------------------------------------------------- wincounts$totals <- c(964076L, 989543L, 1172179L) ## ----update-coldata----------------------------------------------------------- wincounts$sample <- colnames(wincounts) wincounts$treat <- as.factor(c("ctrl", "treat", NA)) colData(wincounts) ## ----plot-densities----------------------------------------------------------- library(extraChIPs) plotAssayDensities(wincounts, colour = "treat", trans = "log1p") ## ----peaks-------------------------------------------------------------------- peaks <- import.bed( system.file("extdata", "peaks.bed.gz", package = "extraChIPs") ) peaks <- granges(peaks) ## ----filtcounts--------------------------------------------------------------- filtcounts <- dualFilter( x = wincounts[, !is.na(wincounts$treat)], bg = wincounts[, is.na(wincounts$treat)], ref = peaks, q = 0.8 # Better to use q = 0.5 on real data ) ## ----plotcpm------------------------------------------------------------------ plotAssayDensities(filtcounts, assay = "logCPM", colour = "treat") plotAssayPCA(filtcounts, assay = "logCPM", colour = "treat", label = "sample") ## ----dims--------------------------------------------------------------------- dim(wincounts) dim(filtcounts) ## ----filt-ranges-------------------------------------------------------------- rowRanges(filtcounts) mean(rowRanges(filtcounts)$overlaps_ref) ## ----voom--------------------------------------------------------------------- v <- voomWeightsFromCPM( cpm = assay(filtcounts, "logCPM"), lib.size = filtcounts$totals, isLogCPM = TRUE ) ## ----add-vals----------------------------------------------------------------- rowRanges(filtcounts)$logCPM <- rowMeans(assay(filtcounts,"logCPM")) rowRanges(filtcounts)$logFC <- rowDiffs(assay(filtcounts,"logCPM"))[,1] rowRanges(filtcounts)$PValue <- 1 - pchisq(rowRanges(filtcounts)$logFC^2, 1) ## ----add-ol------------------------------------------------------------------- res_gr <- mergeByCol(filtcounts, col = "logCPM", pval = "PValue") res_gr$overlaps_ref <- overlapsAny(res_gr, peaks) ## ----ex-mapping1-------------------------------------------------------------- data("ex_genes") data("ex_prom") mapByFeature( res_gr, genes = ex_genes, prom = ex_prom, cols = c("gene", "symbol") ) ## ----ex-mapping2-------------------------------------------------------------- data("ex_hic") res_gr_mapped <- mapByFeature( res_gr, genes = ex_genes, prom = ex_prom, gi = ex_hic, cols = c("gene", "symbol") ) res_gr_mapped ## ----plot-pie----------------------------------------------------------------- res_gr %>% as_tibble() %>% plotPie("overlaps_ref") ## ----plot-dual-pie------------------------------------------------------------ res_gr$Feature <- bestOverlap( res_gr, GRangesList(Promoter = ex_prom), missing = "None" ) res_gr %>% as_tibble() %>% plotPie(x = "Feature", fill = "overlaps_ref") ## ----bwfl--------------------------------------------------------------------- bwfl <- system.file( "extdata", "bigwig", c("ex1.bw", "ex2.bw"), package = "extraChIPs" ) %>% BigWigFileList() %>% setNames(c("ex1", "ex2")) ## ----get-profile-------------------------------------------------------------- pd <- getProfileData(bwfl, res_gr) pd ## ----profile-heatmap---------------------------------------------------------- plotProfileHeatmap(pd, "profile_data") + scale_fill_viridis_c() + labs(fill = "CPM") + theme_bw() ## ----centred-heatmap---------------------------------------------------------- pd <- getProfileData(bwfl, colToRanges(res_gr, "keyval_range")) plotProfileHeatmap(pd, "profile_data") + scale_fill_viridis_c() + labs(fill = "CPM") + theme_bw() ## ----facet-heatmap------------------------------------------------------------ plotProfileHeatmap( pd, "profile_data", facetY = "overlaps_ref", linetype = "overlaps_ref" ) + scale_fill_viridis_c() + scale_colour_manual(values = c("red", "black")) + labs(fill = "CPM") + theme_bw() ## ----plot-empty-hfgc---------------------------------------------------------- data("grch37.cytobands") gr <- range(res_gr) plotHFGC(gr, cytobands = grch37.cytobands) ## ----add-genes---------------------------------------------------------------- data("ex_trans") plotHFGC(gr, genes = ex_trans, cytobands = grch37.cytobands) ## ----plot-trans--------------------------------------------------------------- plotHFGC( gr, genes = ex_trans, genecol = "wheat", collapseTranscripts = FALSE, cytobands = grch37.cytobands, zoom = 1.2 ) ## ----split-trans-------------------------------------------------------------- status_trans <- splitAsList(ex_trans, ex_trans$status) plotHFGC( gr, genes = status_trans, genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"), cytobands = grch37.cytobands, zoom = 1.2 ) ## ----plot-hfgc---------------------------------------------------------------- data("ex_prom") feat_grl <- GRangesList(Promoters = ex_prom) plotHFGC( gr, features = feat_grl, featcol = list(Promoters = "red"), genes = status_trans, genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"), cytobands = grch37.cytobands, zoom = 1.2 ) ## ----plot-with-hic------------------------------------------------------------ plotHFGC( gr, hic = ex_hic, features = feat_grl, featcol = list(Promoters = "red"), genes = status_trans, genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"), cytobands = grch37.cytobands, zoom = 1.2 ) ## ----plot-with-coverage------------------------------------------------------- plotHFGC( gr, hic = ex_hic, features = feat_grl, featcol = list(Promoters = "red"), genes = status_trans, coverage = bwfl, linecol = c(ex1 = "#4B0055", ex2 = "#007094"), genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"), cytobands = grch37.cytobands, zoom = 1.2 ) ## ----plot-with-tracks--------------------------------------------------------- cov_list <- list(TF1 = bwfl) plotHFGC( gr, hic = ex_hic, features = feat_grl, featcol = list(Promoters = "red"), genes = status_trans, coverage = cov_list, linecol = list(TF1 = c(ex1 = "#4B0055", ex2 = "#007094")), genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"), cytobands = grch37.cytobands, zoom = 1.2 ) ## ----cov-annot---------------------------------------------------------------- cov_annot <- splitAsList(res_gr, res_gr$logFC < -1) %>% setNames(c("Unchanged", "Decreased")) %>% endoapply(granges) ## ----plot-annot--------------------------------------------------------------- plotHFGC( gr, hic = ex_hic, features = feat_grl, featcol = list(Promoters = "red"), genes = status_trans, coverage = cov_list, annotation = list(TF1 = cov_annot), annotcol = c(Unchanged = "grey", Decreased = "#3333CC"), linecol = list(TF1 = c(ex1 = "#4B0055", ex2 = "#007094")), genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"), collapseTranscripts = "meta", cytobands = grch37.cytobands, zoom = 1.2 ) ## ----------------------------------------------------------------------------- sessionInfo()