Calculate the test statistics of the AMAP tests

Share:

Description

Calculate the test statistics of the AMAP tests

Usage

1
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
######## 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)