E013 Pho4 chimera activity analysis, new host

suppressPackageStartupMessages(require(plotly))
suppressPackageStartupMessages(require(tidyverse))
suppressPackageStartupMessages(require(cowplot))

Goal

• Analyze the latest batches of flow cytometry data to determine the contribution of different (matching) parts of ScPho4 and CgPho4 to their difference in Pho2 dependence.
• Develop an analysis pipeline to perform QC, correction (if needed) and plotting the results.

Data

Merge the three batches

raw <- list(
  "03/15/22" = read_tsv("../data/20220315-EO-chimera-batch-1/20220329-gated-median-out.txt", show_col_types = FALSE),
  "03/22/22" = read_tsv("../data/20220322-EO-chimera-batch-2/20220327-gated-median-out.txt", show_col_types = FALSE),
  "03/29/22" = read_tsv("../data/20220329-EO-chimera-batch-3/20220330-gated-median-out.txt", show_col_types = FALSE)
)
dat <- bind_rows(raw, .id = "date") %>% 
  #separate(col = well, sep = 1, into = c("row", "col")) %>% 
  separate(col = sample, sep = "-", into = c("plasmid", "host", "rep"))

QC

1. Number of events per sample

dat %>% ggplot(aes(x = n_induction)) + geom_histogram(bins = 30, fill = "forestgreen") + facet_grid(host~date) + theme_bw(base_size = 14) #+ panel_border()

Majority of the experimental wells (not blank) have > 10000 events, with the exception of a few samples

dat %>% filter(host != "blank", n_induction < 10000)

Given the distribution of event counts, we will label all samples with < 7000 counts, which would exclude 2 samples. The rationale is that even though some samples had between 7000 - 10000 events, their mean fluorescence values appear in line with the other replicates of the same genotype. I also removed the blank wells and add a flag to mark the plasmids that Emily remade because the original one (used for this batch of data) contains nonsynonymous mutations.

event.th <- 7000

list.mutations <- as.character(c(215, 209, 211, 210, 217, 218, 220, 221, 223, 224, 222, 227, 230, 229, 231, 233, 240, 241, 250:255, 257, 258, 265, 266))

expt <- dat %>% 
  filter(host != "blank") %>% 
  group_by(date, plasmid, host) %>% 
  summarize(n = n(), n_filter = sum(n_induction >= event.th), .groups = "drop") %>% 
  mutate(flag = ifelse(plasmid %in% list.mutations, "mutation", ""))

dat <- dat %>% 
  mutate(
    low_event_count = n_induction < event.th,
    mutation = ifelse(plasmid %in% list.mutations, TRUE, FALSE)
  ) %>% 
  filter(host != "blank")

How many plasmid-host combination have < 3 replicates after removing samples with low event counts?

dat %>% filter(!low_event_count) %>% 
  count(plasmid, host) %>% 
  filter(n < 3)

2. Background subtraction

Check the background fluorescence levels across different days.

dat %>% filter(host == "156") %>% 
  pivot_longer(BL1.H:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = date, y = intensity)) + geom_point(aes(shape = well)) +
  stat_summary(fun.data = "mean_se", geom = "pointrange", color = "red", position = position_nudge(x = 0.1)) +
  facet_wrap(~parameter, scale = "free_y") +
  theme_bw(base_size = 14)

03/29/22 D9 showed higher background fluorescence than the others. after checking the gating result, I found that NA-156-2 sample in D9 had a slightly wider distribution, and the 2d cluster gate is also wider than the other two samples. The cell size normalized intensity doesn’t show the same pattern, suggesting that at least some of the difference could be due to cell size variation. Try just averaging the three blank wells and perform the background removal.

Subtract the background

bg <- dat %>% 
  filter(host == "156") %>% 
  group_by(date) %>% 
  summarize(across(FSC.H:nRFP, ~ round(mean(.x),1))) %>% 
  column_to_rownames(var = "date")

