ECAM R Tutorial · Lesson 5. In Lesson 4, we saw that some children were sampled more often than others. If we want to compare gut microbiome profiles near day 90, which sample should represent each child? In this lesson, we’ll make that choice once, check who was left out, and save the dataset for the next lessons.

Figure 1. Three example sampling timelines from day 60 to day 120, with the sample closest to day 90 highlighted for each child. This figure was created by AI.

Why use one sample per child?

ECAM followed the same children over time. If we put every record into a simple comparison as though it came from a different child, children with more visits would contribute more observations, and we would ignore the dependence between samples from the same child. Taking one record per child gives us a cross-sectional example near day 90. The age window is a tutorial decision, not a rule from the original ECAM paper.

Prepare the ECAM dataset in R

1. Load the Lesson 3 checkpoint

Lesson 4 saved coverage and count summaries. The matched counts, taxonomy and metadata still live in the Lesson 3 checkpoint, so start there in the same RStudio project. No extra packages or downloads are needed. Run the code blocks in order; the console output beneath each block is for comparison.

path <- file.path("results", "ecam_prepared_lesson03.rds")
if (!file.exists(path)) stop("Open the ECAM project with the Lesson 3 checkpoint.")
full <- readRDS(path)
m <- full$metadata
stopifnot(identical(m$SampleID, rownames(full$counts)),
          identical(full$taxonomy$FeatureID, colnames(full$counts)),
          !anyNA(m[c("SampleID", "studyid", "age_days", "delivery")]),
          !anyDuplicated(m$SampleID))
cat(nrow(m), "records from", length(unique(m$studyid)), "children\n")

Expected console output:

875 records from 43 children

The checks compare sample and feature IDs in order, and stop if key metadata are missing or a sample ID is duplicated. The 875 records come from 43 children because the study sampled children repeatedly.

2. Select the sample closest to day 90

Days 60 and 120 both count. For each child, take the record closest to day 90. If two ages are equally close, take the earlier day; if records still tie, take the first sample ID in sorted order. That last rule makes the selection reproducible without treating one tied sample as biologically preferable.

in_window <- m$age_days >= 60 & m$age_days <= 120
eligible <- in_window & m$delivery %in% c("Cesarean", "Vaginal")
candidates <- m[eligible, , drop = FALSE]
ranked <- order(candidates$studyid, abs(candidates$age_days - 90),
                candidates$age_days, candidates$SampleID, method = "radix")
candidates <- candidates[ranked, , drop = FALSE]
keep <- candidates[!duplicated(candidates$studyid), , drop = FALSE]
keep$delivery <- factor(keep$delivery)
print(table(keep$delivery))
print(tapply(keep$age_days, keep$delivery, range))

Expected console output:

Cesarean  Vaginal
      18       24

$Cesarean
[1] 65 99

$Vaginal
[1] 60 98

The age ranges in the output describe the samples actually selected, not the limits of the window. The final groups contain 18 children born by Cesarean delivery and 24 by vaginal delivery. One sample per child addresses repeat sampling in this exercise; differences in age, feeding or antibiotic exposure can still matter when interpreting group comparisons.

3. Account for every record

A sample can be inside the window and still be left out. The audit gives each input record one reason, so you can trace the difference between a sample that was outside the window and an eligible sample that ranked behind another record from the same child.

audit <- m[c("SampleID", "studyid", "age_days", "delivery")]
audit$decision <- "Outside days 60–120"
audit$decision[in_window & !eligible] <- "Delivery group outside comparison"
audit$decision[eligible] <- "Another record selected for child"
audit$decision[m$SampleID %in% keep$SampleID] <- "Selected"
print(table(audit$decision))
omitted <- sort(setdiff(unique(m$studyid), keep$studyid))
print(data.frame(studyid = omitted,
                 records = vapply(omitted, function(id) sum(m$studyid == id), integer(1)),
                 in_window = vapply(omitted, function(id) sum(m$studyid == id & in_window), integer(1)),
                 row.names = NULL))

Expected console output:

Another record selected for child               Outside days 60–120
                               49                               784
                         Selected
                               42
  studyid records in_window
1      23       2         0

Of the 91 eligible records, 42 were kept and 49 ranked behind another record from the same child. Another 784 were outside the age window; none were excluded because of delivery group. Child 23 had two records, and neither fell in the window.

4. Save the analysis dataset

Use the selected sample IDs to pull the matching count rows in the same order as the metadata. Save those rows with the metadata and taxonomy in an RDS file, and save the record-by-record decisions in a CSV you can open separately.

counts <- full$counts[keep$SampleID, , drop = FALSE]
stopifnot(identical(rownames(counts), keep$SampleID),
          !anyDuplicated(keep$studyid), all(rowSums(counts) > 0))
subset <- list(counts = counts, metadata = keep, taxonomy = full$taxonomy,
               source_commit = full$source_commit,
               selection = list(first_day = 60, target_day = 90, last_day = 120),
               audit = audit)
saveRDS(subset, file.path("results", "ecam_analysis_lesson05.rds"))
write.csv(audit, file.path("results", "lesson05_sample_selection.csv"), row.names = FALSE)
cat(nrow(counts), "selected samples by", ncol(counts), "feature columns\n")

Expected console output:

42 selected samples by 455 feature columns

The new dataset has 42 sample rows and all 455 feature columns. A feature column can be present even if all selected samples have a zero for that feature. This step does not filter features, transform counts or calculate diversity. The original longitudinal checkpoint stays available for Lesson 12.

Next lesson: Calculating Alpha Diversity in R

Lesson 6 will open results/ecam_analysis_lesson05.rds to calculate observed feature richness and Shannon diversity for these same samples. Keep results/ecam_prepared_lesson03.rds for the later longitudinal lesson. If you want to revisit an earlier step, browse the ECAM R Tutorial series.

References

Data source: The Caporaso Lab longitudinal notebooks repository hosts the ECAM teaching files used in the series.