technical_replicates.Rmd

---
output:
  html_document: default
  pdf_document: default
---



***A Case for Technical Replicates in Three Parts***


This post will argue that including technical replicates -- repeated measurements on identical samples -- is worthwhile in NGS experiments, and it will guide you through one way to use them productively. In that endeavor, we will look at a human 16S dataset to get an overall sense of the scale of variability between replicates, and then we will explore how this variability may affect common analyses and what we can do about it. This is a fairly long post, so it will be posted in three parts: in the first part, we will examine measurement error in an important dataset; in the second, we will explore how this error may impact inferential analyses; and in the third, we will present a solution that partially(!) addresses measurement error in this context.

**Part 1: A (Taxonomically) Broad View of Variability Across Technical Replicates**

As with many measurements, the estimates of sample composition we obtain from 16S sequencing are less than perfectly precise. That is, repeatedly 16S sequencing on specimens from the same sample will typically give us different estimates of sample composition, even if we use the same sequencing and bioinformatics protocols. For an example of this, we can turn to data published by the Microbiome Quality Control Consortium Project (MBQC)$^1$, which assessed technical variation in 16S sequencing data generated by many different sequencing laboratories and under many bioinformatics protocols. To start, let's consider technical replicates of a single sample (a stool sample from a healthy 36-year old man) produced by a single sequencing laboratory (wet laboratory B) and a single bioinformatics laboratory (dry laboratory 2).  (For those of you familiar with the MBQC study design, we are only looking at specimens extracted and sequenced at the same laboratory -- samples extracted at a central lab prior to shipment for sequencing are excluded here.) The plot below shows measured phylum-level compositions for each replicate.


```{r, echo = F}
rm(list = ls())
suppressWarnings(suppressMessages(library(tidyr)))
suppressWarnings(suppressMessages(library(magrittr)))
suppressWarnings(suppressMessages(library(ggplot2)))
suppressWarnings(suppressMessages(library(RDS)))
suppressWarnings(suppressMessages(library(dplyr)))
suppressWarnings(suppressMessages(library(nlme)))
suppressWarnings(suppressMessages(library(viridis)))
suppressWarnings(suppressMessages(library(stringi)))
suppressWarnings(suppressMessages(library(RColorBrewer)))

```

