```{r init, include=FALSE}
library(knitr)
library(BiocStyle)
opts_chunk$set(message=FALSE)
```

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```
# Introduction

## K-562

From Wikipedia:

> K562 cells were the first human immortalised myelogenous leukemia cell line to be established. K562 cells are of the erythroleukemia type, and the cell line is derived from a 53-year-old female chronic myelogenous leukemia patient in blast crisis. The cells are non-adherent and rounded, are positive for the bcr:abl fusion gene, and bear some proteomic resemblance to both undifferentiated granulocytes and erythrocytes. In culture they exhibit much less clumping than many other suspension lines, presumably due to the downregulation of surface adhesion molecules by bcr:abl. However, another study suggests that bcr:abl over-expression may actually increase cell adherence to cell culture plastic. K562 cells can spontaneously develop characteristics similar to early-stage erythrocytes, granulocytes and monocytes and are easily killed by natural killer cells[8] as they lack the MHC complex required to inhibit NK activity. They also lack any trace of Epstein-Barr virus and other herpesviruses. In addition to the Philadelphia chromosome they also exhibit a second reciprocal translocation between the long arm of chromosome 15 with chromosome 17. 

[Cellosaurus](https://web.expasy.org/cellosaurus/CVCL_0004) is a great resource for 
finding resources related to cell lines, including derivative lines with specific 
induced mutations.


# Gene Expression Data

## Annotation data

```{r gencodegtf}
download.file("https://www.encodeproject.org/files/gencode.v24.primary_assembly.annotation/@@download/gencode.v24.primary_assembly.annotation.gtf.gz", "gencode.v24.primary_assembly.annotation.gtf.gz")
```

```{r txdb}
library(GenomicFeatures)
txdb = makeTxDbFromGFF('gencode.v24.primary_assembly.annotation.gtf.gz')
```



## Building a SummarizedExperiment

The ENCODE data portal contains links to both raw and processed data. We are going to be working
with relatively standard polyA-selected RNA-seq values. 

[Protocol]

```{r geneExpression}
download.file("https://www.encodeproject.org/files/ENCFF172GIN/@@download/ENCFF172GIN.tsv", "ENCFF172GIN.tsv")
download.file("https://www.encodeproject.org/files/ENCFF768TKT/@@download/ENCFF768TKT.tsv", "ENCFF768TKT.tsv")
```

[Protocol]: https://www.encodeproject.org/documents/6354169f-86f6-4b59-8322-141005ea44eb/@@download/attachment/Long%20RNA-seq%20pipeline%20overview.pdf


# tximport

```{r tximport, warning=FALSE}
library(tximport)
fnames = c("ENCFF172GIN.tsv", "ENCFF768TKT.tsv")
txi = tximport(fnames, txIn = FALSE, txOut = FALSE, type = 'rsem')
```

```{r semake}
library(SummarizedExperiment)
colnames(txi$abundance) = c('rep1','rep2')
colnames(txi$counts)    = c('rep1','rep2')
coldata = DataFrame(cell='K562',rep=c("rep1", "rep2"),row.names = c('rep1', 'rep2'))
se = SummarizedExperiment(list(counts = txi$counts, abundance = txi$abundance), colData = coldata)
```

Now, we limit the gene quantifications to those in our `txdb` object. There are
data for tRNAs and for ERCC control spikeins, but we will simply drop those here.

```{r limitserows, cache=FALSE}
library(GenomicFeatures)
genes_gr = genes(txdb)
common_genes = intersect(names(genes_gr), rownames(se))
se = se[common_genes,]
rowRanges(se) = genes_gr[rownames(se)]
```

# DNAse data

```{r annotationhub}
library(AnnotationHub)
ah = AnnotationHub()
query(ah, c('E123','dnase','narrowpeak'))
k562_dnase = ah[['AH30828']]
```


[narrowPeak format](https://genome.ucsc.edu/FAQ/FAQformat.html#format12)

Note that the k562_dnase includes a ranked list of peaks. The "top" peaks have the highest -log10(p-values)
and generally very high signalValue. We can subset to get the 


```{r fixseqlevels, cache=FALSE}
common_seqlevels = intersect(seqlevels(k562_dnase), seqlevels(se))
keepSeqlevels(se, common_seqlevels, pruning.mode = "coarse")
keepSeqlevels(k562_dnase, common_seqlevels, pruning.mode = "coarse")
ne = distanceToNearest(se, k562_dnase[1:20000])
regioncount = length(k562_dnase)
ne2 = distanceToNearest(se, k562_dnase[(regioncount-20000):regioncount])
```

The `ne` and `ne2` contain the distances to the nearest gene, but for two
distinct sets of peaks. Given that the `ne` variable is for "high confidence"
peaks and that the `ne2` contains distances to the nearest gene for "low confidence"
peaks, are there differences between the two sets of distances?

- Define a good set of peaks and a "bad" set. Look at the location of these peaks with
respect to "promoter" regions of the genes.

- What is the expression value of the genes that have a "good" peak in their promoter
region versus those without?

- How about for the "bad" peak genes?