dat1 <- dat %>% 
  select(date:host, FSC.H:nRFP, low_event_count, mutation) %>%
  mutate(
    BL1.H = BL1.H - bg[date, "BL1.H"],
    YL2.H = YL2.H - bg[date, "YL2.H"],
    nGFP = nGFP - bg[date, "nGFP"],
    nRFP = nRFP - bg[date, "nRFP"],
  )

Double check that the subtraction worked correctly

dat1 %>% filter(host == "156") %>% 
  pivot_longer(BL1.H:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = date, y = intensity)) + geom_point(aes(shape = well)) +
  stat_summary(fun.data = "mean_se", geom = "pointrange", color = "red", position = position_nudge(x = 0.1)) +
  facet_wrap(~parameter, scale = "free_y") +
  theme_bw(base_size = 14)

# remove the yH156 samples
dat1 <- dat1 %>% 
  filter(host != "156") %>% 
  mutate(host = ordered(host, levels = c("555", "373"), labels = c("PHO2", "pho2∆")))

Export the background-subtracted values for later use with the shinyapp

write_tsv(dat1, "../output/20220614-mar22-batch-bg-sub-data.tsv")

3. Consistency across plates.

Check the background-subtracted fluorescence values for the host strains (yH373 and yH555) as well as the positive control strain (pH188 in yH373 or yH555 backgrounds), both of which are present on each plate.

dat1 %>% filter(plasmid %in% c("188", "NA")) %>% 
  pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = plasmid, y = intensity, color = date)) + 
  geom_point(position = position_dodge(0.5)) + 
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scale = "free_y") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(face = 3))

1. The host strain fluorescence levels appear to be consistent across the three days. Additionally, their GFP levels are zero, as expected.
2. The positive control strain has strong BL1.H and correspondingly has strong YL2.H.

The size normalized values show slightly lower variance:

dat1 %>% filter(plasmid %in% c("188", "NA")) %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = plasmid, y = intensity, color = date)) + 
  geom_point(position = position_dodge(0.5)) + 
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scale = "free_y") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(face = 3))

4. Well position

Rationale

• In my original plate design, I placed a positive control strain (envisioned CgPho4-mNeon in PHO2 background for example) in the control columns (1, 5, 9) on every other row (A, C, E, G). The reason for this is to use it to identify any well position effect on the fluorescence readings.
• In Emily’s implementation, she instead put a pair of strains, namely pH188-yH323 and pH188-yH555 in these wells. This means I cannot use the positive control wells exactly the way as I designed them. But I can still use them to spot any trend.

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  ggplot(aes(x = row, y = intensity, group = date)) + 
  geom_point(aes(color = date), size = 1) + 
  #geom_line(aes(color = date), size = 0.4) +
  stat_summary(aes(color = date), geom = "line", fun = mean) +
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

• no obivous trend between the rows (or columns, not shown)
• a clear correlation between YL2.H and BL1.H. at least part of this is due to the cell size differences – see size normalized data below:

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  pivot_longer(c(YL2.H, nRFP), names_to = "parameter", values_to = "intensity") %>% 
  ggplot(aes(x = row, y = intensity, group = date)) + 
  geom_point(aes(color = date), size = 1) + 
  #geom_line(aes(color = date), size = 0.4) +
  stat_summary(aes(color = date), geom = "line", fun = mean) +
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

• What we care about is the ratio between Pho4-mNeon and PHO5pr-mCherry

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  ggplot(aes(x = BL1.H, y = YL2.H)) + 
  geom_point(aes(shape = date, color = row), size = 3) + 
  stat_smooth(method = "lm", se = FALSE, color = "gray50", size = 0.5) +
  scale_color_viridis_d() +
  scale_shape_manual(values = 15:17) +
  #scale_x_log10() + scale_y_log10() +
  #facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

`geom_smooth()` using formula 'y ~ x'

• There is variation in BL1.H and correspondingly in YL2.H, but the ratio between the two are very consistent

Let’s check both ScPho4 (pH194) and CgPho4 (pH188) to see if their behaviors are consistent across the days.

