figure_5_supp_5.Rmd

---
title: "Figure 5 and Supp Fig 5 analyses"
author: "Emily Butka"
date: "2024-06-13"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(Seurat)
library(ggplot2)
library(plyr)
```

## Load data

https://figshare.com/s/396c5f4ebd3c23d20094

```{r load-data}
p21 <- readRDS("~/Downloads/p21_processed.rds")
```

```{r load-data}
p32 <- readRDS("~/Downloads/p32_naive_processed.rds")
```

```{r}
merrick <- readRDS("~/Documents/GitHub/merrick_etal_processed.rds")
```

```{r}
sox_oe <- readRDS("~/Documents/GitHub/p32_sox9_oe_progenitors_processed.rds")
```

```{r}
oe_raw <- read.table("~/OneDrive/MorrisLab/Guillermo_project/Sox9_OE_quant_FACS_2024.csv")$V1
shrna_raw <- read.table("~/OneDrive/MorrisLab/Guillermo_project/Sox9_shrna_quant_facs_2024.csv")$V1
```

```{r}
preadips <- readRDS("~/Downloads/p32_sox9_oe_preadipocytes_processed.rds")
```


## Custom functions

```{r}
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}

summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
```

```{r}
summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
```

## Visualizations

```{r}
p21$new_classification_short <- factor(p21$new_classification_short, levels = c("Int prog", "Trans int prog", "Imm preadip", "Comm preadip", "Multi ID", "Unknown"))
```

```{r}
s5al <- VlnPlot(p21, features = "Dlk1", group.by = "new_classification_short", cols = colpal6[1:4], pt.size = 0) + NoLegend() + theme(axis.title.x = element_blank(), plot.title = element_blank())
s5al
```

```{r}
p32$capy_new_class_short <- factor(p32$capy_new_class_short, levels = c("Prog", "Trans prog", "Imm preadip", "Comm preadip", "Multi ID", "Unknown"))
```

```{r}
s5am <- VlnPlot(p32, features = "Dlk1", group.by = "capy_new_class_short", cols = colpal6, pt.size = 0) + NoLegend() + theme(axis.title.x = element_blank(), plot.title = element_blank())
s5am
```

```{r}
merrick$merrick_class <- gsub("adipocyte", "Adipocyte", merrick$merrick_class)
merrick$merrick_class <- gsub("committed_preAdipocyte", "Committed preadipocyte", merrick$merrick_class)
merrick$merrick_class <- gsub("endothelial", "Endothelial", merrick$merrick_class)
merrick$merrick_class <- gsub("group_3_cd142", "Group 3", merrick$merrick_class)
merrick$merrick_class <- gsub("group_4", "Group 4", merrick$merrick_class)
merrick$merrick_class <- gsub("group_5", "Group 5", merrick$merrick_class)
merrick$merrick_class <- gsub("group_6", "Group 6", merrick$merrick_class)
merrick$merrick_class <- gsub("interstitial_progenitor", "Interstitial progenitor", merrick$merrick_class)
merrick$merrick_class <- gsub("Multi_ID", "Multi ID", merrick$merrick_class)
merrick$merrick_class <- gsub("neural_crest", "Neural crest", merrick$merrick_class)
merrick$merrick_class <- gsub("smooth_muscle", "Smooth muscle", merrick$merrick_class)
```

```{r}
s5ar <- VlnPlot(merrick, features = "Dlk1", group.by = "merrick_class", cols = c("#FFC33B", "#FF6E3A", "#E20134", "#9F0162", "#FFB2FD", "#00C2F9", "#008DF9", "#8400CD", "#FF5AAF", "#009F81", "gray"), pt.size = 0) + NoLegend() + theme(axis.title.x = element_blank(), plot.title = element_blank())
s5ar
```

## EdU proliferation assay 

### Recovered cell numbers

```{r}
facs_counts <- data.frame(experiment = c(rep("KO", 11), rep("OE", 6)), sample = c(rep("control", 6), rep("KO", 5), rep("control", 3), rep("OE", 3)), counts = c(567, 741, 190, 255, 385, 325, 99, 183, 117, 252, 270, 3335, 3239, 4393, 14368, 8144, 15078))
```

```{r}
facs_counts_sum <- ddply(facs_counts, c("sample", "experiment"), .fun = summary_func, "counts")
```

```{r}
facs_counts$mean <- NA
facs_counts$sd <- NA
for(i in 1:nrow(facs_counts)) {
  facs_counts$mean[i] <- facs_counts_sum$mean[which(facs_counts_sum$sample == facs_counts$sample[i] & facs_counts_sum$experiment == facs_counts$experiment[i])]
  facs_counts$sd[i] <- facs_counts_sum$sd[which(facs_counts_sum$sample == facs_counts$sample[i] & facs_counts_sum$experiment == facs_counts$experiment[i])]
}
```

```{r}
pbl <- ggplot(facs_counts[which(facs_counts$experiment == "KO"), ], aes(x = sample, y = counts)) + stat_summary(fun.y=mean,position=position_dodge(width=0.95),geom="bar") + geom_jitter(width = 0.1) +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd),
                width=.2, position=position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Number of recovered cells")
