@@ -4,40 +4,53 @@ Version: 1.2.1

 Date: 2013-06-11

 Author: Joseph Nathaniel Paulson, Mihai Pop,  Hector Corrada Bravo

 Maintainer: Joseph N. Paulson <jpaulson@umiacs.umd.edu>

-Description: metagenomeSeq is designed to determine features (be it Operational Taxanomic Unit (OTU), species, etc.) that are differentially abundant between two or more groups of multiple samples. metagenomeSeq is designed to address the effects of both normalization and under-sampling of microbial communities on disease association detection and the testing of feature correlations.

+Description: metagenomeSeq is designed to determine features (be it Operational

+    Taxanomic Unit (OTU), species, etc.) that are differentially abundant

+    between two or more groups of multiple samples. metagenomeSeq is designed

+    to address the effects of both normalization and under-sampling of

+    microbial communities on disease association detection and the testing of

+    feature correlations.

 License: Artistic-2.0

-Depends: R(>= 3.0), Biobase, limma, matrixStats, methods,

-        RColorBrewer, gplots

-Suggests: annotate

+Depends:

+    R(>= 3.0),

+    Biobase,

+    limma,

+    matrixStats,

+    methods,

+    RColorBrewer,

+    gplots

+Suggests:

+    annotate

 biocViews: Bioinformatics, DifferentialExpression, Metagenomics,

-        Visualization

-Collate: 'zigControl.R' 

-	'cumNorm.R' 

-	'plotOTU.R' 

-	'fitZig.R'

-        'doCountMStep.R' 

-	'doZeroMStep.R' 

-	'doEStep.R' 

-	'getZ.R' 

-	'getPi.R'

-        'getCountDensity.R' 

-	'getNegativeLogLikelihoods.R'

-        'isItStillActive.R' 

-	'getEpsilon.R' 

-	'load_meta.R'

-        'load_phenoData.R' 

-	'exportMat.R' 

-	'exportStats.R'

-        'cumNormStat.R'

-	'plotGenus.R'

-	'aggregateM.R' 

-	'cumNormMat.R'

-        'load_metaQ.R' 

-	'allClasses.R' 

-	'MRtable.R' 

-	'MRcoefs.R'

-        'plotMRheatmap.R'

-	'plotCorr.R'

-	'MRfisher.R'

-	'MRfulltable.R'

+    Visualization

+Collate:

+    'zigControl.R'

+    'cumNorm.R'

+    'plotOTU.R'

+    'fitZig.R'

+    'doCountMStep.R'

+    'doZeroMStep.R'

+    'doEStep.R'

+    'getZ.R'

+    'getPi.R'

+    'getCountDensity.R'

+    'getNegativeLogLikelihoods.R'

+    'isItStillActive.R'

+    'getEpsilon.R'

+    'load_meta.R'

+    'load_phenoData.R'

+    'exportMat.R'

+    'exportStats.R'

+    'cumNormStat.R'

+    'plotGenus.R'

+    'aggregateM.R'

+    'cumNormMat.R'

+    'load_metaQ.R'

+    'allClasses.R'

+    'MRtable.R'

+    'MRcoefs.R'

+    'plotMRheatmap.R'

+    'plotCorr.R'

+    'MRfisher.R'

+    'MRfulltable.R'

 URL: http://cbcb.umd.edu/software/metagenomeSeq

@@ -1,3 +1,46 @@

+#' Table of top-ranked microbial marker gene from linear model fit

+#' 

+#' Extract a table of the top-ranked features from a linear model fit. This

+#' function will be updated soon to provide better flexibility similar to

+#' limma's topTable.

+#' 

+#' 

+#' @param obj A list containing the linear model fit produced by lmFit through

+#' fitZig.

+#' @param by Column number or column name specifying which coefficient or

+#' contrast of the linear model is of interest.

+#' @param coef Column number(s) or column name(s) specifying which coefficient

+#' or contrast of the linear model to display.

+#' @param number The number of bacterial features to pick out.

+#' @param taxa Taxa list.

+#' @param uniqueNames Number the various taxa.

+#' @param adjust.method Method to adjust p-values by. Default is "FDR". Options

+#' include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

+#' "none". See \code{\link{p.adjust}} for more details.

+#' @param group One of three choices, 0,1,2,3. 0: the sort is ordered by a

+#' decreasing absolute value coefficient fit. 1: the sort is ordered by the raw

+#' coefficient fit in decreasing order. 2: the sort is ordered by the raw

+#' coefficient fit in increasing order. 3: the sort is ordered by the p-value

+#' of the coefficient fit in increasing order.

+#' @param eff Restrict samples to have at least eff quantile effective samples.

+#' @param output Name of output file, including location, to save the table.

+#' @return Table of the top-ranked features determined by the linear fit's

+#' coefficient.

+#' @seealso \code{\link{fitZig}} \code{\link{MRtable}}

+#' @examples

+#' 

+#' data(lungData)

