This page was generated on 2020-04-15 12:04:10 -0400 (Wed, 15 Apr 2020).
| EBarrays 2.50.0 Ming Yuan
 
 
| Snapshot Date: 2020-04-14 16:46:13 -0400 (Tue, 14 Apr 2020) |  | URL: https://git.bioconductor.org/packages/EBarrays |  | Branch: RELEASE_3_10 |  | Last Commit: a5415a9 |  | Last Changed Date: 2019-10-29 13:07:27 -0400 (Tue, 29 Oct 2019) |  | malbec1 | Linux (Ubuntu 18.04.4 LTS) / x86_64 | OK | OK | [ OK ] |  |  | 
| tokay1 | Windows Server 2012 R2 Standard / x64 | OK | OK | OK | OK |  | 
| merida1 | OS X 10.11.6 El Capitan / x86_64 | OK | OK | OK | OK |  | 
R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(EBarrays)
Loading required package: Biobase
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':
    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':
    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colnames, dirname, do.call,
    duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
    lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
    tapply, union, unique, unsplit, which, which.max, which.min
Welcome to Bioconductor
    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: lattice
> demo(ebarrays)
	demo(ebarrays)
	---- ~~~~~~~~
> library(EBarrays)
> ## EM algorithm 
> ## Lognormal-Normal Demo
> 
> ## mu10,sigma2,tau are parameters in the LNNB model; pde is the
> ## proportion of differentially expressed genes; n is the
> ## total number of genes; nr1 and nr2 are the number of replicate
> ## arrays in each group.
> 
> lnnb.sim <- function(mu10, sigmasq, tausq, pde, n, nr1, nr2)
+ {
+     de <- sample(c(TRUE, FALSE), size = n, replace = TRUE, prob = c(pde, 1 - pde))
+     x <- matrix(NA, n, nr1)
+     y <- matrix(NA, n, nr2)
+     mu1 <- rnorm(n, mu10, sqrt(tausq))
+     mu2.de <- rnorm(n, mu10, sqrt(tausq))
+     mu2 <- mu1
+     mu2[de] <- mu2.de[de]
+     for(j in 1:nr1) {
+         x[, j] <- rnorm(n, mu1, sqrt(sigmasq))
+     }
+     for(j in 1:nr2) {
+         y[, j] <- rnorm(n, mu2, sqrt(sigmasq))
+     }
+     outmat <- exp(cbind(x, y))
+     list(mu1 = mu1, mu2 = mu2, outmat = outmat, de = de)
+ }
> ## simulating data with
> ##  mu_0 = 2.33, sigma^2 = 0.1, tau^2 = 2
> ##  P(DE) = 0.2
> 
> sim.data1 <- lnnb.sim(2.33, 0.1, 2, 0.2, 2000, nr1 = 3, nr2 = 3)
> de.true1 <- sim.data1$de ## true indicators of differential expression
> sim.data2 <- lnnb.sim(1.33, 0.01, 2, 0.2, 2000, nr1 = 3, nr2 = 3)
> de.true2 <- sim.data2$de ## true indicators of differential expression
> testdata <- rbind(sim.data1$outmat,sim.data2$outmat)
> hypotheses <- ebPatterns(c("1 1 1 1 1 1", "1 1 1 2 2 2")) 
> em.out <- emfit(testdata, family = "LNN", hypotheses,
+                 cluster = 1:5,
+                 type = 2,
+                 verbose = TRUE,
+                 num.iter = 10)
 Checking for negative entries...
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  0.43 seconds user time
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  1.20 seconds user time
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  1.96 seconds user time
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  2.34 seconds user time
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  2.88 seconds user time
> em.out
 EB model fit 
	 Family: LNN ( Lognormal-Normal )
 Model parameter estimates:
              mu_0     sigma.2  tao_0.2
Cluster 1 2.364993 0.098960402 1.981286
Cluster 2 1.279346 0.009695755 2.038273
 Estimated mixing proportions:
          Pattern.1  Pattern.2
Cluster 1 0.4031537 0.10062453
Cluster 2 0.3972889 0.09893293
> post.out <- postprob(em.out, testdata)
> table(post.out$pattern[, 2] > .5, c(de.true1,de.true2))
       
        FALSE TRUE
  FALSE  3191  134
  TRUE     31  644
> table((post.out$cluster[, 2] > .5)+1, c(rep("Cluster 1",2000),rep("Cluster 2",2000)))
   
    Cluster 1 Cluster 2
  1      1874        78
  2       126      1922
> plotMarginal(em.out,testdata)
> par(ask=TRUE)
> plotCluster(em.out,testdata)
> par(ask=FALSE)
> lnnmv.em.out <- emfit(testdata, family = "LNNMV", hypotheses, groupid=c(1,1,1,2,2,2),
+                 verbose = TRUE,
+                 num.iter = 10,
+                 p.init = c(0.95, 0.05))
 Checking for negative entries...
 Generating summary statistics for patterns. 
 This may take a few seconds...
 Starting EM iterations (total 10 ).
 This may take a while
	 Starting iteration 1 ...
	 Starting iteration 2 ...
	 Starting iteration 3 ...
	 Starting iteration 4 ...
	 Starting iteration 5 ...
	 Starting iteration 6 ...
	 Starting iteration 7 ...
	 Starting iteration 8 ...
	 Starting iteration 9 ...
	 Starting iteration 10 ...
 Fit used  0.74 seconds user time
> lnnmv.em.out
 EB model fit 
	 Family: LNNMV ( Lognormal-Normal with modified variances )
 Model parameter estimates:
      mu_0  tao_0.2
1 1.833301 2.297505
 Estimated mixing proportions:
       Pattern.1 Pattern.2
p.temp 0.7877694 0.2122306
> post.out <- postprob(lnnmv.em.out, testdata, groupid=c(1,1,1,2,2,2))
> table(post.out$pattern[, 2] > .5, c(de.true1,de.true2))
       
        FALSE TRUE
  FALSE  3148  125
  TRUE     74  653
There were 50 or more warnings (use warnings() to see the first 50)
> 
> 
> 
> proc.time()
   user  system elapsed 
 11.492   0.108  11.819