p5br <- ggplot(facs_counts[which(facs_counts$experiment == "OE"), ], aes(x = sample, y = counts)) + stat_summary(fun.y=mean,position=position_dodge(width=0.95),geom="bar") + geom_jitter(width = 0.1) +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd),
                width=.2, position=position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Number of recovered cells")
p5bl
p5br
```

```{r}
t.test(facs_counts$counts[which(facs_counts$experiment == "KO" & facs_counts$sample == "KO")], facs_counts$counts[which(facs_counts$experiment == "KO" & facs_counts$sample == "control")])
```

```{r}
t.test(facs_counts$counts[which(facs_counts$experiment == "OE" & facs_counts$sample == "OE")], facs_counts$counts[which(facs_counts$experiment == "OE" & facs_counts$sample == "control")])
```

### EdU absorption

```{r}
preadip_edu <- data.frame(sample = rep(c("preadip_OE_control", "preadip_OE"), each = 3), edu_absorb = c(74.5, 72.6, 74.6, 66.8, 68.3, 66))
```

```{r}
prof_data <- data.frame(sample = c(rep("prog_control", 3), rep("prog_KO", 3), rep("prog_OE", 3)), edu_absorb = c(43.7, 44.4, 41.5, 12.0, 11.0, 6.08, 48.2, 51.9, 49.3))
```

```{r}
preadip_edu2 <- data.frame(sample = rep(c("preadip_KO_control", "preadip_KO"), each = 3), edu_absorb = c(52.8, 52.4, 54.9, 1.05, 0.65, 0.85))
```

```{r}
edu_all <- rbind(prof_data, preadip_edu2, preadip_edu)
```

```{r}
library(plyr)
edu_all_sum <- ddply(edu_all, "sample", .fun = summary_func, "edu_absorb")
```

```{r}
edu_all$mean <- NA
edu_all$sd <- NA
for(i in 1:nrow(edu_all)) {
  edu_all$mean[i] <- edu_all_sum$mean[which(edu_all_sum$sample == edu_all$sample[i])]
  edu_all$sd[i] <- edu_all_sum$sd[which(edu_all_sum$sample == edu_all$sample[i])]
}
```

```{r}
edu_all$mean_normalized <- NA
edu_all$sd_normalized <- NA
for(i in 1:nrow(edu_all)) {
  if(grepl("control", edu_all$sample[i])) {
    edu_all$mean_normalized[i] <- 1
  } else if(grepl("prog", edu_all$sample[i])) {
    edu_all$mean_normalized[i] <- edu_all$mean[i] / edu_all$mean[which(edu_all$sample == "prog_control")][1]
  } else if(grepl("preadip_KO", edu_all$sample[i])) {
    edu_all$mean_normalized[i] <- edu_all$mean[i] / edu_all$mean[which(edu_all$sample == "preadip_KO_control")][1]
  } else if(grepl("preadip_OE", edu_all$sample[i])) {
    edu_all$mean_normalized[i] <- edu_all$mean[i] / edu_all$mean[which(edu_all$sample == "preadip_OE_control")][1]
  }
}