```{r, echo = F, eval = F}

true_value <- c(3,3)

set.seed(0)
axis_1 <-  rep(true_value[1],100) + c(rnorm(50,sd =.03), rnorm(50,sd =.3)) + rep(c(rep(0,25),rep(-.5,25)),2)
axis_2 <-  rep(true_value[2],100) + c(rnorm(50,sd =.03), rnorm(50,sd = .3)) + rep(c(rep(0,25),rep(.75,25)),2)

bias_var <- data.frame(axis_1 = axis_1,
                       axis_2 = axis_2,
                       Bias = rep(rep(c("Less Bias","More Bias"), each = 25),2),
                       Variability = rep(c("Less Variability","More Variability"), each = 50),
                       tv  = "Measurements")


true_value_df <- data.frame(axis_1 = 3, axis_2 = 3, tv = "True Value")

bias_var$Bias <- factor(bias_var$Bias, levels = c("Less Bias","More Bias"))


bias_var %>%
  ggplot() + 
  geom_point(data = true_value_df, aes(x = axis_1, y = axis_2, color = tv), size = 5, alpha = .5) +
  geom_point(aes(x = axis_1, y = axis_2, color = tv), alpha = .25) + 
  facet_grid(Bias~Variability) + 
  scale_color_brewer(palette = "Dark2") + 
  xlab("(Arbitrary) Measurement Axis 1") + 
  ylab("(Arbitrary) Measurement Axis 2") +
  theme_bw() + 
  theme(plot.caption = element_text(hjust = 0, size = 8)) +
  labs(caption = "Simulated example data in four settings: \n - measurements are about equal to the true values of the measurement quantity (low bias, low variability) \n - measurements vary but are on average about equal to the true value  of the measurement quantity (higher variability, low bias) \n - measurements are consist but differ from true value of measurement quantity (higher bias, low variability) \n - measurements vary and on average differ from true value of measurement quantity (higher bias, higher variability)",
       color = "")


```



  
```{r, echo = F, warning = F}

get_lab_data <- function(labname, level = "genus"){
  lab_test <- readRDS(
    paste("./wetlab_files/",
labname,"_test_",level,"_54332",sep = ""))
  lab_train <- readRDS(
    paste("./wetlab_files/",
labname,"_train_",level,"_54332",sep = ""))
  
  return(rbind(lab_test,lab_train))
}

hlb_phylum_data <- get_lab_data("HL-B", "phylum")


hlb_phylum_filtered <- hlb_phylum_data %>%
  filter(dry_lab == "BL-2") %>%
  filter(specimen == "sample4") 

#calculate proportions
hlb_phylum_filtered[,grepl("k__",colnames(hlb_phylum_filtered), fixed = T)] <- 
hlb_phylum_filtered[,grepl("k__",colnames(hlb_phylum_filtered), fixed = T)]  %>%
apply(1, function(x) x/sum(x)) %>% t()

# The palette with grey:

hlb_phylum_filtered %>%
  select(starts_with("Bioinf")|starts_with("k__")) %>%
  mutate(Replicate = as.character(1:nrow(hlb_phylum_filtered))) %>%
  select(-Bioinformatics.ID) %>%
  pivot_longer(cols = starts_with("k__"),
               names_to = "phylum") %>%
  mutate(Phylum = phylum %>%
           sapply(function(x) strsplit(x,"p__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
  filter(value >0) %>%
  ggplot() + 
  geom_bar(aes(x = Replicate,y = value, fill = Phylum),
           stat="identity") +
  ylab("Observed Phylum-level Composition") +
  scale_fill_viridis(discrete = TRUE,option = "C") +
  labs(caption = "Technical replicate measurements of a stool sample from the same healthy 36-year-old man, \nanalyzed by a single sequencing and a single bioinformatics laboratory with sequencing \nand bioinformatics protocol held constant.") +
  theme_bw() + 
   theme(plot.caption = element_text(hjust = 0)) 

```

To some extent, all of these replicates agree -- most 16S reads are attributed to Firmicutes, unsurprisingly. However, we can also see that there is a fair amount of variation between replicates: for example, should we believe that the Firmicutes:Bacteroidetes ratio in this sample is in the neighborhood of 2 or 3 to 1, as replicate 2 might have us think, or is it greater than 10 to 1, as (just eyeballing it), replicate 1 suggests? In any case, we can at a minimum conclude that these replicate measurements cannot all reflect the actual composition of the sample with equal fidelity.

In a moment, we'll look at genus-level data on the same replicates, but first, I'd like to emphasize that because these data are repeat measurements on aliquots from the same sample, if we accept that the MBQC produced uniform aliquots (I do), then **the variation we see here is primarily an artifact of sample preparation, sequencing, and bioinformatics**. That is, in looking at 16S data, we observe the true sample composition plus some noise that I'll call measurement error. Since we're only looking at measurements on a single sample, all (or at least most) of the differences we see between replicates should be due to noise (i.e., measurement error). We can see see this variation because this laboratory sequenced technical replicates, but there is little reason to expect substantially lower noise when technical replicates are not included in a sequencing experiment -- this only prevents us from *observing* variation due to noise.

Now let's take a look at the same replicates at the genus level. To clean things up a bit, the plot below only shows genera that accounted for at least 0.1% of reads in some replicate and that were assigned to a known genus by bioinformatics (i.e., we've thrown out unclassified reads for the time being). 