+#' k = grep("Extraction.Control",pData(lungData)$SampleType)

+#' lungTrim = lungData[,-k]

+#' k = which(rowSums(MRcounts(lungTrim)>0)<10)

+#' lungTrim = lungTrim[-k,]

+#' cumNorm(lungTrim)

+#' smokingStatus = pData(lungTrim)$SmokingStatus

+#' mod = model.matrix(~smokingStatus)

+#' settings = zigControl(maxit=1,verbose=FALSE)

+#' fit = fitZig(obj = lungTrim,mod=mod,control=settings)

+#' head(MRcoefs(fit))

+#' 

 MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FALSE,adjust.method="fdr",group=0,eff=0,output=NULL){

     tb = obj$fit$coefficients

     tx = as.character(taxa);

@@ -43,4 +86,4 @@ MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FALSE,a

     } else{

         return(as.data.frame(mat))

     }

-}

\ No newline at end of file

+}

@@ -1,3 +1,25 @@

+#' Wrapper to run fisher's test on presence/absence of a feature.

+#' 

+#' This function returns a data frame of p-values, odds ratios, lower and upper

+#' confidence limits for every row of a matrix.

+#' 

+#' 

+#' @param obj A MRexperiment object with a count matrix, or a simple count

+#' matrix.

+#' @param cl Group comparison

+#' @param mat logical indicating whether obj is a MRexperiment object or

+#' matrix. Default is a MRexperiment object.

+#' @return NA

+#' @seealso \code{\link{cumNorm}} \code{\link{fitZig}}

+#' @examples

+#' 

+#' data(lungData)

+#' k = grep("Extraction.Control",pData(lungData)$SampleType)

+#' lungTrim = lungData[,-k]

+#' lungTrim = lungTrim[-which(rowSums(MRcounts(lungTrim)>0)<20),]

+#' res = MRfisher(lungTrim,pData(lungTrim)$SmokingStatus);

+#' head(res)

+#' 

 MRfisher<-function(obj,cl,mat=FALSE){

     if(mat==FALSE){

         x = MRcounts(obj)>0;

@@ -1,3 +1,50 @@

+#' Table of top microbial marker gene from linear model fit including sequence

+#' information

+#' 

+#' Extract a table of the top-ranked features from a linear model fit. This

+#' function will be updated soon to provide better flexibility similar to

+#' limma's topTable. This function differs from \code{link{MRcoefs}} in that it

+#' provides other information about the presence or absence of features to help

+#' ensure significant features called are moderately present.

+#' 

+#' 

+#' @param obj A list containing the linear model fit produced by lmFit through

+#' fitZig.

+#' @param by Column number or column name specifying which coefficient or

+#' contrast of the linear model is of interest.

+#' @param coef Column number(s) or column name(s) specifying which coefficient

+#' or contrast of the linear model to display.

+#' @param number The number of bacterial features to pick out.

+#' @param taxa Taxa list.

+#' @param uniqueNames Number the various taxa.

+#' @param adjust.method Method to adjust p-values by. Default is "FDR". Options

+#' include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

+#' "none". See \code{\link{p.adjust}} for more details.

+#' @param group One of three choices, 0,1,2. 0: the sort is ordered by a

+#' decreasing absolute value coefficient fit. 1: the sort is ordered by the raw

+#' coefficient fit in decreasing order. 2: the sort is ordered by the raw

+#' coefficient fit in increasing order. 3: the sort is ordered by the p-value

+#' of the coefficient fit in increasing order.

+#' @param eff Restrict samples to have at least eff quantile effective samples.

+#' @param output Name of output file, including location, to save the table.

+#' @return Table of the top-ranked features determined by the linear fit's

+#' coefficient.

+#' @seealso \code{\link{fitZig}} \code{\link{MRcoefs}} \code{\link{MRtable}}

+#' \code{\link{MRfisher}}

+#' @examples

+#' 

+#' data(lungData)

+#' k = grep("Extraction.Control",pData(lungData)$SampleType)

+#' lungTrim = lungData[,-k]

+#' k = which(rowSums(MRcounts(lungTrim)>0)<10)

+#' lungTrim = lungTrim[-k,]

+#' cumNorm(lungTrim)

+#' smokingStatus = pData(lungTrim)$SmokingStatus

+#' mod = model.matrix(~smokingStatus)

+#' settings = zigControl(maxit=1,verbose=FALSE)

+#' fit = fitZig(obj = lungTrim,mod=mod,control=settings)

+#' head(MRfulltable(fit))

+#' 

 MRfulltable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FALSE,adjust.method="fdr",group=0,eff=0,output=NULL){

     tb = obj$fit$coefficients

@@ -62,4 +109,4 @@ MRfulltable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FAL

     } else{

         return(as.data.frame(mat))

     }

-}

\ No newline at end of file

+}

@@ -1,3 +1,48 @@

+#' Table of top microbial marker gene from linear model fit including sequence

+#' information

+#' 

