himel_code.R

###########################
# Load Required Libraries #
###########################

library(vegan) # TSS Normalization
library(data.table) # fread
# Install MetagenomeSeq (EM_ZIG Model)
library(devtools)
# install_github('HCBravoLab/metagenomeSeq', ref='per_feature_em')
library(metagenomeSeq) #FIT-ZIG

####################################################
# Load A Simple Test Data Generated by sparseDOSSA #
####################################################

data<-fread('third_example.pcl', header = TRUE) # sparseDOSSA output
# You left out converting the data table into a matrix
data<-as.data.frame.matrix(data)
metadata<-data[1:10,-1] # extract continuous metadata
bugs<-data[141:190,-1] # extract spiked-in features
metadata<-data.matrix(metadata) # Make numeric
Y<-data.matrix(bugs) # Make numeric

######################
# TSS Normalization  #
######################
Y<- decostand(Y, method="total", MARGIN=2) # Relative Abundance
apply(Y, MARGIN=2, FUN=sum)  # Checking to See if Each Individual Subject Sums to One

# Prepare Data for MetagenmoseSeq
# You need to still specify the library size estimate
# If you don't want any normFactors, just keep them as one
Y<-newMRexperiment(Y,libSize=colSums(data.matrix(bugs)),normFactors=1)
X<-model.matrix(~t(metadata))
# You provided an unestimable design matrix
X<-X[,1:2]
# I am filtering out those without any zeros in at least two samples -
# need to have code that accounts for that
Y2<- Y[-which(rowSums(MRcounts(Y)!=0)>=(ncol(Y)-1)),]
# Per-feature EM
settings = zigControl(per_feature_zeroModel=TRUE)
fit = fitZig(obj = Y2,mod=X,control=settings,useCSSoffset = FALSE)