```{r, echo = F, fig.height = 4}


hlb_genus_data <- get_lab_data("HL-B", "genus")


hlb_genus_filtered <- hlb_genus_data %>%
  filter(dry_lab == "BL-2") %>%
  filter(specimen == "sample4") 

#calculate proportions
hlb_genus_filtered[,grepl("k__",colnames(hlb_genus_filtered), fixed = T)] <- 
hlb_genus_filtered[,grepl("k__",colnames(hlb_genus_filtered), fixed = T)]  %>%
apply(1, function(x) x/sum(x)) %>% t()

# The palette with grey:

genera <- hlb_genus_filtered %>%
  select(starts_with("Bioinf")|starts_with("k__")) %>%
  mutate(Replicate = as.character(1:nrow(hlb_genus_filtered))) %>%
  select(-Bioinformatics.ID) %>%
  pivot_longer(cols = starts_with("k__"),
               names_to = "genus") %>%
  mutate(Genus = genus %>%
           sapply(function(x) strsplit(x,"g__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
  filter(Genus != "unclassified") %>%
  group_by(Genus) %>%
  summarize(max_prev = max(value)) %>%
  filter(max_prev > .001 )

genera$alpha_guide <- rep(c(0,1),11)



 for_genus_bar <- hlb_genus_filtered %>%
  select(starts_with("Bioinf")|starts_with("k__")) %>%
  mutate(Replicate = as.character(1:nrow(hlb_genus_filtered))) %>%
  select(-Bioinformatics.ID) %>%
  pivot_longer(cols = starts_with("k__"),
               names_to = "genus") %>%
  mutate(Genus = genus %>%
           sapply(function(x) strsplit(x,"g__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
  filter(value >0) %>%
  filter(Genus %in% genera$Genus) %>%
    mutate(Family = genus %>%
           sapply(function(x) strsplit(x,"f__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) strsplit(x,"g__", fixed = TRUE)[[1]][1])  %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
   mutate(alpha_guide = 
            sapply(Genus,
                   function(x) genera$alpha_guide[genera$Genus == x])) %>%
  mutate(taxon = stri_join(Genus, "\n(", Family, ")", sep  = "")) 
 

 for_genus_bar %>%
  ggplot() + 
  geom_bar(aes(x = Replicate,y = value, fill = taxon, 
               alpha = as.factor(alpha_guide)),
           stat="identity", position = "fill") +
   geom_text(aes(x = Replicate, y = c(.07,.25,.15,.18), label = Genus),
            data = for_genus_bar %>% filter(Genus == "Prevotella"),
            color = "white",
            size = 2) +
      geom_text(aes(x = Replicate, y = c(.68,.85,.85,.82), label = Genus),
            data = for_genus_bar %>% filter(Genus == "Anaerostipes"),
            color = "white",
            size = 2) +
  ylab("") +
  guides(fill = guide_legend(title = "Genus (Family)",
                             override.aes = list(size = 2,
                                                 alpha = 
                                                   rep(
                                                     c(.5,1),
                                                     11
                                                   ))),
         alpha = FALSE) +
  scale_fill_viridis(discrete = TRUE,option = "C") +
   scale_alpha_manual(values = c(.5,1)) +
  theme_bw() +
  theme(legend.text=element_text(size=7),
        plot.caption = element_text(hjust = 0)) +
  labs(caption = "Technical replicate measurements of a stool sample from a healthy 36-year-old man, \nanalyzed by a single sequencing and a single bioinformatics laboratory with sequencing \nand bioinformatics protocol held constant.") 


```

Again, we observe substantial variability between replicates: in replicate 1, Prevotella (light peach labeled in white), accounts for less than 10% of reads, whereas in replicate 2 almost half of reads are attributed to Prevotella. This plot also highlights a difficulty of viewing relative abundance data as proportions: errors in a single taxon affect all other taxa because we divide by the sum of all reads in all taxa to obtain an estimated proportion. It's not clear from this plot, for example, whether the variability we see in Anaerostipes (dark purple labeled in white) is due to variability in detection of Anaerostipes or due to variability in detection of Prevotella. That is, perhaps the apparent differences in the proportion of reads attributed to Faecalibacterium across replicates are primarily the result of the large variation in Prevotella.

This issue of propagation of error across taxa is in fact a fundamental part of this kind of data, and not just a peculiarity of this dataset. Since 16S sequencing provides information about relative, but not absolute, abundances, we cannot completely disentangle variation in one taxon from variation in another. That is, we have information about abundances in one taxon relative to abundances in others. This poses a somewhat tricky problem when our measurements contain error -- how should we best view and analyze data when summary statistics for all taxa may be affected by measurement error in a single taxon? For instance, suppose that we knew that replicates 2, 3, and 4 above were contaminated with Prevotella from an outside source -- what effect would this have on how we might compare those replicates with replicate 1? What genera are affected if we do not take this information into account?

