Description Usage Arguments Value References Examples
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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.