## ----style, echo = FALSE, results = 'asis'-------------------------------
BiocStyle::markdown()

## ----,eval=FALSE---------------------------------------------------------
## source('http://www.bioconductor.org/biocLite.R')
## biocLite("XBSeq")

## ----,message = FALSE, warning=FALSE-------------------------------------
library("XBSeq")

## ----,engine='python',eval=FALSE-----------------------------------------
## htseq-count [options] <alignment_file> <gtf_file> > Observed_count.txt
## htseq-count [options] <alignment_file> <gtf_file_bg> > background_count.txt

## ------------------------------------------------------------------------
data(ExampleData)

## ------------------------------------------------------------------------
head(Observed)
head(Background)

## ----,tidy=TRUE----------------------------------------------------------
conditions <- factor(c(rep('C',3), rep('T', 3)))
XB <- XBSeqDataSet(Observed, Background, conditions)

## ----,tidy=TRUE,fig.width=5,fig.height=4---------------------------------
XBplot(XB, Samplenum = 1, unit = "LogTPM", Genelength = genelength[,2])

## ----, tidy=TRUE---------------------------------------------------------
XB <- estimateRealCount(XB)
XB <- estimateSizeFactors(XB)
XB <- estimateSCV( XB, method='pooled', sharingMode='maximum', fitType='local' )

## ----,fig.width=3,fig.height=3-------------------------------------------
plotSCVEsts(XB)

## ------------------------------------------------------------------------
Teststas <- XBSeqTest( XB, levels(conditions)[1L], levels(conditions)[2L], method ='NP')

## ----,fig.width=3,fig.height=3-------------------------------------------
MAplot(Teststas, padj = FALSE, pcuff = 0.01, lfccuff = 1)

## ----,eval=FALSE,tidy=TRUE-----------------------------------------------
## # Alternatively, all the codes above can be done with a wrapper function XBSeq
## Teststats <- XBSeq( Observed, Background, conditions, method='pooled', sharingMode='maximum',
##   fitType='local', pvals_only=FALSE, paraMethod = 'NP' )

## ----,eval=FALSE---------------------------------------------------------
## biocLite("DESeq")

## ----,message=FALSE------------------------------------------------------
library('DESeq')
library('ggplot2')
de <- newCountDataSet(Observed, conditions)
de <- estimateSizeFactors(de)
de <- estimateDispersions(de, method = "pooled", fitType="local")
res <- nbinomTest(de, levels(conditions)[1], levels(conditions)[2])

## ----,warning=FALSE,message=FALSE,tidy=TRUE, fig.width=3,fig.height=3----
DE_index_DESeq <- with(res, which(pval<0.01 & abs(log2FoldChange)>1))
DE_index_XBSeq <- with(Teststas, which(pval<0.01 & abs(log2FoldChange)>1))
DE_index_inters <- intersect(DE_index_DESeq, DE_index_XBSeq)
DE_index_DESeq_uniq <- setdiff(DE_index_DESeq, DE_index_XBSeq)
DE_plot <- MAplot(Teststas, padj = FALSE, pcuff = 0.01, lfccuff = 1, shape=16)
DE_plot + geom_point( data=Teststas[DE_index_inters,], aes(x=baseMean, y=log2FoldChange),
                      color= 'green', shape=16 ) + 
  geom_point( data=Teststas[DE_index_DESeq_uniq,], aes( x=baseMean, y=log2FoldChange ),
              color= 'blue', shape=16 )

## ------------------------------------------------------------------------
sessionInfo()