In the next part of this post, we will look at one way of simplifying the problem of analyzing relative abundance data and consider how measurement error may impact our analyses in this case.

\textbf{References}

1. Rashmi Sinha et al., “Assessment of Variation in Microbial Community Amplicon Sequencing by the Microbiome Quality Control (MBQC) Project Consortium,” Nature Biotechnology 35, no. 11 (2017): pp. 1077-1086, https://doi.org/10.1038/nbt.3981.


***Appendix for the Doubtful***

Should you perhaps be thinking that I have pulled the wool over your eyes by cherry-picking a "bad" sequencing laboratory, here is a multi-laboratory plot of read proportions by phylum for the specimen examined in this post: 

```{r,echo = F, fig.width = 6.5.5}
all_data <- lapply(c(
  "HL-A",
  "HL-B",
  "HL-C",
  "HL-E",
  "HL-F_1",
                "HL-H",
                "HL-I",
                "HL-J",
                "HL-K",
                "HL-L"),
              function(x) get_lab_data(x, "phylum"))

all_data %<>% (function(x) do.call(rbind, x))


phylum_filtered <- all_data %>%
  filter(dry_lab == "BL-2") %>%
  filter(specimen == "sample4") 

#calculate proportions
phylum_filtered[,grepl("k__",colnames(phylum_filtered), fixed = T)] <- 
phylum_filtered[,grepl("k__",colnames(phylum_filtered), fixed = T)]  %>%
apply(1, function(x) x/sum(x)) %>% t()


nb.cols <- 18
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

phylum_filtered %>%
  group_by(blinded_lab) %>%
  select(starts_with("Bioinf")|starts_with("k__"), blinded_lab) %>%
  mutate(Replicate = as.character(1:length(unique(Bioinformatics.ID)))) %>%
  ungroup %>%
  select(-Bioinformatics.ID) %>%
  pivot_longer(cols = starts_with("k__"),
               names_to = "phylum") %>%
  mutate(Phylum = phylum %>%
           sapply(function(x) strsplit(x,"p__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
  filter(value >0) %>%
  mutate(wetlab = sapply(blinded_lab, function(x) substr(x, 4, nchar(x)))) %>%
  ungroup() %>%
  ggplot() + 
  geom_bar(aes(x = as.factor(Replicate),y = value, fill = Phylum),
           stat="identity",
           width = .9) +
  ylab("Observed Phylum-level Composition") +
  scale_fill_manual(values = mycolors) +
  labs(caption = "Technical replicate measurements of a stool sample from the same healthy 36-year-old man, as analyzed by \nmultiple sequencing laboratories and a single bioinformatics laboratory with bioinformatics protocol held constant.") +
  theme_bw() + 
  facet_grid (.~ wetlab, scales = "free_x", space = "free_x") +
   theme(plot.caption = element_text(hjust = 0),
         legend.position = "bottom"
         ) +
  xlab("Sequencing Wet Laboratory and Replicate")
```

By my reckoning, sequencing laboratory B is fairly middling insofar as within-laboratory phylum-level variability in this particular specimen is concerned. However, if you remain unconvinced, we will revisit issues arising from within- and between-laboratory variation in measurements (and beyond!) and continue to examine data across sequencing laboratories in the next two parts of this post. (So stay tuned.)

***Yet Another Appendix***

It occurs to me as well that you may be wondering why I didn't start with laboratory J, which seems to have quite low variation from replicate to replicate on the specimen we've examined in this post. The short answer is that we will return to laboratory J (and the other laboratories) in subsequent posts. A slightly longer answer is that while low techical variation (i.e., low variability between replicate measurements on the same specimen) is _a_ virtue, it is not the _only_ virtue when it comes to measurement. We might, for example, want to know how close our measurements are to the truth (I certainly do). We can't directly ascertain that with the specimen we've been looking at because we don't know its true composition, but the MBQC wisely also chose to include mock specimens of known composition. Here is a plot of measured phylum-level relative abundances (i.e., proportion of reads assigned to each detected phylum) in the fecal mock community, along with the true composition of the fecal mock (indicated by sequencing laboratory "T" at the right hand of the figure):

