###################################################
### chunk number 1: countsTable
###################################################
library("DESeq")
countsTable <- read.delim(system.file("extra/TagSeqExample.tab", package="DESeq"), 
                          header=TRUE, stringsAsFactors=TRUE, row.names="gene")
conds <- factor(c("T", "T", "T", "T2", "N", "N"))


###################################################
### chunk number 2: countsTable
###################################################
head(countsTable)
conds


###################################################
### chunk number 3: demo1
###################################################
cds <- newCountDataSet( countsTable, conds )
cds <- estimateSizeFactors( cds )
cds <- estimateVarianceFunctions( cds )
res <- nbinomTest( cds, "T", "N")


###################################################
### chunk number 4: res
###################################################
head(res)


###################################################
### chunk number 5: 
###################################################
library( DESeq )
exampleFile = system.file( "extra/TagSeqExample.tab", package="DESeq" )
exampleFile


###################################################
### chunk number 6: 
###################################################
countsTable <- read.delim( exampleFile, header=TRUE, stringsAsFactors=TRUE )
head( countsTable )


###################################################
### chunk number 7: 
###################################################
rownames( countsTable ) <- countsTable$gene
countsTable <- countsTable[ , -1 ]


###################################################
### chunk number 8: 
###################################################
conds <- c( "T", "T", "T", "T2", "N", "N" )


###################################################
### chunk number 9: 
###################################################
cds <- newCountDataSet( countsTable, conds )


###################################################
### chunk number 10: 
###################################################
head( counts(cds) )


###################################################
### chunk number 11: 
###################################################
cds <- cds[ ,-1 ]


###################################################
### chunk number 12: 
###################################################
libsizes <- c( T1a=6843583, T1b=7604834, T2=13625570, T3=12291910,
   N1=12872125, N2=10502656 )
sizeFactors(cds) <- libsizes[-1]


###################################################
### chunk number 13: 
###################################################
cds <- estimateSizeFactors( cds )
sizeFactors( cds )


###################################################
### chunk number 14: 
###################################################
cds <- estimateVarianceFunctions( cds )


###################################################
### chunk number 15: 
###################################################
countValue <- 123
baseLevel <- countValue / sizeFactors(cds)["T1b"]
rawVarFuncForGB <- rawVarFunc( cds, "T" )
rawVariance <- rawVarFuncForGB( baseLevel )
fullVariance <- countValue + rawVariance * sizeFactors(cds)["T1b"]^2
sqrt( fullVariance )


###################################################
### chunk number 16: figSCV
###################################################
scvPlot( cds )


###################################################
### chunk number 17: 
###################################################
diagForT <- varianceFitDiagnostics( cds, "T" )
head( diagForT )


###################################################
### chunk number 18: figFit
###################################################
smoothScatter( log10(diagForT$baseMean), log10(diagForT$baseVar) )
lines( log10(fittedBaseVar) ~ log10(baseMean), 
   diagForT[ order(diagForT$baseMean), ], col="red" )


###################################################
### chunk number 19: figECDF
###################################################
par( mfrow=c(1,2 ) )
residualsEcdfPlot( cds, "T" )
residualsEcdfPlot( cds, "N" )


###################################################
### chunk number 20: 
###################################################
rawVarFuncTable( cds )


###################################################
### chunk number 21: 
###################################################
res <- nbinomTest( cds, "N", "T" )
head(res)


###################################################
### chunk number 22: figDE
###################################################
plot( 
   res$baseMean, 
   res$log2FoldChange, 
   log="x", pch=20, cex=.1, 
   col = ifelse( res$padj < .1, "red", "black" ) )


###################################################
### chunk number 23: 
###################################################
resSig <- res[ res$padj < .1, ]


###################################################
### chunk number 24: 
###################################################
head( resSig[ order(resSig$pval), ] )


###################################################
### chunk number 25: 
###################################################
head( resSig[ order( resSig$foldChange, -resSig$baseMean ), ] )


###################################################
### chunk number 26: 
###################################################
head( resSig[ order( -resSig$foldChange, -resSig$baseMean ), ] )


###################################################
### chunk number 27: 
###################################################
resT2vsN <- nbinomTest( cds, "N", "T2" )


###################################################
### chunk number 28: figDE_T2
###################################################
plot( 
   resT2vsN$baseMean, 
   resT2vsN$log2FoldChange, 
   log="x", pch=20, cex=.1, 
   col = ifelse( resT2vsN$padj < .1, "red", "black" ) )


###################################################
### chunk number 29: 
###################################################
cds2 <- cds[ ,c( "T1b", "N1" ) ]


###################################################
### chunk number 30: 
###################################################
cds2 <- estimateVarianceFunctions( cds2, pool=TRUE )


###################################################
### chunk number 31: 
###################################################
res2 <- nbinomTest( cds2, "N", "T" )


###################################################
### chunk number 32: figDE2
###################################################
plot( 
   res2$baseMean,
   res2$log2FoldChange, 
   log="x", pch=20, cex=.1, 
   col = ifelse( res2$padj < .1, "red", "black" ) )


###################################################
### chunk number 33: 
###################################################
table( res_sig = res$padj < .1, res2_sig = res2$padj < .1 )


###################################################
### chunk number 34: 
###################################################
colsN <- conditions(cds) == "N"
colsT <- conditions(cds) == "T"
baseMeansNT <- getBaseMeansAndVariances( 
   counts(cds)[ , colsN|colsT ], 
   sizeFactors(cds)[ colsN|colsT ] )$baseMean
pvals2b <- nbinomTestForMatrices(
   counts(cds)[ ,colsN ],
   counts(cds)[ ,colsT ],
   sizeFactors(cds)[ colsN ],
   sizeFactors(cds)[ colsT ],
   rawVarFunc( cds2, "_pooled" )( baseMeansNT ),
   rawVarFunc( cds2, "_pooled" )( baseMeansNT ) )


###################################################
### chunk number 35: 
###################################################
padj2b <- p.adjust( pvals2b, method="BH" )
table( res_sig = res$padj < .1, res2b_sig = padj2b < .1 )


###################################################
### chunk number 36: 
###################################################
cds3 <- newCountDataSet( countsTable, conds )
cds3 <- estimateSizeFactors( cds3 )
cds3 <- estimateVarianceFunctions( cds3 )


###################################################
### chunk number 37: 
###################################################
vsd <- getVarianceStabilizedData( cds3 )


###################################################
### chunk number 38: 
###################################################
head( vsd )


###################################################
### chunk number 39: 
###################################################
dists <- dist( t( vsd ) )


###################################################
### chunk number 40: figHeatmap
###################################################
heatmap( as.matrix( dists ), symm=TRUE )


###################################################
### chunk number 41: 
###################################################
sessionInfo()