QuasiDE: A new quasi-likelihood method to analyze the RNA-Seq data

In chushugu/QuasiDE: A new quasi-likelihood method to analze RNA-Seq data

Description

A new quasi-likelihood method to analyze the RNA-Seq data

Usage

QuasiDE(data.norm, design1, design2 = NULL, splmodel)

Arguments

data.norm	Raw count data after normalization
design1	Design matrix for the full model
design2	Design matrix for the redcued model
splmodel	The average variance function estimated from the data (spline)

Value

Returns raw p-values

Examples

## Load example data

library(airway)
library(edgeR)
data(airway)
airraw <- assays(airway)$count
air <- airraw[,c(1,3,5,7,2,4,6,8)]
n1 <- 4
n2 <- 4
n <- 8
trt<-c(rep(1,n1),rep(2,n2))
## Filtering
IQR <- apply(air,1,IQR)
airf <- air[IQR != 0 & rowMeans(air)>1,]

## Normalization

libsize = apply(airf[,1:n],2,sum)
nrmfactor <-calcNormFactors(airf[,1:n],method = "TMM")*libsize
air.norm <- as.data.frame(t(t(airf[,1:n])/nrmfactor))*mean(nrmfactor)

## Estimate average variance function

air.norm$mean <- apply(air.norm[,1:n],1,mean)
air.norm$mean1 <- apply(air.norm[,1:n1],1,mean)
air.norm$mean2 <- apply(air.norm[,(n1+1):n],1,mean)
air.norm$var <- apply(air.norm[,1:n],1,var)
air.norm$var1 <- apply(air.norm[,1:n1],1,var)
air.norm$var2 <- apply(air.norm[,(n1+1):n],1,var)
meanest <- c(air.norm$mean1,air.norm$mean2)
varest <- c(air.norm$var1,air.norm$var2)
modnorm <- smooth.spline(meanest[varest != 0],log(varest[varest != 0]))

## Analyze the normalized data

## Use a subset for testing
air.norm <- air.norm[1:3000,]  
dsgn <- model.matrix(~as.factor(trt))
air.norm$QuasiDE.rawp <- QuasiDE(air.norm[,1:n],design1=dsgn,splmodel=modnorm)
air.norm$QuasiDE.adjp <- p.adjust(air.norm$QuasiDE.rawp,method="BH")
sum(air.norm$QuasiDE.adjp<0.05)

chushugu/QuasiDE documentation built on May 18, 2019, 8:11 p.m.