dat1 %>% 
  filter(plasmid %in% c("188", "194")) %>% 
  mutate(`YL2.H/BL1.H` = YL2.H/BL1.H,
  #mutate(`YL2.H/BL1.H` = nRFP/nGFP,
         Pho4 = factor(plasmid, levels = c("194", "188"), 
                       labels = c("ScPho4", "CgPho4"))) %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  ggplot(aes(x = host, y = `YL2.H/BL1.H`)) + 
  geom_point(color = alpha("gray50", 0.9), position = position_jitter(0.2)) + 
  stat_summary(fun.data = "mean_se", geom = "pointrange", color = "red") +
  facet_grid(Pho4~date) +
  theme_bw(base_size = 14)

5. High variance samples

Summarize the background subtracted data by calculating the means and cv for each strain.

dat2 <- dat1 %>% 
  select(-nGFP, -nRFP) %>% 
  pivot_longer(FSC.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  group_by(date, plasmid, host, parameter) %>% 
  summarize(
    n = n(),
    mean = num(mean(intensity), digits = 0),
    cv = num(sd(intensity)/mean(intensity), digits = 2)
  ) %>% 
  arrange(desc(cv))

`summarise()` has grouped output by 'date', 'plasmid', 'host'. You can override using the `.groups` argument.

dat2 %>% 
  filter(plasmid != "NA", parameter != "FSC.H") %>% 
  ggplot(aes(x = cv)) + geom_histogram(aes(y = ..density../20)) + stat_ecdf() + 
  geom_hline(yintercept = 0.8, linetype = 2) +
  ylab("cumulative density") + 
  facet_wrap(~parameter) +
  theme_cowplot()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The histogram’s y-axis is not shown. The line graph represents the empirical CDF, and the dotted horizontal line is at 80%. GFP is more variable than RFP, likely because the absolute values of the former is lower. For both, ~80% of the samples have a CV < 20%.

Do the same for the cell size-normalized values

ndat2 <- dat1 %>% 
  select(-BL1.H, -YL2.H) %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  group_by(date, plasmid, host, parameter) %>% 
  summarize(
    n = n(),
    mean = num(mean(intensity), digits = 0),
    cv = num(sd(intensity)/mean(intensity), digits = 2)
  ) %>% 
  arrange(desc(cv))

`summarise()` has grouped output by 'date', 'plasmid', 'host'. You can override using the `.groups` argument.

ndat2 %>% 
  filter(plasmid != "NA", parameter != "FSC.H") %>% 
  ggplot(aes(x = cv)) + geom_histogram(aes(y = ..density../20)) + stat_ecdf() + 
  geom_hline(yintercept = 0.8, linetype = 2) +
  ylab("cumulative density") + 
  facet_wrap(~parameter) +
  theme_cowplot()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

CV is smaller on the cell size normalized fluorescence values.

Main

Plotting functions

Here I’d like to develop a series of plotting functions that take the names or any part of the Pho4 chimera annotation as input and plot their results in a variety of ways Load the Pho4 plasmid information and merge with the reshaped data

meta <- read_tsv("../data/20220330-chimera-Pho4-makeup.tsv", col_types = "cccc")

Plotting function

# extract ximera names
refs <- c("CCCCC","SSSSS")
ximeras <- setdiff(meta$symbol, refs)
# make a test set
test <- sample(ximeras, 6)
# subset data
tmp <- meta %>% 
  filter(symbol %in% c(refs,test)) %>% 
  inner_join(dat1, by = "plasmid") %>% 
  mutate(symbol = factor(symbol, levels = c(refs, test)))

Plot individual components

tmp %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("nRFP", "nGFP"))) %>% 
  #pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("YL2.H", "BL1.H"))) %>% 
  ggplot(aes(x = symbol, y = intensity, group = host)) + 
  geom_bar(aes(fill = host), width = 0.5, alpha = 0.8,
           stat = "summary", fun = "mean", position = position_dodge(0.5)) +
  geom_point(aes(color = host), size = 1,
             position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_color_manual(values = c("PHO2" = "gray20", "pho2∆" = "gray40")) +
  scale_fill_viridis_d(begin = 0.2, end = 0.6) +
  #stat_summary(fun = "mean", geom = "crossbar", color = "red", width = 0.25,
  #             position = position_dodge(0.75), ) +
  facet_wrap(~parameter, scale = "free_y", ncol = 1) +
  xlab("Pho4 chimera") + expand_limits(y = 0) +
  theme_gray(base_size = 14) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

tmp %>% 
  #mutate(`YL2.H/BL1.H` = YL2.H/BL1.H,
  mutate(`R/G` = nRFP/nGFP) %>% 
  ggplot(aes(x = symbol, y = `R/G`, group = host)) + 
  geom_bar(aes(fill = host), width = 0.5, alpha = 0.8,
           stat = "summary", fun = "mean", position = position_dodge(0.5)) +
  geom_point(aes(color = host), position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_color_manual(values = c("PHO2" = "gray20", "pho2∆" = "gray40")) +
  scale_fill_viridis_d(begin = 0.2, end = 0.6) +
  #stat_summary(fun = "mean", color = "red", geom = "crossbar", width = 0.2,
  #             position = position_dodge(0.75), ) +
  theme_gray(base_size = 14) + xlab("Pho4 chimera") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

Next idea is to plot the ratio between the mean R/G for PHO2 and mean R/G for pho2∆, over the mean R/G for PHO2

First summarize the data.

datsum <- dat1 %>%
  filter(plasmid != "NA") %>% 
  mutate(`R/G` = YL2.H/BL1.H, `nR/G` = nRFP/nGFP) %>% 
  group_by(date, plasmid, host) %>% 
  summarize(across(c(BL1.H, nGFP, YL2.H, nRFP, `R/G`, `nR/G`), mean), .groups = "drop") %>% 
  pivot_wider(names_from = host, values_from = BL1.H:`nR/G`) %>% 
  mutate(`pho2∆/PHO2` = `R/G_pho2∆`/`R/G_PHO2`,
         `n.pho2∆/PHO2` = `nR/G_pho2∆`/`nR/G_PHO2`)

# useful to set a flag for low activity mutants
low.act.th <- datsum %>% filter(plasmid == "194") %>% summarize(across(.cols = where(is.numeric), .fns = mean)) %>% unlist()

Then extract the subset for plotting

tmpsum <- meta %>% 
  filter(symbol %in% c(refs,test)) %>% 
  inner_join(datsum, by = "plasmid") %>% 
  mutate(symbol = factor(symbol, levels = c(refs, test)))

Design the plot

tmpsum %>% 
  mutate(Activity = ifelse(`R/G_PHO2` < 2*low.act.th["R/G_pho2∆"], "low", "pass")) %>% 
  ggplot(aes(x = symbol, y = `pho2∆/PHO2`)) + 
  geom_col(aes(group = date, fill = `R/G_PHO2`), width = 0.75, color = "gray50",
           position = position_dodge(0.9)) +
  #geom_point(aes(color = host), position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_fill_gradient2("Activity") +
  #scale_color_manual(values = c(alpha("black",0), "red3")) +
  #stat_summary(fun = "mean", color = "red", geom = "crossbar", width = 0.2,
  #             position = position_dodge(0.75), ) +
  facet_grid(.~Activity, scales = "free_x", space = "free_x", labeller = "label_both") +
  theme_bw(base_size = 14) + xlab("Pho4 chimera") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

X-Y plot

p3 <- tmpsum %>% 
  mutate(`nR/G_PHO2` = signif(`nR/G_PHO2`, digits = 2),
         `nR/G_pho2∆` = signif(`nR/G_pho2∆`, digits = 2)) %>% 
  ggplot(aes(x = `nR/G_PHO2`, y = `nR/G_pho2∆`, label = symbol)) + 
  geom_point(size = 2.5) + geom_abline(slope = 1) +
  theme_gray(base_size = 14)
ggplotly(p3, tooltip = c("label", "x", "y"))

Warning: `gather_()` was deprecated in tidyr 1.2.0.
Please use `gather()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.