I'm trying to follow a slight variant of the tutorial to get started from ~ page 9: https://www.bioconductor.org/packages/devel/bioc/vignettes/metagenomeSeq/inst/doc/metagenomeSeq.pdf

Essentially, I need to figure out how to go from: (A) my raw Otu counts; and (B) my "class" assignment (e.g. diseased, healthy) to get (C) differentially abundant Otus.

I'm starting with the built in dataset lungData but I'm getting an error when I'm calculating the normalization factors.

I've made sure there are no empty features (Otus) or empty samples. I've documented my code heavily to serve as a tutorial in case anyone needs to learn how to use the pipeline.

Please help me figure out how to calculate the normalization factors

# Load packages library(metagenomeSeq)

# Paths path_datadirectory = system.file("extdata", package = "metagenomeSeq")

# Load in counts data count_data = loadMeta(file.path(path_datadirectory, "CHK_NAME.otus.count.csv")) dim(count_data$counts) # [1] 1000 78 names(count_data) # [1] "counts" "taxa" class(count_data$counts) # [1] "data.frame"

# Load taxonomy data df_taxonomy = read.delim(file.path(path_datadirectory, "CHK_otus.taxonomy.csv"), stringsAsFactors = FALSE) dim(df_taxonomy) # [1] 1000 10 colnames(df_taxonomy) # [1] "OTU" "Taxonomy" "superkingdom" "phylum" "class" "order" "family" "genus" # [9] "species" "strain"

# Loading metadata df_metadata = loadPhenoData(file.path(path_datadirectory, "CHK_clinical.csv"), tran = TRUE) idx_samples_from_counts = colnames(count_data$counts) idx_samples_from_metadata = rownames(df_metadata) idx_samples = intersect(idx_samples_from_counts, idx_samples_from_metadata)

# Reorder count_data$counts = count_data$counts[,idx_samples] df_metadata = df_metadata[idx_samples,] print(all.equal(rownames(df_metadata), colnames(count_data$counts))) # [1] TRUE

# Annotated DataFrames Adf_metadata = AnnotatedDataFrame(df_metadata) # An object of class 'AnnotatedDataFrame' # rowNames: CHK_6467_E3B11_BRONCH2_PREWASH_V1V2 CHK_6467_E3B11_OW_V1V2 ... CHK_6467_E3B09_BAL_A_V1V2 (78 total) # varLabels: SampleType SiteSampled SmokingStatus # varMetadata: labelDescription Adf_taxonomy = AnnotatedDataFrame(df_taxonomy) # An object of class 'AnnotatedDataFrame' # rowNames: 1 2 ... 1000 (1000 total) # varLabels: OTU Taxonomy ... strain (10 total) # varMetadata: labelDescription obj = newMRexperiment( count_data$counts, phenoData=Adf_metadata, featureData=Adf_taxonomy, ) # MRexperiment (storageMode: environment) # assayData: 1000 features, 78 samples # element names: counts # protocolData: none # phenoData # sampleNames: CHK_6467_E3B11_BRONCH2_PREWASH_V1V2 CHK_6467_E3B11_OW_V1V2 ... CHK_6467_E3B09_BAL_A_V1V2 (78 total) # varLabels: SampleType SiteSampled SmokingStatus # varMetadata: labelDescription # featureData # featureNames: 1 2 ... 1000 (1000 total) # fvarLabels: OTU Taxonomy ... strain (10 total) # fvarMetadata: labelDescription # experimentData: use 'experimentData(object)' # Annotation:

# Filter idx_features_f = which(rowSums(obj) >= 10) idx_samples_f = which(colSums(obj) > 0) obj_f = obj[idx_features_f, idx_samples_f] # MRexperiment (storageMode: environment) # assayData: 12 features, 57 samples # element names: counts # protocolData: none # phenoData # sampleNames: CHK_6467_E3B11_OW_V1V2 CHK_6467_E3B08_OW_V1V2 ... CHK_6467_E3B09_BAL_A_V1V2 (57 total) # varLabels: SampleType SiteSampled SmokingStatus # varMetadata: labelDescription # featureData # featureNames: 77 79 ... 978 (12 total) # fvarLabels: OTU Taxonomy ... strain (10 total) # fvarMetadata: labelDescription # experimentData: use 'experimentData(object)' # Annotation:

# Dimensionality n_samples = ncol(obj_f) n_features = nrow(obj_f) library_size = libSize(obj_f) # `libSize(obj)`` is the same as `colSums(count_data$counts)``

# Normalization p = cumNormStatFast(obj_f)

# Error Error in cumNormStatFast(obj_f) : Warning sample with one or zero features