Calculate the test statistics of the AMAP tests

1 |

`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 |

`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 |

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

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)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.