+#' Extract a table of the top-ranked features from a linear model fit. This

+#' function will be updated soon to provide better flexibility similar to

+#' limma's topTable. This function differs from \code{link{MRcoefs}} in that it

+#' provides other information about the presence or absence of features to help

+#' ensure significant features called are moderately present.

+#' 

+#' 

+#' @param obj A list containing the linear model fit produced by lmFit through

+#' fitZig.

+#' @param by Column number or column name specifying which coefficient or

+#' contrast of the linear model is of interest.

+#' @param coef Column number(s) or column name(s) specifying which coefficient

+#' or contrast of the linear model to display.

+#' @param number The number of bacterial features to pick out.

+#' @param taxa Taxa list.

+#' @param uniqueNames Number the various taxa.

+#' @param adjust.method Method to adjust p-values by. Default is "FDR". Options

+#' include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

+#' "none". See \code{\link{p.adjust}} for more details.

+#' @param group One of three choices, 0,1,2. 0: the sort is ordered by a

+#' decreasing absolute value coefficient fit. 1: the sort is ordered by the raw

+#' coefficient fit in decreasing order. 2: the sort is ordered by the raw

+#' coefficient fit in increasing order. 3: the sort is ordered by the p-value

+#' of the coefficient fit in increasing order.

+#' @param output Name of output file, including location, to save the table.

+#' @return Table of the top-ranked features determined by the linear fit's

+#' coefficient.

+#' @seealso \code{\link{fitZig}} \code{\link{MRcoefs}}

+#' @examples

+#' 

+#' data(lungData)

+#' k = grep("Extraction.Control",pData(lungData)$SampleType)

+#' lungTrim = lungData[,-k]

+#' k = which(rowSums(MRcounts(lungTrim)>0)<10)

+#' lungTrim = lungTrim[-k,]

+#' cumNorm(lungTrim)

+#' smokingStatus = pData(lungTrim)$SmokingStatus

+#' mod = model.matrix(~smokingStatus)

+#' settings = zigControl(maxit=1,verbose=FALSE)

+#' fit = fitZig(obj = lungTrim,mod=mod,control=settings)

+#' head(MRtable(fit))

+#' 

 MRtable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FALSE,adjust.method="fdr",group=0,output=NULL){

     tb = obj$fit$coefficients

     tx = as.character(taxa);

@@ -56,4 +101,4 @@ MRtable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,uniqueNames=FALSE,a

     } else{

         return(as.data.frame(mat))

     }

-}

\ No newline at end of file

+}

@@ -1,14 +1,17 @@

-#' Aggregates the counts to a particular classification.

-#'

-#' This function takes an eSet object of data at a particular level with feature information allowing

-#' for aggregation of counts to a particular level. This method assumes taxa begin at the highest level and continue to the current level.

-#'

-#' @param obj An eSet object of count data.

+#' Aggregates counts by a particular classification.

+#' 

+#' This function takes a MRexperiment object of data at a particular level with

+#' feature information allowing for aggregation of counts to a particular

+#' level. This method assumes taxa begin at the highest level and continue to

+#' the current level, reverse assumes taxa begin at the lowest level.

+#' 

+#' 

+#' @param obj A MRexperiment object.

 #' @param lvl The level to go up (numeric, 1,2,3).

+#' @param taxa A vector of taxa annotations with splits

 #' @param split The way character strings in taxa in the obj are split.

-#' @return Updated eSet object with counts aggregated to the various taxanomic levels.

-#'

-#' @name aggregateM

+#' @return Updated object with counts aggregated to the various taxanomic

