DEXseq_example.R

library(DEXSeq)
library(pasilla)

setwd("/path/to/dir")

# Doing the vignette
# Get locations for count files and flattened annotation file in gff formar
inDir = "~/Dropbox/Jacobs_project/relative_exon_usage/exon_counts"
countFiles = list.files(inDir, pattern="exon.txt$", full.names=TRUE)
countFiles
flattenedFile = list.files(inDir, pattern="gff$", full.names=TRUE)
flattenedFile

# Prepare sample table
sampleTable = data.frame(
  row.names = c( "Control_A", "Control_B", "Control_C", 
                 "HF_A", "HF_B", "HF_C" ),
  condition = c("Control", "Control", "Control",  
                "HF", "HF", "HF" ))
  
# Construct DEXSeqDataSet object
dxd = DEXSeqDataSetFromHTSeq(
  countFiles,
  sampleData=sampleTable,
  design= ~ sample + exon + condition:exon,
  flattenedfile=flattenedFile )
dxd

# estimate size factors
dxd = estimateSizeFactors( dxd )

# estimate dispersions
dxd = estimateDispersions( dxd )
par(mar = rep(2, 4))
plotDispEsts( dxd )

# statistical testing
dxd <- testForDEU(dxd)

# estimate exon fold changes
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")

# get results
dxr1 <- DEXSeqResults(dxd)
write.table(dxr1, file="DEXSeq_results.txt", sep="\t", quote=F, col.names = T, row.names = T)
table(dxr1$padj < 0.1)
table ( tapply( dxr1$padj < 0.1, dxr1$groupID, any ) )
png("MA_plot.png", width=800, height = 800)
plotMA(dxr1, cex = 0.8)
dev.off()
sigOnly <- subset(dxr1, padj < 0.1) # select significant differential exon usage events
write.table(sigOnly, file="significant_exon_usage_fdr0.1.txt", sep="\t", quote=F, col.names = T, row.names = T)

# visualize the results differential exon useage events
sigGenes <- rownames(sigOnly)
# select gene ids of significant genes
sigGenes <- sapply(strsplit(sigGenes, split=':', fixed=T), "[[", 1)

plotExonsSingleGene <- function(geneList, DEXSeqRes)
{
  for (i in geneList) 
  {
      pdf(paste(i, "_exon_structure.pdf", sep=""))
      plotDEXSeq( DEXSeqRes, i, displayTranscripts = TRUE,
            legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )
      dev.off()
  }
}

# Build and save differential exon usage figures
setwd("diff_exon_usage_figures/")
plotExonsSingleGene(geneList = sigGenes, DEXSeqRes = dxr1)

par(mfrow=c(1,2))