DEseq_and_heatmap.Rmd

---
title: "DEseq_and_heatmap"
author: "Chengqi(Charley) Wang"
date: "2/1/2020"
output: html_document
---
<br/><br/>

####load the txt count data
```{r}
df      <- read.table('~/Documents/work/paper_write/genome_training/R_training/vehicle_drug_feature_counts.txt',
                 header = T, sep = '\t', row.names = 1)

df      <- df[,6:9]
```
<br/><br/>

####check the distribution
```{r}
par(mfrow = c(2,2), mar = c(4,4,1,1))
print( apply(df, 2, function(x) quantile(as.numeric(x))) )
log1 <- apply(df, 2, function(x) {hist(log2(x), breaks = 100)})
```
<br/><br/>

#####Pick up the gene with expression bigger than 2^4
```{r}
idExp <- apply(df, 1, function(x) any(x > 16))
dfExp <- df[idExp, ]
print(dim(dfExp))
```
<br/><br/>

####build DESeq class
```{r, message = F}
library(DESeq)
```
```{r}
condition<- c('c','c', 't', 't')
cds      <- newCountDataSet( dfExp, condition )
```

<br/><br/>

####Differential expression gene calling
```{r}
##normalize the data, extract the size factor for each assay
cds      <- estimateSizeFactors(cds)
head( counts( cds, normalized=TRUE ) )
##estimate the background gene expression 
cds      <- estimateDispersions( cds )

##estimate the background distribution
plotDispEsts(cds)
##p-value fit
res        <- nbinomTest( cds, 'c', 't' )
plotMA(res)
```
<br/><br/>

####Differential expression gene extraction and visulation
```{r}
## extract the genes with adjP < 0.0001
resSig = res[ which(res$padj < 0.0001), ]
dim(resSig)

norData <- counts( cds, normalized=TRUE )
sigNorData <- norData[which(res$padj < 0.0001),]
save(sigNorData, resSig,
     file = 'deseq.sig.R')


##extract normalized data
norDF           <- counts( cds, normalized=TRUE )
heatMapDf       <- norDF[(rownames(norDF) %in% resSig$id), ]

##heatmap plot
##plot heatmap
library("RColorBrewer")
library("gplots")
hmcol =  colorRampPalette(brewer.pal(9, "GnBu"))(100)
## expression heatmap
heatMapDF_nor <- t( apply(heatMapDf, 1, function(x){(x-mean(x))/sd(x)}) )

colnames(heatMapDF_nor) <- c('control1','control2',
                             'treat1'  , 'treat2')
#pdf('aaaaa.pdf', height = 10, width = 10)
heatmap.2(heatMapDF_nor, col = hmcol, trace="none", margin=c(6, 6),labRow = F)
#dev.off()
```