```

```{r}
edu_all$mean_multiplier <- edu_all$mean_normalized / edu_all$mean
```

```{r}
edu_all$sd_normalized <- edu_all$sd * edu_all$mean_multiplier
edu_all$value_normalized <- edu_all$edu_absorb * edu_all$mean_multiplier
```

```{r}
edu_all$sample_condensed <- ifelse(grepl("control", edu_all$sample), "Control", edu_all$sample)
```

```{r}
edu_all$sample_condensed <- factor(edu_all$sample_condensed, levels = c("Control", "prog_KO", "prog_OE", "preadip_KO", "preadip_OE"))
```

```{r}
edu_all$sd_normalized_condensed <- NA
for(i in 1:nrow(edu_all)) {
  edu_all$sd_normalized_condensed[i] <- sd(edu_all$value_normalized[which(edu_all$sample_condensed == edu_all$sample_condensed[i])])
}
```

```{r}
p5c <- ggplot(edu_all, aes(x = sample_condensed, y = value_normalized)) + geom_hline(yintercept = 1, linetype = "dashed") +
  stat_summary(fun.y = mean, position = position_dodge(width = 0.95), geom="bar") + geom_jitter(width = 0.1) +
  geom_errorbar(aes(ymin = mean_normalized - sd_normalized_condensed, ymax = mean_normalized + sd_normalized_condensed),
                width = .2, position = position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Fold change in EdU absorption relative to matched control")
p5c
```

```{r}
t.test(edu_all$edu_absorb[which(edu_all$sample == "prog_control")], edu_all$edu_absorb[which(edu_all$sample == "prog_KO")])
```

```{r}
t.test(edu_all$edu_absorb[which(edu_all$sample == "prog_control")], edu_all$edu_absorb[which(edu_all$sample == "prog_OE")])
```

```{r}
t.test(edu_all$edu_absorb[which(edu_all$sample == "preadip_KO_control")], edu_all$edu_absorb[which(edu_all$sample == "preadip_KO")])
```

```{r}
t.test(edu_all$edu_absorb[which(edu_all$sample == "preadip_OE_control")], edu_all$edu_absorb[which(edu_all$sample == "preadip_OE")])
```

## S5B????

## Apoptosis Annexin V assay

```{r}
apop_data <- data.frame(sample = rep(c("prog_control", "prog_KO", "preadip_control", "preadip_KO"), each = 3), annexin = c(9.23, 8.338, 12.2, 27.3, 19.4, 19.5, 10.2, 10.5, 13, 26.2, 28.2, 17.8))
```

```{r}
apop_sum <- ddply(apop_data, "sample", .fun = summary_func, "annexin")
```

```{r}
apop_data$mean <- NA
apop_data$sd <- NA
for(i in 1:nrow(apop_data)) {
  apop_data$mean[i] <- apop_sum$mean[which(apop_sum$sample == apop_data$sample[i])]
  apop_data$sd[i] <- apop_sum$sd[which(apop_sum$sample == apop_data$sample[i])]
}
```

```{r}
apop_data$sample <- factor(apop_data$sample, levels = c("prog_control", "prog_KO", "preadip_control", "preadip_KO"))
```

```{r}
s5d <- ggplot(apop_data, aes(x = sample, y = annexin)) + stat_summary(fun.y = mean, position = position_dodge(width = 0.95), geom="bar") + geom_jitter(width = 0.1) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd),
                width = .2, position = position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Annexin absorption (%)") 
