inst/doc/bayNorm.R

In bayNorm: Single-cell RNA sequencing data normalization

## ---- echo=FALSE, results="hide", message=FALSE-------------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE, dpi = 30)
knitr::opts_chunk$set(dev="png",fig.align="center")

## ----library, echo=FALSE------------------------------------------------------
library(bayNorm)
library(BiocStyle)
library(SummarizedExperiment)
library(SingleCellExperiment)

## ----install_Bioconductor, eval=FALSE-----------------------------------------
#  library(BiocManager)
#  BiocManager::install("bayNorm")

## ----install, eval=FALSE------------------------------------------------------
#  library(devtools)
#  devtools::install_github("WT215/bayNorm")

## ----beta, echo=TRUE,eval=TRUE------------------------------------------------
data('EXAMPLE_DATA_list')
#Suppose the input data is a SummarizedExperiment object:
#Here we just used 30 cells for illustration
rse <- SummarizedExperiment::SummarizedExperiment(assays=SimpleList(counts=EXAMPLE_DATA_list$inputdata[,seq(1,30)]))
#SingleCellExperiment object can also be input in bayNorm:
#rse <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=EXAMPLE_DATA_list$inputdata))

BETA_est<-BetaFun(Data=rse,MeanBETA=0.06)
summary(BETA_est$BETA)

## ----intro_bayNorm, echo=TRUE,eval=TRUE---------------------------------------
#Return 3D array normalzied data, draw 20 samples from posterior distribution:
bayNorm_3D<-bayNorm(
    Data=rse,
    BETA_vec = NULL,
    mode_version=FALSE,
    mean_version = FALSE,S=20
    ,verbose =FALSE,
    parallel = TRUE)

#Return 2D matrix normalized data (MAP of posterior):
#Simply set mode_version=TRUE, but keep mean_version=FALSE


#Return 2D matrix normalized data (mean of posterior):
#Simply set mean_version=TRUE, but keep mode_version=FALSE

## ----intro_bayNorm_sup, echo=TRUE,eval=TRUE-----------------------------------
#Now if you want to generate 2D matrix (MAP) using the same prior
#estimates as generated before:
bayNorm_2D<-bayNorm_sup(
    Data=rse,
    PRIORS=bayNorm_3D$PRIORS,
    input_params=bayNorm_3D$input_params,
    mode_version=TRUE,
    mean_version = FALSE,
    verbose =FALSE)

#Or you may want to generate 2D matrix 
#(mean of posterior) using the same prior
#estimates as generated before:
bayNorm_2D<-bayNorm_sup(
    Data=rse,
    PRIORS=bayNorm_3D$PRIORS,
    input_params=bayNorm_3D$input_params,
    mode_version=FALSE,
    mean_version = TRUE,
    verbose =FALSE)

## ----DE_fun, echo=TRUE,eval=FALSE---------------------------------------------
#  library(MAST)
#  library(reshape2)
#  ########SCnorm_runMAST3##########
#  SCnorm_runMAST3 <- function(Data, NumCells) {
#      if(length(dim(Data))==2){
#          resultss<-SCnorm_runMAST(Data, NumCells)
#      } else if(length(dim(Data))==3){
#  
#          library(foreach)
#          resultss<- foreach(sampleind=1:dim(Data)[3],.combine=cbind)%do%{
#              print(sampleind)
#  
#              qq<-SCnorm_runMAST(Data[,,sampleind], NumCells)
#              return(qq$adjpval)
#          }
#      }
#  
#      return(resultss)
#  }
#  
#  SCnorm_runMAST <- function(Data, NumCells) {
#  
#    NA_ind<-which(rowSums(Data)==0)
#    Data = as.matrix(log2(Data+1))
#  
#  
#    G1<-Data[,seq(1,NumCells[1])]
#    G2<-Data[,-seq(1,NumCells[1])]
#  
#    qq_temp<- rowMeans(G2)-rowMeans(G1)
#    qq_temp[NA_ind]=NA
#  
#    numGenes = dim(Data)[1]
#    datalong = melt(Data)
#    Cond = c(rep("c1", NumCells[1]*numGenes), rep("c2", NumCells[2]*numGenes))
#    dataL = cbind(datalong, Cond)
#    colnames(dataL) = c("gene","cell","value","Cond")
#    dataL$gene = factor(dataL$gene)
#    dataL$cell = factor(dataL$cell)
#    vdata = FromFlatDF(dataframe = dataL, idvars = "cell", primerid = "gene", measurement = "value",  id = numeric(0), cellvars = "Cond", featurevars = NULL,  phenovars = NULL)
#  
#    zlm.output = zlm(~ Cond, vdata, method='glm', ebayes=TRUE)
#  
#    zlm.lr = lrTest(zlm.output, 'Cond')
#  
#    gpval = zlm.lr[,,'Pr(>Chisq)']
#    adjpval = p.adjust(gpval[,1],"BH") ## Use only pvalues from the continuous part.
#    adjpval = adjpval[rownames(Data)]
#    return(list(adjpval=adjpval, logFC_re=qq_temp))
#  }
#  
#  
#  
#  

## ----DE_run, echo=TRUE,eval=FALSE---------------------------------------------
#  #Now, detect DE genes between two groups of cells with 15 cells in each group respectively
#  
#  #For 3D array
#  DE_out<-SCnorm_runMAST3(Data=bayNorm_3D$Bay_out, NumCells=c(15,15))
#  med_pvals<-apply(DE_out,1,median)
#  #DE genes are called with threshold 0.05:
#  DE_genes<-names(which(med_pvals<0.05))
#  
#  #For 2D array
#  DE_out<-SCnorm_runMAST3(Data=bayNorm_2D$Bay_out, NumCells=c(15,15))
#  DE_genes<-names(which(DE_out$adjpval<0.05))

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