Supported by Dr. Osamu Ogasawara and  providing  . |
|
Last data update: 2014.03.03 |
Bayesian analysis of multivariate counts data in DNA metabarcoding and ecology.DescriptionThis package enables MCMC-based generalized linear mixed model analysis of multivariate counts data, such as common in DNA metabarcoding and community ecology. The results are summarized and plotted using ggplot2 functions. Details
At the moment, the package handles experimental design with a single multilevel fixed factor or two fully crossed multilevel fixed factors. Any number of scalar (OTU-specific) random factors (i.e. blocking factors) are allowed. By default, it is assumed that variation in total counts per sample is not biologically relevant (refects sequencing or survey effort). See help for the core function mcmc.otu() for more details. Author(s)Mikhail V. Matz <matz@utexas.edu> ReferencesElizabeth A. Green, Sarah W. Davies, Mikhail V. Matz, Monica Medina Quantifying cryptic Symbiodinium diversity within Orbicella faveolata and Orbicella franksi at the Flower Garden Banks, Gulf of Mexico. PeerJ 2014 2:e386 https://peerj.com/articles/386/ Examples# Symbiodinium sp diversity in two coral species at two reefs (banks) data(green.data) # removing outliers goods=purgeOutliers( data=green.data, count.columns=c(4:length(green.data[1,])), zero.cut=0.25 # remove this line for real analysis ) # stacking the data table gs=otuStack( data=goods, count.columns=c(4:length(goods[1,])), condition.columns=c(1:3) ) # fitting the model mm=mcmc.otu( fixed="bank+species+bank:species", data=gs, nitt=3000,burnin=2000 # remove this line for real analysis! ) # selecting the OTUs that were modeled reliably acpass=otuByAutocorr(mm,gs) # calculating effect sizes and p-values: ss=OTUsummary(mm,gs,summ.plot=FALSE) # correcting for mutliple comparisons (FDR) ss=padjustOTU(ss) # getting significatly changing OTUs (FDR<0.05) sigs=signifOTU(ss) # plotting them ss2=OTUsummary(mm,gs,otus=sigs) # bar-whiskers graph of relative changes: # ssr=OTUsummary(mm,gs,otus=signifOTU(ss),relative=TRUE) # displaying effect sizes and p-values for significant OTUs ss$otuWise[sigs] Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 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(MCMC.OTU)
Loading required package: MCMCglmm
Loading required package: Matrix
Loading required package: coda
Loading required package: ape
Loading required package: ggplot2
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MCMC.OTU/MCMC.OTU-package.Rd_%03d_medium.png", width=480, height=480)
> ### Name: MCMC.OTU-package
> ### Title: Bayesian analysis of multivariate counts data in DNA
> ### metabarcoding and ecology.
> ### Aliases: MCMC.OTU-package MCMC.OTU
>
> ### ** Examples
>
>
> # Symbiodinium sp diversity in two coral species at two reefs (banks)
> data(green.data)
>
> # removing outliers
> goods=purgeOutliers(
+ data=green.data,
+ count.columns=c(4:length(green.data[1,])),
+ zero.cut=0.25 # remove this line for real analysis
+ )
[1] "samples with counts below z-score -2.5 :"
[1] "EFAV153"
[1] "zscores:"
EFAV153
-3.884085
[1] "OTUs passing frequency cutoff 0.001 : 10"
[1] "OTUs with counts in 0.25 of samples:"
FALSE TRUE
3 7
>
> # stacking the data table
> gs=otuStack(
+ data=goods,
+ count.columns=c(4:length(goods[1,])),
+ condition.columns=c(1:3)
+ )
>
> # fitting the model
> mm=mcmc.otu(
+ fixed="bank+species+bank:species",
+ data=gs,
+ nitt=3000,burnin=2000 # remove this line for real analysis!
+ )
$PRIOR
$PRIOR$R
$PRIOR$R$V
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0
[5,] 0 0 0 0 1 0 0 0
[6,] 0 0 0 0 0 1 0 0
[7,] 0 0 0 0 0 0 1 0
[8,] 0 0 0 0 0 0 0 1
$PRIOR$R$nu
[1] 7.02
$PRIOR$G
$PRIOR$G$G1
$PRIOR$G$G1$V
[1] 1
$PRIOR$G$G1$nu
[1] 0
$FIXED
[1] "count~otu+bank+species+bank:species+otu:bank+otu:species+otu:bank:species"
$RANDOM
[1] "~sample"
MCMC iteration = 0
Acceptance ratio for liability set 1 = 0.000158
Acceptance ratio for liability set 2 = 0.000404
Acceptance ratio for liability set 3 = 0.000175
Acceptance ratio for liability set 4 = 0.000333
Acceptance ratio for liability set 5 = 0.000263
Acceptance ratio for liability set 6 = 0.000474
Acceptance ratio for liability set 7 = 0.000456
Acceptance ratio for liability set 8 = 0.000579
MCMC iteration = 1000
Acceptance ratio for liability set 1 = 0.211035
Acceptance ratio for liability set 2 = 0.422912
Acceptance ratio for liability set 3 = 0.220947
Acceptance ratio for liability set 4 = 0.435035
Acceptance ratio for liability set 5 = 0.342246
Acceptance ratio for liability set 6 = 0.422088
Acceptance ratio for liability set 7 = 0.367228
Acceptance ratio for liability set 8 = 0.509965
MCMC iteration = 2000
Acceptance ratio for liability set 1 = 0.294263
Acceptance ratio for liability set 2 = 0.424035
Acceptance ratio for liability set 3 = 0.300737
Acceptance ratio for liability set 4 = 0.453596
Acceptance ratio for liability set 5 = 0.379333
Acceptance ratio for liability set 6 = 0.419386
Acceptance ratio for liability set 7 = 0.375825
Acceptance ratio for liability set 8 = 0.521140
MCMC iteration = 3000
Acceptance ratio for liability set 1 = 0.297035
Acceptance ratio for liability set 2 = 0.447509
Acceptance ratio for liability set 3 = 0.318070
Acceptance ratio for liability set 4 = 0.481035
Acceptance ratio for liability set 5 = 0.359351
Acceptance ratio for liability set 6 = 0.411053
Acceptance ratio for liability set 7 = 0.389842
Acceptance ratio for liability set 8 = 0.477263
>
> # selecting the OTUs that were modeled reliably
> acpass=otuByAutocorr(mm,gs)
>
> # calculating effect sizes and p-values:
> ss=OTUsummary(mm,gs,summ.plot=FALSE)
>
> # correcting for mutliple comparisons (FDR)
> ss=padjustOTU(ss)
>
> # getting significatly changing OTUs (FDR<0.05)
> sigs=signifOTU(ss)
>
> # plotting them
> ss2=OTUsummary(mm,gs,otus=sigs)
>
> # bar-whiskers graph of relative changes:
> # ssr=OTUsummary(mm,gs,otus=signifOTU(ss),relative=TRUE)
>
> # displaying effect sizes and p-values for significant OTUs
> ss$otuWise[sigs]
$H35JRAZ01ATSK2
difference
pvalue bankeast:speciesfaveolata bankeast:speciesfranksi
bankeast:speciesfaveolata NA 0.082405883
bankeast:speciesfranksi 0.86400742 NA
bankwest:speciesfaveolata 0.02124985 0.003285617
bankwest:speciesfranksi 0.07235189 0.007556925
difference
pvalue bankwest:speciesfaveolata bankwest:speciesfranksi
bankeast:speciesfaveolata -0.5281631 -0.43595340
bankeast:speciesfranksi -0.6105690 -0.51835928
bankwest:speciesfaveolata NA 0.09220974
bankwest:speciesfranksi 0.8322792 NA
$H35JRAZ03DED1S
difference
pvalue bankeast:speciesfaveolata bankeast:speciesfranksi
bankeast:speciesfaveolata NA -0.31752780
bankeast:speciesfranksi 0.394117505 NA
bankwest:speciesfaveolata 0.003285617 0.02711453
bankwest:speciesfranksi 0.050959699 0.44127550
difference
pvalue bankwest:speciesfaveolata bankwest:speciesfranksi
bankeast:speciesfaveolata -1.2171363 -0.6525417
bankeast:speciesfranksi -0.8996085 -0.3350139
bankwest:speciesfaveolata NA 0.5645946
bankwest:speciesfranksi 0.1760919 NA
$H35JRAZ03DJU1V
difference
pvalue bankeast:speciesfaveolata bankeast:speciesfranksi
bankeast:speciesfaveolata NA -0.681777620
bankeast:speciesfranksi 0.42217817 NA
bankwest:speciesfaveolata 0.02544376 0.001727597
bankwest:speciesfranksi 0.02544376 0.001727597
difference
pvalue bankwest:speciesfaveolata bankwest:speciesfranksi
bankeast:speciesfaveolata 1.190879 1.11141242
bankeast:speciesfranksi 1.872657 1.79319004
bankwest:speciesfaveolata NA -0.07946675
bankwest:speciesfranksi 0.886129 NA
>
>
>
>
>
>
> dev.off()
null device
1
>
|
Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and