s5d
```

```{r}
t.test(apop_data$annexin[which(apop_data$sample == "prog_control")],
       apop_data$annexin[which(apop_data$sample == "prog_KO")])
```

```{r}
t.test(apop_data$annexin[which(apop_data$sample == "preadip_control")],
       apop_data$annexin[which(apop_data$sample == "preadip_KO")])
```

## QPCR data

```{r}
qpcr <- data.frame(apc = rep(c("Prog", "Preadip"), each = 24), 
                   gene = rep(rep(c("Sox9", "AdipoQ", "Plin1", "Pparg"), each = 6), 2), 
                   sample = rep(rep(c("Control", "shRNA"), each = 3), 8), 
                   value = c(1, 1.268228, 1.21750451, # Prog Sox9 control
                             0.0229826, 0.2188674, 0.0824412, # Prog Sox9 shRNA
                             1, 1.01085588, 4.28820553, # Prog AdipoQ control
                             NA, NA, 0.00245533, # Prog AdipoQ shRNA
                             1, 0.97276707, 2.65655492, # Prog Plin1 control
                             NA, NA, NA, # Prog Plin1 shRNA
                             1, 0.83766714, 2.61801842, # Prog Pparg control
                             NA, NA, NA, # Prog Pparg shRNA
                             1, 0.7919938, 0.6604467, # Preadip Sox9 control
                             0.1100334, 0.0818473, 0.1176783, # Preadip Sox9 shRNA
                             1, 0.69614013, 0.60811261, # Preadip AdipoQ control
                             0.00042786, 0.00167799, 0.00126341, # Preadip AdipoQ shRNA
                             1, 0.6749031, 0.6907908, # Preadip Plin1 control
                             0.0005587, 0.0014766, 0.0009472, # Preadip Plin1 shRNA
                             1, 0.709537, 0.7563131, # Preadip Pparg control
                             0.0008457, 0.0018834, 0.000727)) # Preadip Pparg shRNA
```

```{r}
qpcr$sd <- NA
qpcr$mean <- NA
qpcr_summary <- ddply(qpcr, c("apc", "gene", "sample"), .fun = summary_func, "value")
colnames(qpcr_summary)[4] <- "mean"
qpcr_summary$value <- NA
qpcr_plot <- rbind(qpcr, qpcr_summary)
```

```{r}
qpcr_plot$apc <- factor(qpcr_plot$apc, levels = c("Prog", "Preadip"))
qpcr_plot$gene <- factor(qpcr_plot$gene, levels = c("Sox9", "AdipoQ", "Plin1", "Pparg"))
```

```{r}
p5d <- ggplot(qpcr_plot, aes(x = apc, y = value, fill = sample)) + stat_summary(fun.y = mean, position = position_dodge(width=0.95), geom="bar") + geom_jitter(width = 0.1) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd),
                width=.2, position=position_dodge(.9)) +
  facet_wrap(vars(gene), nrow = 1) +
  theme_classic() +
  labs(x = "", y = "Fold change") +
  theme(strip.background = element_blank(), strip.placement = "outside", legend.position = "bottom") + 
  scale_fill_manual(values = c("#4383C4", "#E11D36"))
p5d
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Sox9" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Sox9" & qpcr$sample == "shRNA")])
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Sox9" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Sox9" & qpcr$sample == "shRNA")])
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "AdipoQ" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "AdipoQ" & qpcr$sample == "shRNA")],
       alternative = "two.sided")
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "AdipoQ" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "AdipoQ" & qpcr$sample == "shRNA")])
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Plin1" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Plin1" & qpcr$sample == "shRNA")],
       alternative = "two.sided")
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Plin1" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Plin1" & qpcr$sample == "shRNA")])
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Pparg" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Prog" & qpcr$gene == "Pparg" & qpcr$sample == "shRNA")],
       alternative = "two.sided")
