inst/batch/batch_dexseq_exon.R

In oneChannelGUI: A graphical interface designed to facilitate analysis of microarrays and miRNA/RNA-seq data on laptops

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_dex_counts.txt")
})
counts.files <- as.character(counts.files)
sampleTable <- cbind(counts.files,as.character(samples))
sampleTable <- as.data.frame(sampleTable)
names(sampleTable) <- c("samples","condition")
ecs <- DEXSeqDataSetFromHTSeq(countfiles = as.character(sampleTable$samples), sampleData=sampleTable, flattenedfile = annotationfile)
colData(ecs)$samples <- as.character(targets$Name)
ecs <- estimateSizeFactors(ecs)
sizeFactors(ecs)
ecs <- estimateDispersions(ecs)
#ecs <- fitDispersionFunction(ecs)
ecs <- testForDEU(ecs)
ecs <- estimateExonFoldChanges(ecs, fitExpToVar="condition")
res.exons <- DEXSeqResults(ecs)
res.exons <- res.exons[!is.na(res.exons$padj),]
save(res.exons, file="resexons_time.rda")
save(ecs, file="ecs_time.rda")
proc.time() - ptm
q(save="no")