CTA1 induction under phosphate starvation in S. cerevisiae

suppressPackageStartupMessages(library(tidyverse))

Warning message:
In if (charToRaw(x) < 20) paste("\\u", toupper(format(as.hexmode(as.integer(charToRaw(x))),  :
  the condition has length > 1 and only the first element will be used

suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(GEOquery))
suppressPackageStartupMessages(library(SummarizedExperiment))

Background

Jinye has discovered that the sole catalase gene CTA1 is induced under phosphate starvation in C. glabrata, and is possibly regulated by the general stress response TF CgMsn4. We want to determine if the same are true in S. cerevisiae, i.e. if either ScCTA1 or ScCTT1 is induced under phosphate starvation and whether ScMsn4 is responsible for its induction (although Jinye has also shown that ScMsn4 DOES NOT localize into the nucleus upon phosphate starvation, or at least not in sufficient amount for it to be detectable by fluorescent imaging).

To answer the first question, the idea is to reanalyze two transcriptome-profiling dataset, one from Zhou and O’Shea 2011 (PMID: 21700227) and the other from Gurvich et al. 2017 (PMID: 29236696).

Zhou and O’Shea 2011

Data

GEO#	Description	Reference
GSE23580	Microarray expression analyses of S. cerevisiae wt, pho4∆, pho2∆ in rich or no Pi media	Zhou & O’Shea 2011

Zhou X, O’Shea EK. 2011. Integrated approaches reveal determinants of genome-wide binding and function of the transcription factor Pho4. Mol Cell 42:826–836.

Analysis

• GEOquery workshop offered by Jason Ratcliff at the Iowa Institute of Human Genetics

getGEO() will return an ExpressionSet object, which is then converted into a “SummarizedExperiment” object, which is a more modern data structure that is easier to deal with.

Sys.setenv("VROOM_CONNECTION_SIZE" = 131072*10) # increase the local cache size
gse <- getGEO(GEO = "GSE23580", destdir = "data/Zhou2011")

Found 1 file(s)
GSE23580_series_matrix.txt.gz
Using locally cached version: data/Zhou2011/GSE23580_series_matrix.txt.gz

Rows: 6314 Columns: 67
── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr  (1): ID_REF
dbl (66): GSM578408, GSM578409, GSM578410, GSM578411, GSM578412, GSM578413, GSM578414, GSM578415, GSM578416, GSM578417, GSM578418, GSM57841...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Using locally cached version of GPL10786 found here:
data/Zhou2011/GPL10786.soft 

gse23580 <- as(gse$GSE23580_series_matrix.txt.gz, "SummarizedExperiment")

The experimental information is stored in the colData fields. The ones we are interested in are:

colData(gse23580) %>% 
  as_tibble() %>%
  filter(grepl("Wild type no vs high Pi conditions", title) | grepl("Comparison 1$", title)) %>% 
  select(title, geo_accession, organism = organism_ch1, strain = characteristics_ch1, 
         condition_ch1 = characteristics_ch1.2, label_ch1, condition_ch2 = characteristics_ch2.2, label_ch2)

The first three are part of the mutant cycle while the latter four are not. The latter four are said to have incorporated a dye swap, although I can’t tell how the swap was done from the table above.

Separately extract the two subsets and examine them separately.

gse1 <- gse23580[, grepl("Comparison 1$", gse23580$title)]
gse2 <- gse23580[, grepl("Wild type no vs high Pi conditions", gse23580$title)]

The probe annotation is stored as row metadata:

#gse23580 %>%  rowData() %>% as_tibble() %>% DT::datatable()

Check to see if the assay data is indeed normalized and log2 transformed, by plotting the values as boxplots:

cbind(assay(gse1, 1), assay(gse2, 1)) %>% boxplot()

Looks like they are.

Now we are ready to check CTA1 expression. First we need to find out the probe ID for CTA1, and then we can extract and plot the log2 transformed ratio. For comparison, we will include two other genes, one PHO target and one housekeeping.

Mutant Cycle

genes <- c("CTA1", "CTT1", "HSP12", "ALG9", "PHO8", "SPL2")
probes <- rowData(gse1) %>% 
  as_tibble() %>% 
  filter(Name %in% genes) %>% 
  select(ID, Name)

res1 <- assay(gse1)[probes$ID,] %>% as_tibble(rownames = "ID") %>% 
  left_join(probes, by = "ID") %>% 
  pivot_longer(cols = starts_with("GSM"), names_to = "Acc", values_to = "log2") %>% 
  left_join(colData(gse1) %>% as_tibble() %>% select(Sample = title, Acc = geo_accession), by = "Acc") %>% 
  mutate(Sample = gsub(" Comparison 1||Mutant Cycle ", "", Sample)) %>% 
  select(ID, Gene = Name, Acc, Sample, log2)

ggplot(res1, aes(x = Sample, y = log2, fill = Gene)) + 
  geom_col(position = position_dodge()) +
  ylab("log2 No Pi vs High Pi") + theme_cowplot()

ggsave("output/20210920-Zhou2011-mutant-cycle-data-plot.png")

Saving 7.29 x 4.51 in image

Dye swap

res2 <- assay(gse2)[probes$ID,] %>% as_tibble(rownames = "ID") %>% 
  left_join(probes, by = "ID") %>% 
  pivot_longer(cols = starts_with("GSM"), names_to = "Acc", values_to = "log2") %>% 
  left_join(colData(gse2) %>% as_tibble() %>% select(Sample = title, Acc = geo_accession), by = "Acc") %>% 
  mutate(Sample = gsub("Wild type no vs high Pi conditions ", "", Sample)) %>% 
  select(ID, Gene = Name, Acc, Sample, log2)

ggplot(res2, aes(x = Sample, y = log2, fill = Gene)) + geom_col(position = position_dodge()) +
  ylab("log2 No Pi vs High Pi") + theme_cowplot()

ggsave("output/20210920-Zhou2011-dye-swap-data-plot.png")

Saving 7.29 x 4.51 in image

Gurvich 2017

Data

Gurvich Y, Leshkowitz D, Barkai N. 2017. Dual role of starvation signaling in promoting growth and recovery. PLOS Biology 15:e2002039.

The data is extracted from the supplementary table 2, which records the data behind all the figures in the paper. The raw sequencing data were deposited in SRA, but not submitted to GEO. I decided to use their processed data, which is in log2 ratios (compared with the 0 min time point in most cases).

The extracted data is stored in an excel file with two sheets, corresponding to the 0 mM Pi and 0.06 mM Pi timecourse.

gurvich.file <- "data/Gurvich2017/20210921-extracted-supp-data-Fig1C.xlsx"
excel_sheets(gurvich.file)

[1] "Fig1c-0Pi"    "Fig3C-0.06Pi" "Fig1C-0.06Pi"

gd00 <- read_excel(gurvich.file, sheet = 1, na = "NaN") %>% 
  pivot_longer(!Group:Gene, names_to = "Time", values_to = "Log2.ratio") %>% 
  mutate(Time = as.numeric(Time))
gd06 <- read_excel(gurvich.file, sheet = 3, na = "NaN") %>% 
  pivot_longer(!Group:Gene, names_to = "Time", values_to = "Log2.ratio") %>% 
  mutate(Time = as.numeric(Time))

Analysis

I selected four genes in the “stress” category based on Hasan et al. 2002 Mol. Microbiol. “The control of the yeast H2O2 response by the Msn2/4 transcription factors”.

We first look at the 0 mM Pi data

genes <- c("CTT1", "TSL1", "HSP12", "HSP42", "PHO84", "SPL2", "PHO5")
p <- gd00 %>% 
  filter(Gene %in% genes, Time < 20) %>% 
  mutate(Fold = 2^Log2.ratio) %>% 
  ggplot(aes(x = Time, y = Fold, color = Gene)) + geom_point() + geom_line(size = 0.5) +
  xlab("Time (hours after -> 0 mM Pi)") + ylab("Fold induction relative to 0 min") +
  facet_wrap(~ Group)
p + theme_cowplot()

Notice how the stress genes show an early peak at the 2 hour time point, after which it drops to a lower level before reaching a second higher peak after quite some time (>8 hrs). What could be the explanation for this? The 2 hour time point induction shouldn’t be a result of switching media, which requires centrifugation or filtering, both of which could induce Msn2/4 translocation (?). If that’s not a likely reason, then we have to assume that these putative Msn2/4 targets are induced because the TFs responded to Pi starvation, likely at around 1 hour time point.

Next let’s look at the 0.06 mM Pi data. Note that because

p <- gd06 %>% 
  filter(Gene %in% genes, Time < 20) %>% 
  mutate(Fold = 2^Log2.ratio) %>% 
  ggplot(aes(x = Time, y = Fold, color = Gene)) + geom_point() + geom_line(size = 0.5, linetype = 2) +
  geom_hline(yintercept = 2, linetype = 3) +
  xlab("Time (hours after -> 0.06 mM Pi)") + ylab("Fold induction relative to 0 min") +
  facet_wrap(~ Group, nrow = 2)
p + theme_cowplot()

Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 1 row(s) containing missing values (geom_path).

The dotted horizontal line indicates 2 fold induction.

What’s different in this plot from the previous one is that at the intermediate level of Pi, the cells exhibit a strongly biphasic induction of the PHO genes, with ~6.5 hours as the separation point. Before 6.5 hours, PHO5 and PHO84 are mildly induced and after that time point, both are more strongly induced. SPL2 is an early gene that gets induced to a higher level even in the early time point (it is also less dependent on Pho2). For the stress genes, they really aren’t induced until after 6.5 hrs when Pi starts to become growth limiting. This suggests that the signaling to Msn2/4 is probably via the sensing of nutrient levels (PKA? TOR)