library(RColorBrewer)
library(limma)
library(edgeR)

setwd("~/Projects/2014-04-08_GBP_pig/Analysis/")

counttable_all <- read.table("../RawData/featurecounts.merge", header=T, row.names=1)

timepoints = rep(c(0,15,30,60,90),times=3)
design <- model.matrix(~timepoints)

counts <- counttable_all[,1:15]
dge <- DGEList(counts=counts)
dge <- calcNormFactors(dge)
dge <- dge[ rowSums(cpm(dge) > 1)>2, ]

v <- voom(dge,design,plot=TRUE)

fit  <- lmFit(v, design)
fit <- eBayes(fit)

tt_genes_full<-topTable(fit,coef="timepoints",p.value=1,n=200000,sort="p")

table(tt_genes_full$adj.P.Val < .05)

## 
## FALSE 
## 14191

par(mfrow=c(3,3))

for (gene in rownames(tt_genes_full)[1:9]) {

  boxplot(v$E[gene,]~timepoints,ylab=gene,col=brewer.pal(length(unique(timepoints)),"Accent"),
          main=paste("Adj. p-val:",prettyNum(tt_genes_full[gene,]$adj.P.Val)),cex=2)

}

sessionInfo()

## R version 3.0.2 (2013-09-25)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] edgeR_3.4.2        limma_3.18.13      RColorBrewer_1.0-5
## [4] knitr_1.5.21      
## 
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.3 formatR_0.10.4 stringr_0.6.2  tools_3.0.2