```{r, echo = F, fig.width = 8.5}

fecal_truth <- data.frame(family = 
                            c("anaeroplasmataceae",
                              "bifidobacteriaceae",
                              "coriobacteriaceae",
                              "desulfovibrionaceae",
                              "enterobacteraceae",
                              "fusobacteraceae",
                              "lachnospiraceae",
                              "lactobacillaceae",
                              "oxalobacteraceae",
                              "paenibacillaceae",
                              "porphyromonadaceae",
                              "bacteroideaceae",
                              "propionibacteriaceae",
                              "rikenellaceae",
                              "ruminococcaceae",
                              "streptococcaceae",
                              "synergistetes"),
                          Phylum = 
                            c("Tenericutes",
                              "Actinobacteria",
                              "Actinobacteria",
                              "Proteobacteria",
                              "Proteobacteria",
                              "Fusobacteria",
                              "Firmicutes",
                              "Firmicutes",
                              "Proteobacteria",
                              "Firmicutes",
                              "Bacteroidetes",
                              "Bacteroidetes",
                              "Actinobacteria",
                              "Bacteroidetes",
                              "Firmicutes",
                              "Firmicutes",
                              "Synergistetes"),
                          value = c(3.5,
                                    6,
                                    6,
                                    2.25,
                                    6,
                                    12,
                                    12,
                                    7.5,
                                    6,
                                    3.5,
                                    6,
                                    6,
                                    6,
                                    3.5,
                                    1.5,
                                    6,
                                    2))
                             
fecal_truth$value %<>% (function(x) x/sum(x))  
fecal_truth$wetlab <- "T"
fecal_truth$Phylum %<>% as.character()
fecal_truth$Replicate <- as.character(1)
fecal_truth %<>% select(-family)

phylum_filtered <- all_data %>%
  filter(dry_lab == "BL-2") %>%
  filter(specimen_type_collapsed == "Fecal artificial colony") 

phylum_filtered[,grepl("k__",colnames(phylum_filtered), fixed = T)] <- 
phylum_filtered[,grepl("k__",colnames(phylum_filtered), fixed = T)]  %>%
apply(1, function(x) x/sum(x)) %>% t()


phylum_filtered %>%
  group_by(blinded_lab) %>%
  select(starts_with("Bioinf")|starts_with("k__"), blinded_lab) %>%
  mutate(Replicate = as.character(1:length(unique(Bioinformatics.ID)))) %>%
  ungroup %>%
  select(-Bioinformatics.ID) %>%
  pivot_longer(cols = starts_with("k__"),
               names_to = "phylum") %>%
  mutate(Phylum = phylum %>%
           sapply(function(x) strsplit(x,"p__", fixed = TRUE)[[1]][2]) %>%
           sapply(function(x) substr(x,1,nchar(x) -1))) %>%
  filter(value >0) %>%
  mutate(wetlab = sapply(blinded_lab, function(x) substr(x, 4, nchar(x)))) %>%
  ungroup() %>% 
  select(Phylum,value, wetlab, Replicate) %>%  
  (function(x) rbind(x, fecal_truth)) %>%
  ggplot() + 
  geom_bar(aes(x = as.factor(Replicate),y = value, fill = Phylum),
           stat="identity",
           width = .9) +
  ylab("Observed Phylum-level Composition") +
  scale_fill_manual(values = mycolors) +
  labs(caption = "Technical replicate measurements on a fecal mock community, as analyzed by multiple sequencing laboratories and a single \nbioinformatics laboratory with bioinformatics protocol held constant. Laboratory \"T\" (far right-hand column) gives true \nrelative abundances in the specimen.") +
  theme_bw() + 
  facet_grid (.~ wetlab, scales = "free_x", space = "free_x") +
   theme(plot.caption = element_text(hjust = 0),
         legend.position = "bottom"
         ) +
  xlab("Sequencing Wet Laboratory and Replicate")

```

I will let you draw your own conclusions here, but I hope we can agree that the relationship between measure and true relative abundances is, while far from arbitrary, not particularly straightforward.