## ---- echo=FALSE, message=FALSE, warning=FALSE--------------------------- library(diffloop) library(diffloopdata) library(ggplot2) library(GenomicRanges) ## ----load_data----------------------------------------------------------- bed_dir <- system.file("extdata", "esc_jurkat", package="diffloopdata") bed_dir samples <- c("naive_esc_1", "naive_esc_2", "primed_esc_1", "primed_esc_2", "jurkat_1", "jurkat_2") ## ----read_data----------------------------------------------------------- full <- loopsMake(bed_dir, samples, 1500) celltypes <- c("naive", "naive", "primed", "primed", "jurkat", "jurkat") full <- updateLDGroups(full, celltypes) head(full, 3) ## ----pca_prenorm, fig.width=7, fig.height=4, fig.cap ="PCA plot of**un-normalized** ChIA-PET loop data"---- pcaPlot(full) + geom_text(aes(label=samples)) ## ----sizeFactors--------------------------------------------------------- full <- calcLDSizeFactors(full) head(full,3) ## ----qc------------------------------------------------------------------ loopMetrics(full) ## ----qc2----------------------------------------------------------------- full_5kb <- filterLoops(full, width=5000, nreplicates=3, nsamples = 1) dim(full_5kb) valid_full <- filterLoops(full, width=5000, nreplicates=3, nsamples = 2) dim(valid_full) ## ----qc3----------------------------------------------------------------- dim(intrachromosomal(valid_full)) dim(interchromosomal(valid_full)) ## ----qc4----------------------------------------------------------------- numLoops(valid_full,1:5) numAnchors(valid_full) ## ----pca_postnorm, fig.width=7, fig.height=4, fig.cap = "PCA plot of normalized ChIA-PET loop data"---- pcaPlot(valid_full) + geom_text(aes(label=samples)) + expand_limits(x = c(-20, 34)) ## ----geneAssoc----------------------------------------------------------- # Fit edgeR Models Pairwise between cell types fit <- loopFit(valid_full) jn_res <- loopTest(fit,coef=2) jp_res <- loopTest(fit,coef=3) np_res <- loopTest(fit,contrast=c(0,-1,1)) head(np_res) ## ----quickAssoc---------------------------------------------------------- np <- valid_full[,1:4] summary(np[1:3,]) np_res2 <- quickAssoc(np) ## ----assocResults-------------------------------------------------------- jnp_res <- loopTest(fit,coef=2:3) head(jnp_res) ## ----swAssoc------------------------------------------------------------- sw_jn <- slidingWindowTest(jn_res,200000,200000) head(sw_jn) ## ----featAssoc----------------------------------------------------------- chr1 <- getHumanGenes(c("1")) ft_jnp <- featureTest(jnp_res,chr1) head(ft_jnp,10) ## ----pairwiseSum, message=FALSE, warning=FALSE--------------------------- np_sig_res <- topLoops(np_res, FDR = 0.2) dim(np_sig_res) dim(topLoops(jp_res, FDR = 0.2)) dim(topLoops(jn_res, FDR = 0.2)) summary(np_sig_res) ## ----pca_sigfix, fig.width=7, fig.height=4, fig.cap = "PCA plot of differential loops"---- pcaPlot(topLoops(jnp_res, FDR = 0.05)) + geom_text(aes(label=samples)) ## ----dloplots2,fig.width=7, fig.height=10-------------------------------- regA <- GRanges(c("1"),ranges=IRanges(c(36000000),c(36300000))) full.plot <- loopPlot(jn_res, regA)