ptm <- proc.time()
library(DEXSeq)
annotationfile <- gff
targets <- read.table(target.file, sep="\t", header=T)
samples.tmp <- strsplit(as.character(targets$Target), "_")
samples <- sapply(samples.tmp, function(x)x[1]) 
type <- sapply(samples.tmp, function(x)x[2]) 
type <- type[!is.na(type)]
counts.files <- sapply(as.character(targets$Name), function(x){
	file.path(getwd(), x, "accepted_hits_counts.txt")
})
counts.files <- as.character(counts.files)
ecs <- read.HTSeqCounts(countfiles = counts.files, design = samples, flattenedfile = annotationfile)
sampleNames(ecs) <- as.character(targets$Name)
gene.table <- geneCountTable(ecs)
tmp.sum <- apply(gene.table, 1,sum) 
gene.table <- gene.table[which(tmp.sum != 0),]
save(gene.table, file="time_genetable.rda")
ecs <- estimateSizeFactors(ecs)
sizeFactors(ecs)
ecs <- estimateDispersions(ecs)
ecs <- fitDispersionFunction(ecs)
ecs <- testForDEU(ecs)
ecs <- estimatelog2FoldChanges(ecs)
res.exons <- DEUresultTable(ecs)
save(res.exons, file="time_resexons.rda")
save(ecs, file="time_ecs.rda")
proc.time() - ptm
q(save="no")