```

```{r}
t.test(qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Pparg" & qpcr$sample == "Control")],
       qpcr$value[which(qpcr$apc == "Preadip" & qpcr$gene == "Pparg" & qpcr$sample == "shRNA")])
```

## Sox9 OE progenitors

```{r}
sox_oe_df <- data.frame(table(sox_oe$replicate, sox_oe$celltag, sox_oe$capy_new_class))
colnames(sox_oe_df) <- c("replicate", "celltag", "celltype", "Freq")
```

```{r}
sox_oe_df$rel_abundance <- NA
for(i in 1:nrow(sox_oe_df)) {
  sox_oe_df$rel_abundance[i] <- sox_oe_df$Freq[i] / sum(sox_oe_df$Freq[which(sox_oe_df$replicate == sox_oe_df$replicate[i] & sox_oe_df$celltag == sox_oe_df$celltag[i])])
}
```

```{r}
sox_oe_df$FC_rel_host <- NA
for(i in 1:nrow(sox_oe_df)) {
  sox_oe_df$FC_rel_host[i] <- sox_oe_df$rel_abundance[i] / sox_oe_df$rel_abundance[which(sox_oe_df$replicate == sox_oe_df$replicate[i] & sox_oe_df$celltype == sox_oe_df$celltype[i] & sox_oe_df$celltag == "None")]
}
```

```{r}
ctcapy_control_sum <- ddply(sox_oe_df[which(sox_oe_df$celltype != "Unknown"), ], c("celltype", "celltag"), .fun = summary_func, "rel_abundance")
```

```{r}
ctcapy_control_sum$celltype <- factor(ctcapy_control_sum$celltype, levels = c("progenitor", "transitioning_progenitor", "immature_preadipocyte", "committed_preadipocyte", "Multi_ID"))
```

```{r}
#ctcapy <- data.frame(table(p32$celltag, p32$capy_new_class_short))
p5fl <- ggplot(ctcapy_control_sum[ctcapy_control_sum$celltag == "None", ], aes(x="", y=mean, fill=celltype)) + geom_bar(stat="identity", width=1, color="white") + theme_void() + ggtitle("Host cells") + theme(legend.position="none", text = element_text(size = 8), plot.title = element_text(hjust = 0.5)) + labs(fill = "Cell type") + scale_fill_manual(name = "Cell type", values = colpal6)
p5fr <- ggplot(ctcapy_control_sum[ctcapy_control_sum$celltag == "OE", ], aes(x="", y=mean, fill=celltype)) + geom_bar(stat="identity", width=1, color="white") + theme_void() + ggtitle("OE CellTagged cells") + scale_fill_manual(name = "Cell type", values = colpal6) + theme(text = element_text(size = 8), legend.position = "none", plot.title = element_text(hjust = 0.5))
p5f <- p5fl + p5fr
p5f
```

```{r}
sox_oe_df <- data.frame(table(sox_oe$replicate, sox_oe$celltag, sox_oe$capy_new_class))
colnames(sox_oe_df) <- c("replicate", "celltag", "celltype", "Freq")
```

```{r}
sox_oe_df$rel_abundance <- NA
for(i in 1:nrow(sox_oe_df)) {
  sox_oe_df$rel_abundance[i] <- sox_oe_df$Freq[i] / sum(sox_oe_df$Freq[which(sox_oe_df$replicate == sox_oe_df$replicate[i] & sox_oe_df$celltag == sox_oe_df$celltag[i])])
}
```

```{r}
sox_oe_df$FC_rel_host <- NA
for(i in 1:nrow(sox_oe_df)) {
  sox_oe_df$FC_rel_host[i] <- sox_oe_df$rel_abundance[i] / sox_oe_df$rel_abundance[which(sox_oe_df$replicate == sox_oe_df$replicate[i] & sox_oe_df$celltype == sox_oe_df$celltype[i] & sox_oe_df$celltag == "None")]
}
```

```{r}
ctcapy_control_sum <- ddply(sox_oe_df[which(sox_oe_df$celltype != "Unknown"), ], c("celltype", "celltag"), .fun = summary_func, "FC_rel_host")
```

```{r}
ctcapy_control_sum$celltype <- factor(ctcapy_control_sum$celltype, levels = c("progenitor", "transitioning_progenitor", "immature_preadipocyte", "committed_preadipocyte", "Multi_ID"))
```

Compare to previous P32 data:

```{r}
adi_df <- data.frame(table(p32$dataset, p32$celltag, p32$capy_new_class))
colnames(adi_df) <- c("replicate", "celltag", "celltype", "Freq")
```

```{r}
adi_df$rel_abundance <- NA
for(i in 1:nrow(adi_df)) {
  adi_df$rel_abundance[i] <- adi_df$Freq[i] / sum(adi_df$Freq[which(adi_df$replicate == adi_df$replicate[i] & adi_df$celltag == adi_df$celltag[i])])
}
```

```{r}
adi_df$FC_rel_host <- NA
for(i in 1:nrow(adi_df)) {
  adi_df$FC_rel_host[i] <- adi_df$rel_abundance[i] / adi_df$rel_abundance[which(adi_df$replicate == adi_df$replicate[i] & adi_df$celltype == adi_df$celltype[i] & adi_df$celltag == "none")]
}
```

```{r}
adi_sum <- ddply(adi_df[which(adi_df$celltype != "Unknown"), ], c("celltype", "celltag"), .fun = summary_func, "FC_rel_host")
```

```{r}
adi_sum$celltype <- factor(adi_sum$celltype, levels = c("progenitor", "transitioning_progenitor", "immature_preadipocyte", "committed_preadipocyte", "Multi_ID"))
```

```{r}
adi_sum$celltag <- factor(adi_sum$celltag, levels = c("v2", "v1", "none"))
```

```{r}
prog_interest <- rbind(ctcapy_control_sum, adi_sum)
```

```{r}
p5g <- ggplot(prog_interest[which(prog_interest$celltag %in% c("v2", "OE")), ], aes(x = celltag, group = celltype, fill = celltype, y = mean)) + geom_hline(aes(yintercept = 1), linetype = "dashed") + geom_bar(stat="identity", color="black", width = 0.75, position = position_dodge(width = 0.9)) + geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width=.2, position=position_dodge(.9)) + theme_classic() + theme(legend.position = "none") + xlab("CellTagged populations") + ylab("Fold change relative to corresponding host proportions") + scale_fill_manual(values = colpal6)
p5g
```

## FACS data

```{r}
oe <- data.frame(t(data.frame(sapply(oe_raw, strsplit, ","))))
rownames(oe) <- NULL
colnames(oe) <- oe[1, ]
oe <- oe[-1, ]
oe$Mean <- as.double(oe$Mean)
oe$Stdev <- as.double(oe$Stdev)