+#' levels.

 aggregateM <-

 function(obj,taxa,lvl,split=";"){

@@ -11,6 +11,33 @@ setMethod("[", "MRexperiment", function (x, i, j, ..., drop = FALSE) {

         obj

 })

+

+

+#' Create a MRexperiment object

+#' 

+#' This function creates a MRexperiment object from a matrix or data frame of

+#' count data.

+#' 

+#' See \code{\link{MRexperiment-class}} and \code{eSet} (from the Biobase

+#' package) for the meaning of the various slots.

+#' 

+#' @param counts A matrix or data frame of count data. The count data is

+#' representative of the number of reads annotated for a feature (be it gene,

+#' OTU, species, etc). Rows should correspond to features and columns to

+#' samples.

+#' @param phenoData An AnnotatedDataFrame with pertinent sample information.

+#' @param featureData An AnnotatedDataFrame with pertinent feature information.

+#' @param libSize libSize, library size, is the total number of reads for a

+#' particular sample.

+#' @param normFactors normFactors, the normalization factors used in either the

+#' model or as scaling factors of sample counts for each particular sample.

+#' @return an object of class MRexperiment

+#' @author Joseph N Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#' cnts = matrix(abs(rnorm(1000)),nc=10)

+#' obj <- newMRexperiment(cnts)

+#' 

 newMRexperiment <- function(counts, phenoData=NULL, featureData=NULL,libSize=NULL, normFactors=NULL) {

     counts= as.matrix(counts)

@@ -19,9 +46,46 @@ newMRexperiment <- function(counts, phenoData=NULL, featureData=NULL,libSize=NUL

     if( is.null( phenoData ) )

       phenoData   <- annotatedDataFrameFrom(counts, byrow=FALSE)

     if( is.null( libSize ) )

+

+

+#' Access sample depth of coverage from MRexperiment object

+#' 

+#' The libSize vector represents the column (sample specific) sums of features,

+#' i.e. the total number of reads for a sample. It is used by

+#' \code{\link{fitZig}}.

+#' 

+#' 

+#' @name libSize

+#' @aliases libSize,MRexperiment-method libSize

+#' @docType methods

+#' @param obj a \code{MRexperiment} object.

+#' @author Joseph N. Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#'    data(lungData)

+#'    head(libSize(lungData))

+#' 

       libSize <- as.matrix(colSums(counts))

       rownames(libSize) = colnames(counts)

     if( is.null( normFactors ) ){

+

+

+#' Access the normalization factors in a MRexperiment object

+#' 

+#' Function to access the scaling factors, aka the normalization factors, of

+#' samples in a MRexperiment object.

+#' 

+#' 

+#' @name normFactors

+#' @aliases normFactors,MRexperiment-method normFactors

+#' @docType methods

+#' @param obj a \code{MRexperiment} object.

+#' @author Joseph N. Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#'    data(lungData)

+#'    head(normFactors(lungData))

+#' 

       normFactors <- as.matrix(rep( NA_real_, length(libSize) ))

       rownames(normFactors) = rownames(libSize)

     }

@@ -45,6 +109,26 @@ setValidity( "MRexperiment", function( object ) {

     TRUE

 } )

+

+

+#' Accessor for the counts slot of a MRexperiment object

+#' 

+#' The counts slot holds the raw count data representing (along the rows) the

+#' number of reads annotated for a particular feature and (along the columns)

+#' the sample.

+#' 

+#' 

+#' @name MRcounts

+#' @aliases MRcounts,MRexperiment-method MRcounts

+#' @docType methods

+#' @param cnts a \code{MRexperiment} object.

+#' @param norm logical indicating whether or not to return normalized counts.

+#' @author Joseph N. Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#'    data(lungData)

+#'    head(MRcounts(lungData))

+#' 

 MRcounts <- function( obj ,norm=FALSE) {

    stopifnot( is( obj, "MRexperiment" ) )

    if(!norm){

@@ -57,6 +141,23 @@ MRcounts <- function( obj ,norm=FALSE) {

    }

 }

+

+

+#' Access the posterior probabilities that results from analysis

+#' 

+#' Accessing the posterior probabilities following a run through

+#' \code{\link{fitZig}}

+#' 

+#' 

+#' @name posterior.probs

+#' @aliases posterior.probs,MRexperiment-method posterior.probs

+#' @docType methods

+#' @param obj a \code{MRexperiment} object.

+#' @author Joseph N. Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#' # see vignette

+#' 

 posterior.probs <- function( obj ) {

    stopifnot( is( obj, "MRexperiment" ) )

    assayData(obj)[["z"]]

@@ -74,7 +175,26 @@ libSize<-function(obj){

    ls

 }

+

+

+#' Access MRexperiment object experiment data

+#' 

+#' The expSummary vectors represent the column (sample specific) sums of

+#' features, i.e. the total number of reads for a sample, libSize and also the

+#' normalization factors, normFactor.

+#' 

+#' 

+#' @name expSummary

+#' @aliases expSummary,MRexperiment-method expSummary

+#' @docType methods

+#' @param obj a \code{MRexperiment} object.

+#' @author Joseph N. Paulson, jpaulson@@umiacs.umd.edu

+#' @examples

+#' 

+#' data(mouseData)

+#' expSummary(mouseData)

+#' 

 expSummary<-function(obj){

   stopifnot( is( obj, "MRexperiment" ) )

   pData(obj@expSummary$expSummary)

-}

\ No newline at end of file

+}

@@ -1,13 +1,19 @@

 #' Cumulative sum scaling factors.

-#'

-#' Calculates each column's quantile and calculates the sum up to and including that quantile.

-#'

-#' @param jobj An eSet object.

+#' 

+#' Calculates each column's quantile and calculates the sum up to and including

+#' that quantile.

+#' 

+#' 

+#' @param obj An MRexperiment object.

 #' @param p The pth quantile.

 #' @return Vector of the sum up to and including a sample's pth quantile

-#'

-#' @name cumNorm

-#' @seealso \code{\link{fitZig}}

+#' @seealso \code{\link{fitZig}} \code{\link{cumNormStat}}

+#' @examples

+#' 

+#' data(mouseData)

+#' cumNorm(mouseData)

+#' head(normFactors(mouseData))

+#' 

 cumNorm <-

 function(obj,p=cumNormStat(obj)){

 	x = MRcounts(obj)

@@ -1,13 +1,19 @@

 #' Cumulative sum scaling factors.

-#'

-#' Calculates each column's quantile and calculates the sum up to and including that quantile.

-#'

-#' @param jobj An eSet object.

+#' 

+#' Calculates each column's quantile and calculates the sum up to and including

+#' that quantile.

+#' 

+#' 

+#' @param obj A MRexperiment object.

 #' @param p The pth quantile.

-#' @return jobj An updated eSet object with normalized counts

-#'

-#' @name cumNormMat

+#' @return Returns a matrix normalized by scaling counts up to and including

+#' the pth quantile.

 #' @seealso \code{\link{fitZig}} \code{\link{cumNorm}}

+#' @examples

+#' 

+#' data(mouseData)

+#' head(cumNormMat(mouseData))

+#' 

 cumNormMat <-

 function(obj,p= cumNormStat(obj)){

 ####################################################################################

@@ -1,12 +1,22 @@

-#' Cumulative normalization statistic.

-#'

-#' @param obj An eSet object.

-#' @param pFlag Whether or not to plot the reference.

-#' @param rel Relative difference of rel percent.

-#' @return P-value for which to cumulative normalize.

-#'

-#' @name cumNormStat

+#' Cumulative sum scaling percentile selection

+#' 

+#' Calculates the percentile for which to sum counts up to and scale by.

+#' 

+#' 

+#' @param obj A list with count data

+#' @param pFlag Plot the median difference quantiles

+#' @param rel Cutoff for the relative difference from one median difference

+#' from the reference to the next

+#' @param qFlag Flag to either calculate the proper percentile using a

+#' step-wise or triangular approximation of the sample count distribution.

+#' @param ... Applicable if pFlag == TRUE. Extra plotting parameters.

+#' @return Percentile for which to scale data

 #' @seealso \code{\link{fitZig}} \code{\link{cumNorm}}

+#' @examples

+#' 

+#' data(mouseData)

+#' p = round(cumNormStat(mouseData,pFlag=FALSE),digits=2)

+#' 

 cumNormStat <-

 function(obj,pFlag = FALSE,rel=.1,qFlag = TRUE, ...){

@@ -42,4 +52,4 @@ function(obj,pFlag = FALSE,rel=.1,qFlag = TRUE, ...){

     x = which(abs(diff(diffr2))/diffr2[-1]>rel)[1] / length(diffr2)

     obj@expSummary$cumNormStat = x;

 	return(x)

-}

\ No newline at end of file

+}

@@ -1,22 +1,26 @@

 #' Compute the Maximization step calculation for features still active.

-#'

-#' Maximization step is solved by weighted least squares. The function also computes counts residuals.

-#'

-#' Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $\deta_{ij}$ = 1 if $y_{ij}$

-#' is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_{ij} = \pi_j(S_j) \cdot f_{0}(y_{ij}) 

-#' +(1-\pi_j (S_j))\cdot f_{count}(y_{ij};\mu_i,\sigma_i^2)$.

-#' The log-likelihood in this extended model is

-#' $(1−\delta_{ij}) \log f_{count}(y;\mu_i,\sigma_i^2 )+\delta_{ij} \log \pi_j(s_j)+(1−\delta_{ij})\log (1−\pi_j (sj))$.

-#' The responsibilities are defined as $z_{ij} = pr(\delta_{ij}=1 | data)$.

-#'

-#' @param z Matrix (m x n) of estimate responsibilities (probabilities that a count comes from a spike distribution at 0).

+#' 

+#' Maximization step is solved by weighted least squares.  The function also

+#' computes counts residuals.

+#' 

+#' Maximum-likelihood estimates are approximated using the EM algorithm where

+#' we treat mixture membership $delta_ij$ = 1 if $y_ij$ is generated from the

+#' zero point mass as latent indicator variables. The density is defined as

+#' $f_zig(y_ij = pi_j(S_j)*f_0(y_ij) +(1-pi_j (S_j)) *

+#' f_count(y_ij;mu_i,sigma_i^2)$. The log-likelihood in this extended model is

+#' $(1-delta_ij) log f_count(y;mu_i,sigma_i^2 )+delta_ij log

+#' pi_j(s_j)+(1-delta_ij)log (1-pi_j (s_j))$. The responsibilities are defined

+#' as $z_ij = pr(delta_ij=1 | data)$.

+#' 

+#' @param z Matrix (m x n) of estimate responsibilities (probabilities that a

+#' count comes from a spike distribution at 0).

 #' @param y Matrix (m x n) of count observations.

 #' @param mmCount Model matrix for the count distribution.

-#' @param stillActive Boolean vector of size M, indicating whether a feature converged or not.

+#' @param stillActive Boolean vector of size M, indicating whether a feature

+#' converged or not.

 #' @param fit2 Previous fit of the count model.

-#' @return Update matrix (m x n) of estimate responsibilities (probabilities that a count comes from a spike distribution at 0).

-#'

-#' @name doCountMStep

+#' @return Update matrix (m x n) of estimate responsibilities (probabilities

+#' that a count comes from a spike distribution at 0).

 #' @seealso \code{\link{fitZig}}

 doCountMStep <-

 function(z, y, mmCount, stillActive,fit2=NULL){

@@ -1,20 +1,22 @@

 #' Compute the Expectation step.

-#'

-#' Estimates the responsibilities $z_{ij} = \frac{\pi_j \cdot I_{0}(y_{ij}}{\pi_j \cdot I_{0}(y_{ij} + (1-\pi_j) \cdot f_{count}(y_{ij}}

-#'

-#' Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $\deta_{ij}$ = 1 if $y_{ij}$

-#' is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_{ij} = \pi_j(S_j) \cdot f_{0}(y_{ij}) 

-#' +(1-\pi_j (S_j))\cdot f_{count}(y_{ij};\mu_i,\sigma_i^2)$.

-#' The log-likelihood in this extended model is

-#' $(1−\delta_{ij}) \log f_{count}(y;\mu_i,\sigma_i^2 )+\delta_{ij} \log \pi_j(s_j)+(1−\delta_{ij})\log (1−\pi_j (sj))$.

-#' The responsibilities are defined as $z_{ij} = pr(\delta_{ij}=1 | data)$.

-#'

+#' 

+#' Estimates the responsibilities $z_ij = fracpi_j cdot I_0(y_ijpi_j cdot

+#' I_0(y_ij + (1-pi_j) cdot f_count(y_ij

+#' 

+#' Maximum-likelihood estimates are approximated using the EM algorithm where

+#' we treat mixture membership $delta_ij$ = 1 if $y_ij$ is generated from the

+#' zero point mass as latent indicator variables. The density is defined as

+#' $f_zig(y_ij = pi_j(S_j) cdot f_0(y_ij) +(1-pi_j (S_j))cdot

+#' f_count(y_ij;mu_i,sigma_i^2)$. The log-likelihood in this extended model is

+#' $(1-delta_ij) log f_count(y;mu_i,sigma_i^2 )+delta_ij log

+#' pi_j(s_j)+(1-delta_ij)log (1-pi_j (sj))$. The responsibilities are defined

+#' as $z_ij = pr(delta_ij=1 | data)$.

+#' 

 #' @param countResiduals Residuals from the count model.

 #' @param zeroResiduals Residuals from the zero model.

 #' @param zeroIndices Index (matrix m x n) of counts that are zero/non-zero.

-#' @return Updated matrix (m x n) of estimate responsibilities (probabilities that a count comes from a spike distribution at 0).

-#'

-#' @name doEStep

+#' @return Updated matrix (m x n) of estimate responsibilities (probabilities

+#' that a count comes from a spike distribution at 0).

 #' @seealso \code{\link{fitZig}}

 doEStep <-

 function(countResiduals,  zeroResiduals, zeroIndices)

@@ -1,22 +1,25 @@

 #' Compute the zero Maximization step.

-#'

-#' Performs Maximization step calculation for the mixture components. Uses least squares to fit the parameters of the mean of the logistic distribution.

-#' $$

-#' pi_j = \sum_i^M \frac{1}{M}z_{ij}

-#' $$

-#' Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $\deta_{ij}$ = 1 if $y_{ij}$

-#' is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_{ij} = \pi_j(S_j) \cdot f_{0}(y_{ij}) 

-#' +(1-\pi_j (S_j))\cdot f_{count}(y_{ij};\mu_i,\sigma_i^2)$.

-#' The log-likelihood in this extended model is

-#' $(1−\delta_{ij}) \log f_{count}(y;\mu_i,\sigma_i^2 )+\delta_{ij} \log \pi_j(s_j)+(1−\delta_{ij})\log (1−\pi_j (sj))$.

-#' The responsibilities are defined as $z_{ij} = pr(\delta_{ij}=1 | data)$.

-#'

-#' @param z Matrix (m x n) of estimate responsibilities (probabilities that a count comes from a spike distribution at 0).

+#' 

+#' Performs Maximization step calculation for the mixture components. Uses

+#' least squares to fit the parameters of the mean of the logistic

+#' distribution. $$ pi_j = sum_i^M frac1Mz_ij $$ Maximum-likelihood estimates

+#' are approximated using the EM algorithm where we treat mixture membership

+#' $delta_ij$ = 1 if $y_ij$ is generated from the zero point mass as latent

+#' indicator variables. The density is defined as $f_zig(y_ij = pi_j(S_j) cdot

+#' f_0(y_ij) +(1-pi_j (S_j))cdot f_count(y_ij;mu_i,sigma_i^2)$. The

+#' log-likelihood in this extended model is $(1-delta_ij) log

+#' f_count(y;mu_i,sigma_i^2 )+delta_ij log pi_j(s_j)+(1-delta_ij)log (1-pi_j

+#' (sj))$. The responsibilities are defined as $z_ij = pr(delta_ij=1 | data)$.

+#' 

+#' 

+#' @param z Matrix (m x n) of estimate responsibilities (probabilities that a

+#' count comes from a spike distribution at 0).

 #' @param zeroIndices Index (matrix m x n) of counts that are zero/non-zero.

-#' @param mmZero The zero model, the model matrix to account for the change in the number of OTUs observed as a linear effect of the depth of coverage.

-#' @return List of the zero fit (zero mean model) coefficients, variance - scale parameter (scalar), and normalized residuals of length sum(zeroIndices).

-#'

-#' @name doZeroMStep

+#' @param mmZero The zero model, the model matrix to account for the change in

+#' the number of OTUs observed as a linear effect of the depth of coverage.

+#' @return List of the zero fit (zero mean model) coefficients, variance -

+#' scale parameter (scalar), and normalized residuals of length

+#' sum(zeroIndices).

 #' @seealso \code{\link{fitZig}}

 doZeroMStep <-

 function(z, zeroIndices, mmZero)

@@ -1,18 +1,18 @@

 #' export the normalized eSet dataset as a matrix.

-#'

-#' This function allows the user to take the normalized dataset or counts and output

-#' the dataset to the user's workspace as a tab-delimited file, etc.

-#'

-#' @param jobj An eSet object with count data.

-#' @param output Output file name.

+#' 

+#' This function allows the user to take a dataset of counts and output the

+#' dataset to the user's workspace as a tab-delimited file, etc.

+#' 

+#' 

+#' @aliases exportMatrix exportMat

+#' @param mat A matrix of values (normalized, or otherwise)

+#' @param output Output file name

 #' @return NA

-#'

-#' @name export_mat

-#' @aliases exportMatrix

 #' @seealso \code{\link{cumNorm}}

 #' @examples

-#' export_mat(jobj,output="~/Desktop/normMatrix.tsv");

-

+#' 

+#' # see vignette

+#' 

 exportMat <-

 function(mat,output="~/Desktop/matrix.tsv"){

 	matrix = mat;

@@ -1,17 +1,19 @@

 #' Various statistics of the count data.

-#'

-#' A matrix of values for each sample. The matrix consists of sample ids, the sample scaling factor, quantile value, and the number of number of features.

-#'

-#' @param obj An eSet object with count data.

-#' @param p Quantile value to calculate the scaling factor and quantiles for the various samples.

+#' 

+#' A matrix of values for each sample. The matrix consists of sample ids, the

+#' sample scaling factor, quantile value, and the number of number of features.

+#' 

+#' 

+#' @param obj A MRexperiment object with count data.

+#' @param p Quantile value to calculate the scaling factor and quantiles for

+#' the various samples.

 #' @param output Output file name.

 #' @return None.

-#'

-#' @name export_stats

 #' @seealso \code{\link{cumNorm}} \code{\link{quantile}}

 #' @examples

-#' export_stats(obj,p=1,output="~/Desktop/obj-stats.tsv")

-

+#' 

+#' # see vignette

+#' 

 exportStats <-

 function(obj,p= cumNormStat(obj),output="~/Desktop/res.stats.tsv"){

@@ -1,27 +1,41 @@

 #' Computes the weighted fold-change estimates and t-statistics.

-#'

-#' Wrapper to actually run the Expectation-maximization algorithm and estimate $f_{count}$ fits.

-#' Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $\deta_{ij}$ = 1 if $y_{ij}$

-#' is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_{ij} = \pi_j(S_j) \cdot f_{0}(y_{ij}) 

-#' +(1-\pi_j (S_j))\cdot f_{count}(y_{ij};\mu_i,\sigma_i^2)$.

-#' The log-likelihood in this extended model is

-#' $(1−\delta_{ij}) \log f_{count}(y;\mu_i,\sigma_i^2 )+\delta_{ij} \log \pi_j(s_j)+(1−\delta_{ij})\log (1−\pi_j (sj))$.

-#' The responsibilities are defined as $z_{ij} = pr(\delta_{ij}=1 | data)$.

-#'

-#' @param obj An eSet object with count data.

+#' 

+#' Wrapper to actually run the Expectation-maximization algorithm and estimate

+#' $f_count$ fits.  Maximum-likelihood estimates are approximated using the EM

+#' algorithm where we treat mixture membership $delta_ij = 1$ if $y_ij$ is

+#' generated from the zero point mass as latent indicator variables. The

+#' density is defined as $f_zig(y_ij = pi_j(S_j)*f_0(y_ij) +(1-pi_j (S_j)) *

+#' f_count(y_ij; mu_i, sigma_i^2)$. The log-likelihood in this extended model

+#' is: $(1-delta_ij) log f_count(y;mu_i,sigma_i^2 )+delta_ij log

+#' pi_j(s_j)+(1-delta_ij) log (1-pi_j (s_j))$. The responsibilities are defined

+#' as $z_ij = pr(delta_ij=1 | data)$.

+#' 

+#' 

+#' @param obj A MRexperiment object with count data.

 #' @param mod The model for the count distribution.

-#' @param s95 A vector of size M of the scaling values to be included in the model.

-#' @param zeroMod The zero model, the model to account for the change in the number of OTUs observed as a linear effect of the depth of coverage.

-#' @param useS95offset Boolean, whether to include the default scaling parameters in the model or not.

-#' @param control The settings for fitZig. 

-#' @param s The raw total counts for the various samples.

-#' @return The fits, posterior probabilities, posterior probabilities used at time of convergence for each feature, ebayes (limma object) fit, among other data. 

-#'

-#' @name fitZig

+#' @param zeroMod The zero model, the model to account for the change in the

+#' number of OTUs observed as a linear effect of the depth of coverage.

+#' @param useS95offset Boolean, whether to include the default scaling

+#' parameters in the model or not.

+#' @param control The settings for fitZig.

+#' @return The fits, posterior probabilities, posterior probabilities used at

+#' time of convergence for each feature, ebayes (limma object) fit, among other

+#' data.

+#' @export

 #' @seealso \code{\link{cumNorm}} \code{\link{zigControl}}

 #' @examples

-#' model = model.matrix(~1+type+log2(s95/1000+1))

-#' res = fitZig(obj = obj,mod=mod,useS95offset=FALSE)

+#' 

+#' data(lungData)

+#' k = grep("Extraction.Control",pData(lungData)$SampleType)

+#' lungTrim = lungData[,-k]

+#' k = which(rowSums(MRcounts(lungTrim)>0)<30)

+#' cumNorm(lungTrim)

+#' lungTrim = lungTrim[-k,]

+#' smokingStatus = pData(lungTrim)$SmokingStatus

+#' mod = model.matrix(~smokingStatus)

+#' settings = zigControl(maxit=1,verbose=FALSE)

+#' fit = fitZig(obj = lungTrim,mod=mod,control=settings)

+#' 

 fitZig <-

 function(obj,mod,zeroMod=NULL,useS95offset=TRUE,control=zigControl()){

@@ -109,4 +123,4 @@ function(obj,mod,zeroMod=NULL,useS95offset=TRUE,control=zigControl()){

 		dat = list(fit=fit$fit,countResiduals=fit$residuals,

 				   z=z,eb=eb,taxa=rownames(obj),counts=y,zeroMod =mmZero,stillActive=stillActive,stillActiveNLL=stillActiveNLL,zeroCoef=zeroCoef)

 		return(dat)

-	}

\ No newline at end of file

+	}

@@ -1,18 +1,20 @@

-#' Compute the value of the count density function from the count model residuals.

-#'

-#' Calculate density values from a normal: $f(x) = 1/(\sqrt (2 \pi ) \sigma ) e^-((x - \mu )^2/(2 \sigma^2))$.

-#' Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $\deta_{ij}$ = 1 if $y_{ij}$

-#' is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_{ij} = \pi_j(S_j) \cdot f_{0}(y_{ij}) 

-#' +(1-\pi_j (S_j))\cdot f_{count}(y_{ij};\mu_i,\sigma_i^2)$.

-#' The log-likelihood in this extended model is

-#' $(1−\delta_{ij}) \log f_{count}(y;\mu_i,\sigma_i^2 )+\delta_{ij} \log \pi_j(s_j)+(1−\delta_{ij})\log (1−\pi_j (sj))$.

-#' The responsibilities are defined as $z_{ij} = pr(\delta_{ij}=1 | data)$.

-#'

+#' Compute the value of the count density function from the count model

+#' residuals.

+#' 

+#' Calculate density values from a normal: $f(x) = 1/(sqrt (2 pi ) sigma )

+#' e^-((x - mu )^2/(2 sigma^2))$.  Maximum-likelihood estimates are

+#' approximated using the EM algorithm where we treat mixture membership

+#' $deta_ij$ = 1 if $y_ij$ is generated from the zero point mass as latent

+#' indicator variables. The density is defined as $f_zig(y_ij = pi_j(S_j) cdot

+#' f_0(y_ij) +(1-pi_j (S_j))cdot f_count(y_ij;mu_i,sigma_i^2)$. The

+#' log-likelihood in this extended model is $(1-delta_ij) log

+#' f_count(y;mu_i,sigma_i^2 )+delta_ij log pi_j(s_j)+(1-delta_ij)log (1-pi_j

+#' (sj))$. The responsibilities are defined as $z_ij = pr(delta_ij=1 | data)$.

+#' 

+#' 

 #' @param residuals Residuals from the count model.

 #' @param log Whether or not we are calculating from a log-normal distribution.

 #' @return Density values from the count model residuals.

-#'

-#' @name getCountDensity

 #' @seealso \code{\link{fitZig}}

 getCountDensity <-

 function(residuals, log=FALSE){

@@ -1,15 +1,19 @@

-#' Calculate the relative difference between iterations of the negative log-likelihoods.