figure_3_b_twi.Rmd

``` {r setup, echo=FALSE, message=FALSE, include=FALSE, error=FALSE}
library(xtable)
library(ggplot2)
library(plyr)
library(BSgenome.Dmelanogaster.UCSC.dm3)
library(parallel)
library(rtracklayer)

# Output folder for this document
options(knitr.figure_dir = "figure_3_b_twi_output")

source("shared_code/knitr_common.r")
source("shared_code/granges_common.r")
source("shared_code/samples.r")
source("shared_code/exo_metapeak.r")
source("shared_code/profiles_common.r")
source("shared_code/ggplot_theme.r")
```

# Figure 3B: Twist ChIP-nexus

``` {r header_child, child="child_docs/header_child.Rmd"}
```

## Overview

We will compare the profiles of the following samples around the top 200 Twist `CABATG` motifs:

``` {r max_samples_table, results="asis"}
twi.df <- subset(samples.df, sample == "dmel_embryo_twi_chipnexus_01")
html_table(twi.df)
```

## Motif

E-BOX motif: `CABATG`

Window size centered at motif: 51 bp

``` {r calculate_profile_reads, include=FALSE}

reads.list <- cache("reads.list.rds", function() {
  twi.motifs <- c("CACATG", "CATGTG", "CAGATG", "CATCTG", "CATATG")
  twi.gr <- filter_chrs(do.call(c, lapply(twi.motifs, vmatchPattern, Dmelanogaster, max.mismatch=0, fixed=TRUE)))
  twi.gr <- twi.gr[strand(twi.gr) == "+"]

  checked_mclapply(twi.df$sample, process_sample, twi.gr, n=200, region_size=51, mc.cores=2)
})

```

## Average profile (top 200)

``` {r plots_per_sample, fig.cap="", fig.width=9, fig.height=6, dpi=100}
plots.list <- lapply(reads.list, build_plot)
nothing <- lapply(plots.list, print)
```

``` {r zoomed_plot, warning=FALSE, fig.cap="", fig.width=9, fig.height=6, dpi=100}
reads.df <- reads.list[[1]]

motif.box <- data.frame(xmin=0, 
                        xmax=reads.df$motif_width[1]-1,
                        ymin=-Inf,
                        ymax=Inf)
                          
g <- ggplot(reads.df, aes(x=tss_distance, y=reads, color=strand)) +
     geom_line(size=1.2) +
     geom_vline(xintercept=0:5, color="gray50") +
     scale_colour_manual(name="Strand", values=c("+"="red", "-"="darkblue")) +
     geom_rect(show_guide=FALSE, inherit.aes=FALSE, data=motif.box, 
               aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
               alpha=0.25, fill="gray80") +
     theme_manuscript() +
     scale_x_continuous(breaks=c(-20, -15, -10, -5, 0, 5, 10, 15, 20, 25), 
                        labels=c(-20, -15, -10, -5, 0, 0,  5, 10, 15, 20),
                        limits=c(-30, 35)) +
     labs(x="Distance to motif edge", y="Average ChIP-nexus reads", title="CABATG")
g
```

``` {r save_zoomed_pdf, include=FALSE}
pdf(figure_path("twi_CABATG_profile.pdf"), width=9, height=6)
print(g)
dev.off()
```

## Heatmap

``` {r collect_heatmap_reads, include=FALSE}

beds <- list.files(figure_path(), "bed")
names(beds) <- gsub("_top_200_motifs.bed", "", beds)

sample.list <- mclapply(beds, heatmap_reads, mc.cores=4)
```

``` {r plot_heatmaps, fig.cap="", fig.width=6, fig.heigth=6, dpi=100}
nothing <- lapply(names(sample.list), function(n) { draw_exo_heatmap(sample.list[[n]], n, 75:125)})
```

``` {r create_heatmap_pdf, include=FALSE}
pdf(figure_path("twi_heatmaps.pdf"), width=6, height=6)
nothing <- lapply(names(sample.list), function(n) { draw_exo_heatmap(sample.list[[n]], n, 75:125)})
dev.off()
```


``` {r session_info_child, child="child_docs/session_info_child.Rmd"}
```