lncRNA-heatmap.rmd

---
title: "lncRNA_heatmap"
author: "Ming Tang"
date: "July 28, 2016"
output: html_document
---

Read in the bigwig files for each mark. bigwig files were generated by Deeptools from bam files.
```{r}
library(EnrichedHeatmap)
library(rtracklayer)
library(GenomicFeatures)


TCGA.27.1831.H3K27me3<- import("/Users/mtang1/projects/GBM_ChIP-seq/results/GBM_normalized_bigwigs/TCGA-27-1831-K27me3.bw", format="bigWig")

TCGA.27.1831.H3K79me2<-import("/Users/mtang1/projects/GBM_ChIP-seq/results/GBM_normalized_bigwigs/TCGA-27-1831-K79me2.bw", format="bigWig")

normal.brain.mean.H3K27me3<- import("/Users/mtang1/projects/GBM_ChIP-seq/results/GBM_normalized_bigwigs/normalBrain-mean-K27me3.bw", format="bigWig")
```


```{r}
lncRNA<- readRDS("~/Desktop/tier1_lncs_extended_annotations_gencode_v19.rds")

lncRNA.gr<- GRanges(seqnames=lncRNA$seqid, ranges=IRanges(start=lncRNA$start, end=lncRNA$end), strand=lncRNA$strand)
lncRNA.gr$gene_name<- lncRNA$gene_name
```

change to UCSC chromosome format because the bigwig chromosome names are with "chr"

```{r}
library(BSgenome.Hsapiens.UCSC.hg19)
seqlevelsStyle(lncRNA.gr) <- "UCSC"
seqlevels(lncRNA.gr)<- seqlevels(BSgenome.Hsapiens.UCSC.hg19)
seqlengths(lncRNA.gr)<- seqlengths(BSgenome.Hsapiens.UCSC.hg19)

## get the 1bp of the TSS 

lncRNA.TSS<- resize(lncRNA.gr, width = 1)

``` 

Calculate the matrix, extend TSS to upstream and downstream 10kb, and use 100bp as bin size
```{r}
lncRNA.H3K27me3.TSS.mat<- normalizeToMatrix(TCGA.27.1831.H3K27me3, lncRNA.TSS, value_column = "score", extend=10000, mean_mode="w0", w=100)

## H3K27me2 is enriched in active gene body, use lncRNA whole gene body as target
lncRNA.H3K79me2.TSS.mat<- normalizeToMatrix(TCGA.27.1831.H3K79me2, lncRNA.gr, value_column = "score", extend=10000,
                                            mean_mode="w0", w=100, target_ratio = 0.5)
## normal brain 
lncRNA.normal.H3K27me3.TSS.mat<- normalizeToMatrix(normal.brain.mean.H3K27me3, lncRNA.TSS, value_column = "score", extend=10000, mean_mode="w0", w=100)

```

To visualize the matrix, check the range of the the data first, and use a color mapping function for robust visulization (to outliers)

```{r}
quantile(lncRNA.H3K27me3.TSS.mat, probs = c(0.005, 0.5,0.995))
col_fun_H3K27me3<- circlize::colorRamp2(c(0,15), c("white", "red"))

quantile(lncRNA.H3K79me2.TSS.mat, probs = c(0.005, 0.5,0.995))
col_fun_H3K79me3<- circlize::colorRamp2(c(0,27), c("white", "red"))

quantile(lncRNA.normal.H3K27me3.TSS.mat, probs = c(0.005, 0.5,0.995))
col_fun_H3K27me3.normal<- circlize::colorRamp2(c(0,24), c("white", "red"))

## draw heatmap, try kmeans clustering for rows, or you can disable it
lncRNA.ht.H3K79me2<-EnrichedHeatmap(lncRNA.H3K79me2.TSS.mat, col= col_fun_H3K79me3, name = "H3K79me2 at lncRNA gene body",
                axis_name_rot = 0, use_raster = FALSE, km=2, column_title = "H3K79me2")

lncRNA.ht.H3K27me3<-EnrichedHeatmap(lncRNA.H3K27me3.TSS.mat, col= col_fun_H3K27me3, name = "H3K27me3 at lncRNA TSS",
                axis_name_rot = 0, use_raster = FALSE, column_title = "H3K27me3" ) 

lncRNA.ht.H3K27me3.normal<- EnrichedHeatmap(lncRNA.normal.H3K27me3.TSS.mat, col= col_fun_H3K27me3.normal, 
                                             name = "H3K27me3 at lncRNA TSS for normal",
                axis_name_rot = 0, use_raster = FALSE)  

## draw heatmap separately   
draw(lncRNA.ht.H3K27me3)
draw(lncRNA.ht.H3K79me2)
draw(lncRNA.ht.H3K27me3.normal)

## or draw them together and set H3K79me2 as the main heatmap, and the other heatmap will split rows according to the main heatmap

draw(lncRNA.ht.H3K79me2 + lncRNA.ht.H3K27me3 )
```