XBSeq.R
In XBSeq: Test for differential expression for RNA-seq data

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

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("XBSeq")

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

## ----eval=FALSE---------------------------------------------------------------
#  library(devtools)
#  install_github('liuy12/XBSeq')

## ----eval=FALSE---------------------------------------------------------------
#  fc_signal <- featureCounts(files = bamLists, annot.ext = gtf_file, isGTFAnnotationFile = TRUE)
#  fc_bg <- featureCounts(files = bamLists, annot.ext = gtf_file_bg, isGTFAnnotationFile = TRUE)

## ----eval=FALSE---------------------------------------------------------------
#  features_signal <- import(gtf_file)
#  features_signal <- split(features_signal, mcols(features_signal)$gene_id)
#  so_signal <- summarizeOverlaps(features_signal, bamLists)
#  
#  ## for background noise
#  features_bg <- import(gtf_file_bg)
#  features_bg <- split(features_bg, mcols(features_bg)$gene_id)
#  so_bg <- summarizeOverlaps(features_bg, bamLists)

## ----eval=FALSE---------------------------------------------------------------
#  apaStats <- apaUsage(bamTreatment, bamControl, apaAnno)

## -----------------------------------------------------------------------------
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---------------------------------------------------------------
#  BiocManager::install("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()