test.AMAP: Calculate the test statistics of the AMAP tests

In AMAP.Seq: Compare Gene Expressions from 2-Treatment RNA-Seq Experiments

Description

Calculate the test statistics of the AMAP tests

Usage

test.AMAP(data, MGN, del.lim = NULL, FC = NULL, print.steps = FALSE, Integration = "MC",nMC=NULL)

Arguments

data	RNA-seq data standardized by function RNASeq.Data()
MGN	The joint distribution, pi(lambda,delta), in form of Mixture Gamam-Normal
del.lim	An interval, for example del.lim=c(-1,1), that is the null space for delta
FC	A number >=1 so that the test detects genes with fold-changes greater than FC. If to detect DE genes, FC=1.
print.steps	Print the process when calculating the test statistics
Integration	Value can be "grid" or "MC". If Integration="grid", then the integration is done by dividing the 2-D space into grids. If Integration="MC", then the integration is done by Monte Carlo sampling.
nMC	number of data points randomly drawn from MGN distribution by Monte Carlo simulation,the default is 50000

Value

stat	test statistics of the AMAP tests, in logarithm scale
prob	posterior probability of the null hypothesis, equal to exp(stat)
fdr	estiamted FDR level if the cut-off is chosen at the gene

References

Yaqing Si and Peng Liu (2012), An Optimal Test with Maximum Average Power While Controlling FDR with Application to RNA-seq Data

Examples

######## Please read the help instribution above and the manuscript to 
######## CHOOSE PROPER PARAMETERS LIKE nK, iter.max, nK0, FC and nMC for best use of the function 
set.seed(100)
data("SimuHapMap")  # a matrix 'counts' storing simulated data with 10000 genes, two treatments, of which each has 5 replicates
head(cbind(counts,del.true))
counts=counts[1:200,] ### use data for only 200 genes to save time for testing example
                      ### the computation usually requires tens of minutes for 10000 genes
group=rep(1:2,each=5)

### standardize the RNA-seq data

size=Norm.GMedian(counts)   ## normalizing factor using Geometric Median
mydata=RNASeq.Data(counts=counts,size=size,group=group,model="nbinom")

### test DE genes

decom.est=MGN.EM(mydata,nK=3,iter.max=3,nK0=3,nMC=100)
s1=test.AMAP(mydata,MGN=decom.est$MGN,FC=1.0,nMC=100)
head(s1)

### test for FC>1.1

decom.est=MGN.EM(mydata,nK=3,iter.max=3,nK0=0,nMC=100)
s2=test.AMAP(mydata,MGN=decom.est$MGN,FC=1.1,nMC=100)
head(s2)

AMAP.Seq documentation built on May 2, 2019, 6:45 a.m.