Practical Microbiome Analysis in R · Lesson 3
Lesson 2 matched the microbial count matrix, taxonomy lookup and sample metadata. Its final audit contained 875 sample records, 455 features and 43 children, with no unmatched identifiers. We now use those matched objects to prepare age and categorical variables and examine how sampling records are distributed across children and days.
This lesson keeps every record and every original metadata column. It adds five clearly named columns and saves a new checkpoint. The count matrix and taxonomy remain unchanged. No package installation, download, archive extraction or abundance import is needed.
Before drawing a graph or comparing groups, check that ages can be used in calculations and that sample records are distinguished from children and sampling days.
Lessons in this series
- Getting Started with R and RStudio for Microbiome Analysis
- Importing and Checking ECAM Microbiome Data in R — previous lesson.
- Preparing ECAM Metadata for Analysis in R — current lesson.
Before you begin
Keep the MicrobesWorld_ECAM project and R session open after Lesson 2. Save the companion script as scripts/lesson03.R, open it in RStudio and click Source. The script uses ecam_lesson02, the object returned by the complete Lesson 2 script. It also accepts bundle, the equivalent object created by following Lesson 2's teaching blocks individually.
Use one execution route: Source the companion script, or work through the R code blocks below in order. The script wraps the steps in a local block and returns the successful result as ecam_lesson03. Running the teaching blocks individually creates prepared. Both routes write the same new checkpoint. Copy only the commands in an R code box, not expected output, console prompts or formatting marks. On small screens, scroll code boxes and tables sideways to see complete lines. Stop at the first error rather than executing a later success message separately.
If you have restarted R, restore the checkpoint that Lesson 2 actually created with this one recovery command:
Recovery code — only after restarting R
ecam_lesson02 <- readRDS(file.path("results", "ecam_matched_lesson02.rds"))Skip that command when you are continuing directly from Lesson 2. It restores previously saved objects; it does not repeat the source-data import. If the file is missing, stop and check that Lesson 2 completed and that you opened the correct project.
Step 1. Continue with the matched objects
The first block selects the existing Lesson 2 object and checks the expected source version and table relationships. These are short checks of the handoff between lessons, not another import or a repeat of the complete identifier audit.
R code
if (exists("ecam_lesson02", envir = .GlobalEnv, inherits = FALSE)) {
previous <- get("ecam_lesson02", envir = .GlobalEnv)
} else if (exists("bundle", envir = .GlobalEnv, inherits = FALSE)) {
previous <- get("bundle", envir = .GlobalEnv)
} else {
stop("The Lesson 2 objects are missing. Restore results/ecam_matched_lesson02.rds before continuing.")
}
stopifnot(is.list(previous),
all(c("counts", "taxonomy", "metadata", "source_commit") %in% names(previous)),
identical(previous$source_commit, "563c8a7e39b1704baea0b35ac3d4774914559af3"),
is.matrix(previous$counts), is.numeric(previous$counts),
identical(dim(previous$counts), c(875L, 455L)),
is.data.frame(previous$metadata), ncol(previous$metadata) == 12L,
is.data.frame(previous$taxonomy),
identical(rownames(previous$counts), previous$metadata$SampleID),
identical(colnames(previous$counts), previous$taxonomy$FeatureID),
all(c("studyid", "day_of_life", "delivery", "sex", "diet") %in% names(previous$metadata)))
metadata_raw <- previous$metadata
metadata <- metadata_raw
cat("Continuing from the matched Lesson 2 objects.\n")Expected output
Continuing from the matched Lesson 2 objects.
The working metadata are a copy of the matched data, not the earlier unmatched or alphabetically sorted table. In the complete script, this copy stays inside the local execution block, so the objects from Lessons 1 and 2 are preserved.
Step 2. Look for missing and blank entries
R's missing value NA and an empty text field are different things. A field containing spaces can also look empty on screen. The following audit counts NA values separately from empty or whitespace-only text in each of the 12 original columns. It reports what is present without replacing values.
R code
missing_audit <- data.frame(
variable = names(metadata_raw),
missing_NA = vapply(metadata_raw, function(x) sum(is.na(x)), integer(1)),
blank_text = vapply(metadata_raw, function(x) {
sum(!is.na(x) & trimws(as.character(x)) == "")
}, integer(1)),
row.names = NULL
)
print(missing_audit, row.names = FALSE)
stopifnot(!anyNA(metadata$SampleID), !anyNA(metadata$studyid),
all(trimws(metadata$SampleID) != ""),
all(trimws(metadata$studyid) != ""),
anyDuplicated(metadata$SampleID) == 0L)Expected output
variable missing_NA blank_text
SampleID 0 0
antiexposedall 0 0
day_of_life 0 0
delivery 0 0
diet 0 0
diet_3 0 0
mom_child 0 0
month 0 0
month_of_life 0 0
sample_summary 0 0
sex 0 0
studyid 0 0All counts are zero for this source version. This result covers NA and blank text only. It does not prove that every source annotation is complete or correct: a source-specific code could represent an unknown or unavailable value without being blank. We therefore inspect the observed category labels in Step 4 rather than inventing missing-value definitions.
The SampleID and studyid checks are prerequisites for the new grouping operations. Missing participant labels must not be silently dropped by a grouped summary.
Step 3. Add age as a numeric variable
Lesson 1 deliberately imported the metadata as text to preserve their source representation. That means day_of_life is still a character column. We need a numeric age column for arithmetic and later plots, while keeping the original text available for comparison.
For this particular source, the code requires nonmissing whole-day strings before conversion. It then checks that the resulting values are finite, nonnegative whole numbers. A decimal, a negative value or an unexpected text label stops this source-specific check rather than being silently accepted.
R code
age_text <- trimws(as.character(metadata$day_of_life))
stopifnot(!anyNA(age_text), all(grepl("^[0-9]+$", age_text)))
metadata$age_days <- as.numeric(age_text)
stopifnot(all(is.finite(metadata$age_days)), all(metadata$age_days >= 0),
all(metadata$age_days == floor(metadata$age_days)))
print(data.frame(
source_class = class(metadata$day_of_life),
new_class = class(metadata$age_days),
minimum_days = min(metadata$age_days),
maximum_days = max(metadata$age_days)
), row.names = FALSE)Expected output
source_class new_class minimum_days maximum_days
character numeric 0 864age_days now supports numerical calculations, and day_of_life remains character. Do not convert a factor directly with as.numeric() to recover age: that can return internal level codes rather than the displayed values. Here, conversion starts from checked character strings.
The observed maximum is 864 days. That is a fact about this processed file, not an instruction to truncate follow-up at two years or to relabel day values as months. We preserve all recorded ages at this stage.
Step 4. Represent categories explicitly
A participant identifier is a label, even when it looks like a number. Treating studyid as a numeric predictor could incorrectly impose a numerical relationship between children. We add child_id as a factor and make categorical copies of delivery, sex and diet. The original columns and their codes remain intact.
Factor levels follow the order in which labels occur in the matched metadata. That order is an implementation choice, not a biological ranking or a prespecified reference group for a model. A later modeling lesson must choose and explain its reference categories.
R code
metadata$child_id <- factor(metadata$studyid, levels = unique(metadata$studyid))
metadata$delivery_group <- factor(metadata$delivery, levels = unique(metadata$delivery))
metadata$sex_group <- factor(metadata$sex, levels = unique(metadata$sex))
metadata$diet_code <- factor(metadata$diet, levels = unique(metadata$diet))
stopifnot(nlevels(metadata$child_id) == 43L,
identical(as.character(metadata$child_id), metadata$studyid),
identical(as.character(metadata$diet_code), metadata$diet))
category_variables <- c("delivery", "sex", "diet", "diet_3", "antiexposedall", "mom_child")
category_counts <- do.call(rbind, lapply(category_variables, function(v) {
z <- table(metadata_raw[[v]], useNA = "ifany")
data.frame(variable = v, source_value = names(z),
sample_records = as.integer(z), stringsAsFactors = FALSE)
}))
rownames(category_counts) <- NULL
print(category_counts, row.names = FALSE)Expected output
variable source_value sample_records
delivery Cesarean 349
delivery Vaginal 526
sex Female 316
sex Male 559
diet bd 636
diet fd 239
diet_3 bd 94
diet_3 eb 542
diet_3 fd 239
antiexposedall n 569
antiexposedall y 306
mom_child C 875These frequencies count sample records, not children. For example, 636 records carry diet code bd, but that is not a count of 636 participants. Children contributed different numbers of samples.
The frequency table retains source codes such as bd, fd, eb, n and y. Turning a column into a factor does not establish the full definition or timing of its annotations. In particular, diet and diet_3 are separate source variables; do not merge them because some labels overlap. We have not recoded diet histories or interpreted antiexposedall as exposure at a particular visit.
The new categorical copies retain their source labels. This makes the R representation explicit while leaving biological definitions and any later eligibility rules open to documented source review.
Step 5. Flag repeated child-day combinations
A unique SampleID does not guarantee a unique sampling day for a child. Two records can have different sample identifiers while sharing studyid and age_days. We inspect that relationship before calculating longitudinal summaries.
The code uses the pair of columns directly, rather than building an ambiguous concatenated identifier. duplicated() detects later occurrences; checking from both directions flags every member of a repeated child-day group, including its first record.
R code
child_day <- metadata[c("studyid", "age_days")]
repeated_child_day <- duplicated(child_day) | duplicated(child_day, fromLast = TRUE)
repeated_records <- metadata[repeated_child_day,
c("SampleID", "studyid", "day_of_life", "age_days"), drop = FALSE]
repeated_summary <- data.frame(
measure = c("Unique child-day combinations", "Repeated child-day combinations",
"Records in repeated child-day groups", "Records beyond one per child-day"),
value = c(nrow(unique(child_day)),
nrow(unique(child_day[repeated_child_day, , drop = FALSE])),
sum(repeated_child_day), sum(duplicated(child_day)))
)
print(repeated_summary, row.names = FALSE)Expected output
measure value
Unique child-day combinations 724
Repeated child-day combinations 149
Records in repeated child-day groups 300
Records beyond one per child-day 151There are 724 distinct child-day pairs. Of those pairs, 149 have multiple records, containing 300 records altogether. The figure 151 counts records beyond one per pair: 875 minus 724. It is not a recommended deletion count.
The output does not identify why a child-day pair occurs more than once. It does not establish whether those records are separate specimens, technical replicates or another kind of source record. We preserve every SampleID and save the flagged records for inspection. Neither dropping the later row nor summing its counts is justified by matching ages alone.
This distinction matters for later analysis: repeated observations from the same child must not be treated automatically as independent participants, and multiple records on a day must not be mistaken for multiple distinct time points.
Step 6. Separate record counts from observed sampling days
For each child, count both the number of records and the number of distinct recorded ages. The first can exceed the second. The summary also reports the earliest and latest observed day; those endpoints do not imply continuous observation between them.
split() groups row positions by child_id. Within each group, unique() counts distinct ages, and min() and max() describe the observed time range. We print only the first six child summaries; the saved CSV contains all 43.
R code
child_rows <- split(seq_len(nrow(metadata)), metadata$child_id, drop = TRUE)
child_summary <- do.call(rbind, lapply(child_rows, function(i) {
d <- metadata$age_days[i]
data.frame(studyid = metadata$studyid[i[1L]], sample_records = length(i),
distinct_days = length(unique(d)), first_day = min(d), last_day = max(d),
stringsAsFactors = FALSE)
}))
rownames(child_summary) <- NULL
print(head(child_summary, 6), row.names = FALSE)
print(data.frame(
children = nrow(child_summary),
children_with_one_distinct_day = sum(child_summary$distinct_days == 1L),
minimum_distinct_days = min(child_summary$distinct_days),
maximum_distinct_days = max(child_summary$distinct_days)
), row.names = FALSE)Expected output
studyid sample_records distinct_days first_day last_day
1 32 30 0 855
2 22 19 0 856
4 15 15 0 366
5 24 21 1 855
7 28 25 1 853
8 18 17 1 795
children children_with_one_distinct_day minimum_distinct_days
43 1 1
maximum_distinct_days
30Child 1 illustrates the difference: 32 records correspond to 30 distinct days, with observed ages from day 0 to day 855. Across the dataset, one child has records on only one distinct day, even though Lesson 1 showed that every child has at least two records. Counting rows alone would miss that distinction.
These findings describe the sampling structure. They do not by themselves select an analysis cohort or justify excluding a child. A later longitudinal method may require multiple distinct times, but that requirement must be stated for the specific analysis. No child is removed here.
The printed child order follows the factor levels inherited from the matched table, not a numeric ranking of participant identifiers. Creating the summary does not reorder the working metadata or the count matrix.
Step 7. Save the prepared objects
The new checkpoint retains the count matrix, taxonomy and source-version details from Lesson 2. Its metadata contain the original 12 columns plus age_days, child_id, delivery_group, sex_group and diet_code. metadata_raw preserves the original matched metadata separately.
Before saving, the code confirms that every original metadata column and the sample order are unchanged. After saving, it reloads the RDS file and compares the complete object. The final message appears only after those checks pass.
R code
stopifnot(identical(metadata[names(metadata_raw)], metadata_raw),
identical(metadata$SampleID, rownames(previous$counts)),
nrow(metadata) == 875L, ncol(metadata) == 17L,
dir.exists("results"))
prepared <- previous
prepared$metadata <- metadata
prepared$metadata_raw <- metadata_raw
prepared$metadata_audit <- list(missing = missing_audit, categories = category_counts,
repeated_summary = repeated_summary,
repeated_records = repeated_records,
children = child_summary)
prepared$lesson03_note <- paste(
"Added age_days and four factors; preserved all 12 source metadata columns.",
"No records removed or combined; no counts transformed; no models fitted."
)
saveRDS(prepared, file.path("results", "ecam_prepared_lesson03.rds"))
write.csv(metadata, file.path("results", "lesson03_metadata_prepared.csv"), row.names = FALSE)
write.csv(missing_audit, file.path("results", "lesson03_missing_audit.csv"), row.names = FALSE)
write.csv(category_counts, file.path("results", "lesson03_category_counts.csv"), row.names = FALSE)
write.csv(repeated_summary, file.path("results", "lesson03_repeated_day_summary.csv"), row.names = FALSE)
write.csv(repeated_records, file.path("results", "lesson03_repeated_child_day_records.csv"), row.names = FALSE)
write.csv(child_summary, file.path("results", "lesson03_child_sampling_summary.csv"), row.names = FALSE)
writeLines(capture.output(sessionInfo()), file.path("results", "lesson03_sessionInfo.txt"))
stopifnot(identical(readRDS(file.path("results", "ecam_prepared_lesson03.rds")), prepared),
identical(prepared$counts, previous$counts),
identical(prepared$taxonomy, previous$taxonomy))
cat("Saved Lesson 3: 875 records, 17 metadata columns; counts and taxonomy unchanged.\n")Expected output
Saved Lesson 3: 875 records, 17 metadata columns; counts and taxonomy unchanged.
The script creates eight files in results:
| File | Purpose |
|---|---|
| ecam_prepared_lesson03.rds | Prepared R objects, original metadata and audit tables. |
| lesson03_metadata_prepared.csv | 875 rows and 17 metadata columns for inspection. |
| lesson03_missing_audit.csv | Missing-value and blank-text counts for the original columns. |
| lesson03_category_counts.csv | Observed category labels and sample-record frequencies. |
| lesson03_repeated_day_summary.csv | The four sampling-day summary measures. |
| lesson03_repeated_child_day_records.csv | All 300 records in repeated child-day groups. |
| lesson03_child_sampling_summary.csv | One descriptive row for each of the 43 children. |
| lesson03_sessionInfo.txt | The R and package environment record. |
Use the RDS checkpoint to continue later lessons: it retains factor classes and the full object structure. CSV files are useful inspection copies, but a fresh CSV import does not automatically preserve factor levels. Rerunning Lesson 3 replaces its generated outputs; it does not overwrite the Lesson 2 checkpoint or either original source file.
What is ready for the next lesson?
The prepared dataset still has 875 records, 455 taxonomic features and 43 children. Age is available as a checked numeric column, participant identity is explicitly categorical, and the sampling-day audit is saved alongside the data. This is preparation for descriptive plots, not evidence for a dietary effect or a completed statistical analysis.
The next lesson will use ecam_lesson03 to examine sampling coverage and the distribution of processed-table count totals. Those totals must be described as counts retained in this processed source, not absolute microbial loads. We will also distinguish a graph describing sampling from a graph testing a biological hypothesis.
Troubleshooting
The Lesson 2 objects are missing: if R was restarted, use the single recovery command near the beginning. If its checkpoint cannot be found, check the project and the successful Lesson 2 save step. Do not redownload the source files as a first response.
A source-version or alignment check fails: stop and use the matched Lesson 2 objects. An arbitrary metadata frame or a previously modified bundle is not an equivalent input.
The age check stops: inspect day_of_life for missing values or unexpected strings. Do not suppress the error or replace a failed conversion with zero.
The output shows repeated child-days: this is an audit finding, not a script error. Inspect the exported record list; no automatic removal is performed.
The Environment does not show a separate metadata object after Source: that is expected. The complete script returns its results in ecam_lesson03. The prepared metadata are ecam_lesson03$metadata; the original matched metadata are ecam_lesson03$metadata_raw.
Previous lesson: Importing and Checking ECAM Microbiome Data in R
Next lesson: Exploring ECAM Sampling Coverage and Count Totals in R — coming soon.
Lesson script
The companion file contains the same seven steps in one execution block. Save it as scripts/lesson03.R and click Source after completing Lesson 2 in the same R session. You only need to download it once.
References
Bokulich, N. A., Dillon, M. R., Zhang, Y., Rideout, J. R., Bolyen, E., Li, H., Albert, P. S., & Caporaso, J. G. (2018). q2-longitudinal: Longitudinal and paired-sample analyses of microbiome data. mSystems, 3(6), e00219-18. https://doi.org/10.1128/mSystems.00219-18
Caporaso Lab. (n.d.). longitudinal-notebooks [Data and analysis repository; commit 563c8a7e39b1704baea0b35ac3d4774914559af3]. GitHub. https://github.com/caporaso-lab/longitudinal-notebooks/tree/563c8a7e39b1704baea0b35ac3d4774914559af3
R Core Team. (n.d.-a). Determine duplicate elements. R documentation. https://stat.ethz.ch/R-manual/R-devel/library/base/html/duplicated.html
R Core Team. (n.d.-b). Divide into groups and reassemble. R documentation. https://stat.ethz.ch/R-manual/R-devel/library/base/html/split.html
R Core Team. (n.d.-c). Factors. R documentation. https://stat.ethz.ch/R-manual/R-devel/library/base/html/factor.html
R Core Team. (n.d.-d). Numeric vectors. R documentation. https://stat.ethz.ch/R-manual/R-devel/library/base/html/numeric.html
R Core Team. (n.d.-e). Serialization interface for single objects. R documentation. https://stat.ethz.ch/R-manual/R-devel/library/base/html/readRDS.html
No comments:
Post a Comment