shrna <- data.frame(t(data.frame(sapply(shrna_raw, strsplit, ","))))
rownames(shrna) <- NULL
colnames(shrna) <- shrna[1, ]
shrna <- shrna[-1, ]
shrna$Mean <- as.double(shrna$Mean)
shrna$Stdev <- as.double(shrna$Stdev)
```

```{r}
s5g <- ggplot(shrna, aes(x=Treatment, y=Mean)) +
  geom_bar(stat="identity", color="black", fill = "#1f77b4", position=position_dodge()) +
  geom_errorbar(aes(ymin=Mean-Stdev, ymax=Mean+Stdev),
                width=.2, position=position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Percent total APCs") + 
  facet_wrap(~Cell_Type, nrow = 1)
s5h <- ggplot(oe, aes(x=Treatment, y=Mean)) +
  geom_bar(stat="identity", color="black", fill = "#1f77b4", position=position_dodge()) +
  geom_errorbar(aes(ymin=Mean-Stdev, ymax=Mean+Stdev),
                width=.2, position=position_dodge(.9)) +
  theme_classic() +
  labs(x = "", y = "Percent total APCs") + 
  facet_wrap(~Cell_Type, nrow = 1)
s5g
s5h
```

## Sox9 OE Preadipocytes

```{r}
ctcapysox23 <- data.frame(table(preadips$celltag, preadips$capy_new_class_short))
ctcapysox23$Var2 <- gsub("_", " ", ctcapysox23$Var2)
ctcapysox23$Var2 <- factor(ctcapysox23$Var2, levels = c("Prog", "Trans prog", "Imm preadip", "Comm preadip", "Multi ID", "Unassigned"))
s5jl <- ggplot(ctcapysox23[ctcapysox23$Var1 == "Ctrl", ], aes(x="", y=Freq, fill=Var2)) + geom_bar(stat="identity", position = "fill") + ggtitle("Control Preadipocytes\n(n = 365)") + theme(legend.position="bottom") + labs(fill = "Cell type") + scale_fill_manual(values = colpal6) + theme_void()
s5jm <- ggplot(ctcapysox23[ctcapysox23$Var1 == "Sox9_OE", ], aes(x="", y=Freq, fill=Var2)) + geom_bar(stat="identity", position = "fill") + ggtitle("Sox OE Preadipocytes\n(n = 538)") + scale_fill_manual(values = colpal6) + theme_void()
s5jr <- ggplot(ctcapysox23[ctcapysox23$Var1 == "Dpp4", ], aes(x="", y=Freq, fill=Var2)) + geom_bar(stat="identity", position = "fill") + ggtitle("Control Progenitors\n(n = 266)") + scale_fill_manual(values = colpal6) + theme_void()
s5j <- grid.arrange(p5jl + theme(legend.position="none"), p5jm + theme(legend.position="none"), p5jr + theme(legend.position="bottom"), nrow=1)
s5j
```

```{r}
preadips$celltag[which(is.na(preadips$celltag))] <- "None"
```

```{r}
preadips_df <- data.frame(table(preadips$celltag, preadips$capy_new_class))
colnames(preadips_df) <- c("celltag", "celltype", "Freq")
```

```{r}
preadips_df$rel_abundance <- NA
for(i in 1:nrow(preadips_df)) {
  preadips_df$rel_abundance[i] <- preadips_df$Freq[i] / sum(preadips_df$Freq[which(preadips_df$celltag == preadips_df$celltag[i])])
}
```

```{r}
preadips_df$FC_rel_host <- NA
for(i in 1:nrow(preadips_df)) {
  preadips_df$FC_rel_host[i] <- preadips_df$rel_abundance[i] / preadips_df$rel_abundance[which(preadips_df$celltype == preadips_df$celltype[i] & preadips_df$celltag == "None")]
}
```

```{r}
ctcapy_control_sum_2 <- ddply(preadips_df[which(preadips_df$celltype != "Unknown"), ], c("celltype", "celltag"), .fun = summary_func, "FC_rel_host")
```

```{r}
ctcapy_control_sum_2$celltype <- gsub("interstitial_", "", ctcapy_control_sum_2$celltype)
ctcapy_control_sum_2$celltype <- factor(ctcapy_control_sum_2$celltype, levels = c("progenitor", "transitioning_progenitor", "immature_preadipocyte", "committed_preadipocyte", "Multi_ID"))
```

Compare to previous P32:

```{r}
preadip_interest <- rbind(ctcapy_control_sum_2, adi_sum)
```

```{r}
s5k <- ggplot(preadip_interest[which(preadip_interest$celltag %in% c("v1", "Sox9_OE")), ], aes(x = celltag, group = celltype, fill = celltype, y = mean)) + geom_hline(aes(yintercept = 1), linetype = "dashed") + geom_bar(stat="identity", color="black", width = 0.75, position = position_dodge(width = 0.9)) + geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width=.2, position=position_dodge(.9)) + theme_classic() + theme(legend.position = "bottom") + xlab("CellTagged populations") + ylab("Fold change relative to corresponding host proportions") + scale_fill_manual(values = colpal6) #+ scale_y_cut(breaks=3, which=1, scales=0.1) 
s5k
```

Stats:

```{r}
# test Sox9 OE against host, preadipocytes
randomization_test(preadips@meta.data, "capy_new_class", c("Sox9_OE", "None"), "celltag", "SeuratProject", "orig.ident", 10000)
```

## Analysis of GO terms generated from list of features co-enriched in P32 preadipocytes (scRNA) and Sox9 OE progenitors (bulk)

Analysis performed with Panther

```{r}
go_res_whole <- read.csv("~/OneDrive/MorrisLab/Guillermo_project/23_11_21_bulk/go_analysis/analysis.txt")
go_res_c <- go_res_whole[-(1:5), ]
```

```{r}
go_res <- data.frame(t(data.frame(strsplit(go_res_c, "\t"))))
rownames(go_res) <- NULL
colnames(go_res) <- go_res[1, ]
go_res <- tail(go_res, -1)
rownames(go_res) <- NULL
```

```{r}
tissue_or_celltype_specific_rows <- c(4, 5, 7, 9, 11, 14, 18, 19, 20, 27, 31, 35, 51, 54, 58, 59, 64, 65, 68, 70, 71, 72, 73, 76, 77, 80, 84, 85, 89, 93, 99, 100, 102, 107, 108, 110, 111, 113, 116, 117, 118, 119, 121, 125, 129, 133, 134, 136, 138, 140, 142, 150, 152, 153, 157, 158, 159, 163, 167, 180, 182, 184, 187, 189, 193, 210, 212, 219, 220, 222, 225, 230, 233, 234, 241, 246, 248, 250, 254, 255, 258, 262, 267, 268, 296, 297, 328, 331, 333, 341, 353, 356)
length(tissue_or_celltype_specific_rows)
```

```{r}
go_res_sub <- go_res[-tissue_or_celltype_specific_rows, ]
```

```{r}
go_categories <- data.frame(term = go_res_sub$`GO biological process complete`)
```

```{r}
go_categories$Proliferation <- ifelse(grepl("mitosis|spindle|mitotic|cell cycle|proliferation|replication|chromatid|kinetochore|division|chromosome separation|metaphase", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Signaling <- ifelse(grepl("signal|transport|import|export|communication", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Development <- ifelse(grepl("differentiation|development|morphogenesis", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Growth <- ifelse(grepl("growth", go_res_sub$`GO biological process complete`), 1, 0)
```

```{r}
go_categories$Metabolism <- ifelse(grepl("metabol", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Death <- ifelse(grepl("apopto|death", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Phosphorylation <- ifelse(grepl("phosphorylation", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Adhesion <- ifelse(grepl("adhesion", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Migration <- ifelse(grepl("migration|locomotion|motility", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Organization <- ifelse(grepl("organization", go_res_sub$`GO biological process complete`), 1, 0)
go_categories$Immune <- ifelse(grepl("immune|defense|viral|host|lymphocyte|presentation|antigen|killer|leukocyte|immunity|interleukin|inflammatory|chemotaxis|class|cytotoxicity|cytokine", go_categories$term), 1, 0)
```

```{r}
go_categories_plot <- data.frame(category = colnames(go_categories[, -1]), num_terms = apply(go_categories[, -1], 2, sum))
go_categories_plot <- go_categories_plot[order(go_categories_plot$num_terms), ]
go_categories_plot$category <- factor(go_categories_plot$category, levels = go_categories_plot$category)
p5h <- ggplot(go_categories_plot, aes(x = num_terms, y = category)) + geom_col(color = "black", fill = "#1f77b4") + theme_classic() + xlab("Number of GO terms significantly co-enriched in preadipocytes and Sox9-OE progenitors") + ylab("GO term category")
p5h
```