OSHetero2022_Supplemental_Figures_ReviewerRequests.Rmd

---
title: "Supplemental Figures - Reviewer Requests"
author: "Emily Franz, Matthew Cannon, Ryan Roberts"
date: "`r format(Sys.time(), '%m/%d/%Y')`"
output: 
  html_document:
    toc: true
    toc_float: true
    toc_depth: 5
    number_sections: false
    code_folding: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  tidy = TRUE,
  echo = TRUE,
  cache = TRUE,
  collapse = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  cache.lazy = FALSE,
  fig.show='hide'
)
```

```{r lib, cache=FALSE}
#source("scSeurat.R")
source("Downstream.v2.R")
set.seed(888)
# loading libraries
library(Seurat)
library(future)
library(ggplot2)
library(msigdbr)
library(clusterProfiler)
library(SingleR)
library(dplyr)
library(pheatmap)
library(RColorBrewer)
library(viridis) # inferno color palette
library(grid)
library(ggsci)
library(data.table)
library(ggvenn)
library(tibble)
library(tidyr)
library(tidyverse)
library(stringr)
library(reshape2)
library(pzfx)
library(ggridges)
library(effsize)
library(perm)
```

# Reviewer Requests

## S3: Add downregulated pathways

*6. Figure 1C, 3D, 5F— please include the pathway enrichment analysis on the downregulated genes, or please clarify the reason for only focusing on the upregulated genes for pathway enrichment analysis.*

### From Fig 1C

```{r 1c3d5fdown}
source("Downstream.v3.R")
# Prep directory
#Save as PDF
if(!dir.exists("Data/Supplement/")) {
  dir.create("Data/Supplement/")
}
if(!dir.exists("Data/Supplement/requests")) {
  dir.create("Data/Supplement/requests")
}

#### 1C ####
load("Data/os17_cx_CCR.RData")
load(file ="Data/os7_cx_CCR.RData")
# Pathway enrichment analysis
# Display adjust p-values
B.list <- list(OS17 = os17_cx,
               NCHOS7 = os7_cx)

em.hm.list <- list()

for (i in 1:length(B.list)) {
  em.hm.list[[i]] <- DGEA(B.list[[i]],
                          direction = "down")
}

for(i in 1:(length(em.hm.list))) {
  em.hm.list[[i]] <- setDT(em.hm.list[[i]],
                           keep.rownames = TRUE)[]
}

temp1 <- em.hm.list[[1]]
temp2 <- em.hm.list[[2]]
cx.pd <- full_join(temp1,
                   temp2,
                   by = "rn")

cx.pd <- column_to_rownames(cx.pd,
                            var = "rn")
cx.pd[is.na(cx.pd)] <- 1

#remove "HALLMARK_"
c <- rownames(cx.pd)
c <- gsub("HALLMARK_",
          "",
          c)

#remove "_" by removing special characters
c <- gsub("_",
          " ",
          c)
rownames(cx.pd) <- c

# Take the -log10
cx.pd.log <- -log10(cx.pd)
colnames(cx.pd.log) <- c("A0",
                         "A1",
                         "A2",
                         "A3",
                         "B0",
                         "B1",
                         "B2",
                         "B3")

# To prevent values that show up as -0.00, we take the absolute values
cx.pd.log <- abs(cx.pd.log)

figure_1c_down <- cx.pd.log %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column(var = "Sample") %>%
  pivot_longer(c(-Sample),
               names_to = "Pathway",
               values_to = "pval",
               names_repair = "minimal") %>%
  mutate(Signif = pval >= (-1 * log10(0.05))) %>%
  ggplot(.,
         aes(x = Sample,
             y = Pathway,
             fill = Signif)) +
  geom_tile() +
  geom_text(aes(label = sprintf("%0.2f",
                                pval))) + 
  scale_fill_manual(values = c("gray",
                               "#56B4E9")) +
  scale_y_discrete(limits = rev) +
  theme(legend.position = "right",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_blank()) +
  ylab("")
figure_1c_down

#Save as PDF
if(!dir.exists("Figures/supplement/")) {
  dir.create("Figures/supplement/")
}
if(!dir.exists("Figures/supplement/requests")) {
  dir.create("Figures/supplement/requests")
}

pdf("Figures/supplement/requests/figure_1c_down.pdf",
    width = 8,
    height = 8)
figure_1c_down
dev.off()
```

### From Fig 3D

```{r Fig3_down}
load("Data/OS_list_CCR.RData")

OS_list$OS17_TL <- OS_list$OS17_TL %>%
    FindClusters(resolution = 0.2)
OS_list$`143B_TL` <- OS_list$`143B_TL` %>%
    FindClusters(resolution = 0.15)
OS_list$NCHOS2_TL <- OS_list$NCHOS2_TL %>%
    FindClusters(resolution = 0.15)
OS_list$NCHOS7_TL <- OS_list$NCHOS7_TL %>%
    FindClusters(resolution = 0.2)

# Cluster distribution in each sample
# proportion of cells in lung in each cluster
lung_prop_tbl <- tibble()
for (s_obj in c(OS_list$OS17_TL,
                OS_list$`143B_TL`,
                OS_list$NCHOS2_TL,
                OS_list$NCHOS7_TL)) {
    s_obj$src <- factor(s_obj$src, levels = c("Tibia", "Lung"))
    lung_prop_tbl <-
        s_obj@meta.data %>%
        as_tibble() %>%
        select(model, src, seurat_clusters) %>%
        group_by(model, seurat_clusters) %>%
        summarize(lung_perc = sum(grepl("Lung", src)) / n() * 100,
                  .groups = "drop") %>%
        dplyr::rename(sample = model,
                      cluster = seurat_clusters) %>%
        bind_rows(lung_prop_tbl)
}

write_tsv(lung_prop_tbl, "Data/Supplement/requests/Fig3d_lung_prop_tbl.tsv")
```

```{r Fig3d_down_process, dependson='Fig3_down'}
library(data.table)
library(msigdbr)
library(clusterProfiler)
source("Downstream.v3.R")

#library(pheatmap)
#Create list of objects used in Figure 3 heatmap
B.list <- list(OS17 = OS_list$`OS17_TL`,
               t143b = OS_list$`143B_TL`,
               OS2 = OS_list$`NCHOS2_TL`,
               OS7 = OS_list$`NCHOS7_TL`)
# Use DGEA function to find hallmark pathway expression (downregulated)
em.hm.list <- list()
set.seed(1337)
for (i in seq_len(length(B.list))) {
  em.hm.list[[i]] <- DGEA(B.list[[i]],
                          direction = "down") %>%
    setDT(keep.rownames = TRUE) %>%
    as_tibble() %>%
    rename_with(.cols = where(is.numeric), ~ str_c(names(B.list[i]), "_", .x))
}

cx.pd <- inner_join(inner_join(em.hm.list[[1]],
                               em.hm.list[[2]],
                               by = "rn"),
                    inner_join(em.hm.list[[3]],
                               em.hm.list[[4]],
                               by = "rn"),
                    by = "rn")

head(cx.pd)
##Write table to edit rownames to easily convert to dataframe
write_tsv(cx.pd, file = "Data/Supplement/requests/cxpd.tsv")
```

```{r Fig3d_down}
lung_percent <-
    read_tsv("Data/Supplement/requests/Fig3d_lung_prop_tbl.tsv") %>%
    mutate(sample_order = rank(sample) * 1000 - lung_perc,
           x_label = paste0(sample, " c", cluster) %>%
                reorder(sample_order))
lung_percent$sample <- factor(lung_percent$sample,
                              levels = c("OS-17",
                                         "143B",
                                         "NCH-OS-2",
                                         "NCH-OS-7"))
cx.pd <-
    read_tsv("Data/Supplement/requests/cxpd.tsv",
             show_col_types = FALSE) %>%
    pivot_longer(c(-rn),
                 names_to = "sample",
                 values_to = "pval") %>%
    mutate(cluster = str_remove(sample, ".+_") %>%
                as.numeric(),
           sample = str_replace(sample, "OS2", "NCH-OS-2") %>%
                str_replace("OS7", "NCH-OS-7") %>%
                str_replace("t143b", "143B") %>%
                str_replace("OS17", "OS-17") %>%
                str_remove("_.+"),
            rn = str_remove(rn, "HALLMARK_") %>%
                str_replace("_", " ")) %>%
    left_join(lung_percent) %>%
    mutate(pval = -log10(pval),
           sample_order = rank(sample) * 1000 - lung_perc,
           signif_bool = pval >= (-1 * log10(0.05)),
           signif = if_else(signif_bool == TRUE, "p-value < 0.05", "p-value > 0.05"))

cx.pd$sample <- factor(cx.pd$sample,
                              levels = c("OS-17",
                                         "143B",
                                         "NCH-OS-2",
                                         "NCH-OS-7"))
P1 <-
  cx.pd %>%
  ggplot(., aes(x = x_label, y = rn, fill = signif)) +
    geom_tile() +
    geom_text(aes(label = sprintf("%0.2f", pval)), size = 3) +
    scale_fill_manual(values = c("#56B4E9", "gray"), name = "") +
    scale_y_discrete(limits = rev, position = "right") +
    theme(panel.background = element_blank(),
          strip.background = element_blank(),
          strip.text.x = element_blank(),
          plot.margin = unit(c(0, 0, 0, 0), "cm"),
          panel.border = element_rect(colour = "black", fill = NA),
          axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust = 0.4,
                                     size = 12),
          axis.text.y = element_text(size = 12)) +
    labs(x = "",
         y = "") +
    facet_grid(. ~ sample,
               scales = "free_x",
               space = "free")

P2 <-
  lung_percent %>%
  ggplot(., aes(x = x_label, y = lung_perc, fill = sample)) +
    geom_bar(stat = "identity") +
    ylab("Percent\nin lung") +
    xlab("") +
    labs(fill = "Model",
         y = "Percent in\nmetastases(%)") +
    scale_fill_brewer(palette = "Set2") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          panel.background = element_blank(),
          strip.background = element_blank(),
          strip.text.x = element_blank(),
          plot.margin = unit(c(0, 0, -0.2, 0), "cm"),
          panel.border = element_rect(colour = "black", fill = NA)) +
    facet_grid(. ~ sample,
               scales = "free_x",
               space = "free")

figure_3d_down <- cowplot::plot_grid(P2,
                   P1,
                   ncol = 1,
                   align = "v",
                   rel_heights = c(2, 10))
figure_3d_down

pdf("Figures/supplement/requests/figure_3d_down.pdf",
       width = 14,
       height = 7)
figure_3d_down
dev.off()
```

### From Fig 5F

```{r Fig5f_down}
#source("scSeurat.R")
library(rrrSingleCellUtils)
library(data.table)
library(msigdbr)
library(clusterProfiler)
source("Downstream.v3.R")

load("Data/lung_sub_enr.RData")

B.list <- list(OS17 = lung_sub_enr)

# Perform differential gene expression analysis for hallmark pathways
set.seed(100)
em.hm.list <- list()
for (i in 1:length(B.list)) {
  em.hm.list[[i]] <- DGEA(B.list[[i]],
                          direction = "down")
}

# Achieve -log10 of the p-value
temp1 <- em.hm.list[[1]]
temp1 <- -log10(temp1)

# Remove unnecessary characters
c <- rownames(temp1)
c <- gsub("_", " ", c)
c <- gsub("HALLMARK ", "", c) 
rownames(temp1) <- c

# Remove negative 0 values due to number of decimal points displayed
temp1 <- abs(temp1)
temp1 <- temp1["Enriched clones"]

# Subset to only pathways with significant pvalue (-1 * log10(0.05)) = 1.30103)
temp1["Enriched clones"] >= 1.30103 -> temp1$logical
colnames(temp1) <- c("pvalue",
                     "significant")
logical <- temp1$significant
temp1 <- temp1[logical, ]
temp1 <- rownames_to_column(temp1,
                            var = "Pathway")

# Order by -log10pvalue for plotting
temp1$Pathway <- factor(temp1$Pathway,
                        levels = temp1$Pathway[order(temp1$pvalue)])

# Plot 
figure_5f_down <- ggplot(temp1,
                    aes(x = pvalue,
                        y = Pathway)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  ylab("") +
  xlab("-log10(p-value)") +
  geom_vline(xintercept = 1.5,
               color = "gray",
               linetype = 2,
               size = 1) +
  ggtitle("Downregulated Pathways in Enriched Clonal Populations")

figure_5f_down

pdf("Figures/supplement/requests/figure_5f_down.pdf",
        width = 8,
        height = 3)
figure_5f_down
dev.off()
```

## S33: Test if clones tend to be enriched within a single cluster

```{r bootstrapCloneEnrich, eval=TRUE}
load("Data/os17.RData")
n_clones <- 30
top_clones <-
  tibble(src = os17$src,
         clone = os17$lt,
         cluster = paste0("c_", os17$seurat_clusters)) %>%
  na.omit() %>%
  group_by(src, clone, cluster) %>%
  mutate(n_cells_in_cluster = n()) %>%
  ungroup() %>%
  group_by(src, clone) %>%
  mutate(n_cells_in_clone = n()) %>%
  ungroup() %>%
  group_by(src) %>%
  arrange(desc(n_cells_in_clone), clone, .by_group = TRUE) %>%
  unique() %>%
  pivot_wider(id_cols = c(src, clone, n_cells_in_clone),
              names_from = cluster,
              values_from = n_cells_in_cluster) %>%
  filter(n_cells_in_clone >= 5) %>%
  #slice_head(n = n_clones) %>%
  ungroup() %>%
  pivot_longer(cols = starts_with("c_"),
               names_to = "cluster",
               values_to = "n_cells_in_cluster") %>%
  na.omit() %>%
  mutate(clone = fct_reorder(clone, n_cells_in_clone))
cell_counts_to_test <-
  top_clones$n_cells_in_clone %>%
  unique() %>%
  sort()
names(cell_counts_to_test) <- paste0("c_", cell_counts_to_test)
# Get the number of clusters represented by n randomly selected cells
max_per_cluster_n_cells <- function(sobj, n_cells) {
  selected_cell_indexes <-
    sample(seq_len(ncol(sobj)),
           size = n_cells,
           replace = FALSE)
  max_n_per_cluster <-
    sobj$seurat_clusters[selected_cell_indexes] %>%
    table() %>%
    max()
  return(max_n_per_cluster)
}
# bootstrap the data n_perms times to get a distribution
perm_n_times <- function(sobj, n_cells, n_perms = 1000) {
  message("Testing ", n_cells, " cells")
  n_clusters <- replicate(n_perms, max_per_cluster_n_cells(sobj, n_cells))
  return(n_clusters)
}
# Do this for all observed cells/clone
n_perms <- 100000
perm_table <-
  sapply(cell_counts_to_test,
         function(cell_num) perm_n_times(os17,
                                         cell_num,
                                         n_perms = n_perms),
         USE.NAMES = TRUE) %>%
  as.data.frame()

contin <- top_clones %>%
  rowwise() %>%
  mutate(p_value = (
    sum(perm_table[[paste0("c_", n_cells_in_clone)]] >=
          n_cells_in_cluster) /
      n_perms
  ),
  signif = p_value <= 0.05,
  cluster = str_remove(cluster, "c_")) %>%
  ggplot(aes(x = cluster,
             y = clone,
             fill = n_cells_in_cluster)) +
  geom_tile() +
  geom_text(aes(label = n_cells_in_cluster,
                color = signif)) +
  scale_color_manual(values = c("gray", "red")) +
  facet_grid(src ~ ., scales = "free_y", space = "free") +
  labs(x = "Cluster",
       y = "Clone",
       fill = "Cells in Cluster",
       color = "Enriched Cluster") +
  theme_bw()
contin

pdf("Figures/supplement/requests/contingencytable.pdf",
    height = 10,
    width = 6)
contin
dev.off()
```

## S4-5: Replicates in other cell lines (mouse)

### K7M2

```{r load_k7m2}
# k7m2
#mm10 prefix not present 
## Culture
k7m2.cx.raw <- tenx_load_qc( "/gpfs0/home2/gdrobertslab/lab/Counts/S0201/filtered_feature_bc_matrix/")

k7m2.cx.raw <- subset(k7m2.cx.raw,
                      subset = nFeature_RNA > 2500 &
                        nCount_RNA < 40000 &
                        percent.mt < 10)

k7m2.cx.raw$src <- "Culture"
k7m2.cx.raw$model <- "K7M2"

## Tibia
k7m2.tib.raw <- tenx_load_qc( "/gpfs0/home2/gdrobertslab/lab/Counts/S0105/filtered_feature_bc_matrix/")
k7m2.tib.raw <- subset(k7m2.tib.raw,
                       subset = nFeature_RNA > 800 &
                         nCount_RNA < 50000 &
                         percent.mt < 15)

k7m2.tib.raw$src <- "Tibia"
k7m2.tib.raw$model <- "K7M2"

## Lung
k7m2.lung.raw <- tenx_load_qc( "/gpfs0/home2/gdrobertslab/lab/Counts/S0075-Balb-C-K7M2-Epcam-enriched/outs/filtered_feature_bc_matrix/")

k7m2.lung.raw <- subset(k7m2.lung.raw,
                        subset = nFeature_RNA > 1000 &
                          nCount_RNA < 55000 &
                          percent.mt < 15)

k7m2.lung.raw$src <- "Lung"
k7m2.lung.raw$model <- "K7M2"
```

```{r k7m2_1}
set.seed(100)
# Merge into a single Seurat object
k7m2_all <- merge(k7m2.cx.raw,
              y = c(k7m2.tib.raw,
                    k7m2.lung.raw),
              add.cell.ids = c("Culture",
                               "Tibia",
                               "Lung"),
              project = "LineageTracing")

# Process and cluster
k7m2_all <- NormalizeData(k7m2_all) %>%
  FindVariableFeatures(selection.method = "vst") %>%
  ScaleData() %>%
  RunPCA(pc.genes = k7m2_all@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

if (!dir.exists("Data/Supplement/k7m2_f420")) {
  dir.create("Data/Supplement/k7m2_f420")
}
save(k7m2_all, file = "Data/Supplement/k7m2_f420/k7m2_all.RData")
```

```{r k7m2_2}
load("Data/Supplement/k7m2_f420/k7m2_all.RData")
set.seed(100)

DimPlot(k7m2_all,
        reduction = "umap",
        pt.size = 1,
        label = T) +
  coord_fixed() +
  ggtitle("K7M2 Clusters") +
  scale_color_npg(alpha = 0.7)

# Select clusters showing broad expression of all of the following OS gene markers
FeaturePlot(k7m2_all,
	    features = c("Col1a1","Col1a2","Spp1"),
	    split.by = "src")

set.seed(100)
k7m2 <- subset(k7m2_all, subset = RNA_snn_res.0.3  == c(1, 2, 3))

k7m2 <- RunPCA(k7m2, pc.genes = k7m2@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

# CCR Attempt to regress out the effects of cell cycle on these tumor cells
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

## Convert mouse genes to human genes
library(nichenetr)
s.genes <- convert_human_to_mouse_symbols(s.genes)
g2m.genes <- convert_human_to_mouse_symbols(g2m.genes)

k7m2 <- CellCycleScoring(object = k7m2 ,
                         s.features = s.genes,
                         g2m.features = g2m.genes,
                         set.ident = TRUE)

k7m2 <- ScaleData(object = k7m2,
                  vars.to.regress = c("S.Score",
                                      "G2M.Score"),
                  features = rownames(x = k7m2))

k7m2 <- RunPCA(k7m2, pc.genes = k7m2@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

## Find the optimum resolution for clustering
k7m2.nC <- nRes(k7m2,
                res = seq(from = 0.1,
                          to = 0.3,
                          by = 0.01))

k7m2 <- RunPCA(k7m2, pc.genes = k7m2@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.1)

save(k7m2, file = "Data/Supplement/k7m2_f420/k7m2_tumor.RData")
```

```{r k7m2_3}
load("Data/Supplement/k7m2_f420/k7m2_tumor.RData")

DimPlot(k7m2,
        reduction = "umap",
        pt.size = 1,
        label = T) +
  coord_fixed() +
  ggtitle("K7M2 Clusters") +
  scale_color_npg(alpha = 1)

DimPlot(k7m2,
        reduction = "umap",
        pt.size = 1,
	group.by = "Phase") +
  coord_fixed() +
  ggtitle("K7M2 by Phase") +
  scale_color_npg(alpha = 1)


# Plot the data
set.seed(100)
k7m2$src <- factor(k7m2$src,
                      levels = c("Culture",
				 "Tibia",
                                 "Lung"))
k7m2_umap <- DimPlot(k7m2,
        reduction = "umap",
        group.by = "src",
        pt.size = 1,
        label = F) +
  coord_fixed() +
  ggtitle("K7M2 by Source") +
  scale_color_npg()
k7m2_umap

if (!dir.exists("Figures/supplement/requests/k7m2_f420")) {
  dir.create("Figures/supplement/requests/k7m2_f420")
}
pdf("Figures/supplement/requests/k7m2_f420/k7m2_umap.pdf",
        width = 5,
        height = 5)
k7m2_umap
dev.off()
```

### F420

```{r load_F420}
# F420
#mm10 already taken out
## Culture
F420.cx.raw <- tenx_load_qc( "/gpfs0/home2/gdrobertslab/lab/Counts/S0200/filtered_feature_bc_matrix/")

F420.cx.raw <- subset(F420.cx.raw,
                      subset = nFeature_RNA > 2500 &
                        nCount_RNA < 45000 &
                        percent.mt < 10)

F420.cx.raw$src <- "Culture"
F420.cx.raw$model <- "F420"

## Tibia
F420.tib.raw <- tenx_load_qc( "/gpfs0/home2/gdrobertslab/lab/Counts/S0104/filtered_feature_bc_matrix/")
F420.tib.raw <- subset(F420.tib.raw,
                       subset = nFeature_RNA > 800 &
                         nCount_RNA < 50000 &
                         percent.mt < 15)

F420.tib.raw$src <- "Tibia"
F420.tib.raw$model <- "F420"

## Lung
## double check usage throughout rest of code - Emily (fixed issue where coming from wrong file)
## (Epcam postive:negative 1:1)
F420.lung.raw <- tenx_load_qc("/gpfs0/home2/gdrobertslab/lab/Counts/S0067/filtered_feature_bc_matrix/")

F420.lung.raw <- subset(F420.lung.raw,
                        subset = nFeature_RNA > 1500 &
                          nCount_RNA < 55000 &
                          percent.mt < 20)

F420.lung.raw$src <- "Lung"
F420.lung.raw$model <- "F420"

# Second sample from same batch (Epcam postive:negative 1:1)
F420.lung.raw2 <- tenx_load_qc("/gpfs0/home2/gdrobertslab/lab/Counts/S0068/filtered_feature_bc_matrix/")

F420.lung.raw2 <- subset(F420.lung.raw2,
                        subset = nFeature_RNA > 1500 &
                          nCount_RNA < 45000 &
                          percent.mt < 20)

F420.lung.raw2$src <- "Lung"
F420.lung.raw2$model <- "F420"
```

```{r F420_2}
# Merge into a single Seurat object
F420_all <- merge(F420.cx.raw,
              y = c(F420.tib.raw,
                    F420.lung.raw,
                    F420.lung.raw2),
              add.cell.ids = c("Culture",
                               "Tibia",
                               "Lung",
                               "Lung"),
              project = "LineageTracing")

# Process and cluster
F420_all <- NormalizeData(F420_all) %>%
  FindVariableFeatures(selection.method = "vst") %>%
  ScaleData() %>%
  RunPCA(pc.genes = F420_all@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

DimPlot(F420_all,
        reduction = "umap",
        pt.size = 1,
        label = T) +
  coord_fixed() +
  ggtitle("F420 Clusters") +
  scale_color_npg(alpha = 0.7)

# Plot the data
set.seed(100)
cell.ids <- sample(colnames(F420_all))

DimPlot(F420_all,
        reduction = "umap",
        group.by = "src",
        pt.size = 1,
        label = F,
        order = cell.ids) +
  coord_fixed() +
  ggtitle("F420 by Source") +
  scale_color_npg()

save(F420_all, file = "Data/Supplement/k7m2_f420/F420_all.RData")
```

```{r f420_3}
load("Data/Supplement/k7m2_f420/F420_all.RData")

DimPlot(F420_all,
        reduction = "umap",
        pt.size = 1,
        label = T) +
  coord_fixed() +
  ggtitle("F420 Clusters") 

# Select clusters showing broad expression of all of the following OS gene markers
FeaturePlot(F420_all,
	    features = c("Col1a1","Col1a2","Spp1"),
	    split.by = "src")

set.seed(100)
F420 <- subset(F420_all, subset = RNA_snn_res.0.3  == c(0, 3))

F420 <- RunPCA(F420 , pc.genes = F420 @var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

# CCR Attempt to regress out the effects of cell cycle on these tumor cells
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

# Convert mouse genes to human genes
library(nichenetr)
s.genes <- convert_human_to_mouse_symbols(s.genes)
g2m.genes <- convert_human_to_mouse_symbols(g2m.genes)

F420 <- CellCycleScoring(object = F420,
                         s.features = s.genes,
                         g2m.features = g2m.genes,
                         set.ident = TRUE)

F420 <- ScaleData(object = F420,
                  vars.to.regress = c("S.Score",
                                      "G2M.Score"),
                  features = rownames(x = F420))

F420 <- RunPCA(F420, pc.genes = F420@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.3)

## Find the optimum resolution for clustering
F420.nC <- nRes(F420,
                res = seq(from = 0.1,
                          to = 0.3,
                          by = 0.01))

F420 <- RunPCA(F420, pc.genes = F420@var.genes, npcs = 20) %>%
  RunUMAP(reduction = "pca", dims = 1:20) %>%
  FindNeighbors(reduction = "pca", dims = 1:20) %>%
  FindClusters(resolution = 0.1)
save(F420, file = "Data/Supplement/k7m2_f420/F420_tumor.RData")
```

```{r f420_4}
DimPlot(F420,
        reduction = "umap",
        pt.size = 1,
        label = T) +
  coord_fixed() +
  ggtitle("F420 Clusters") +
  scale_color_npg(alpha = 1)

DimPlot(F420,
        reduction = "umap",
        pt.size = 1,
	group.by = "Phase") +
  coord_fixed() +
  ggtitle("F420 by Phase") +
  scale_color_npg(alpha = 1)


# Plot the data
set.seed(100)
F420$src <- factor(F420$src,
                   levels = c("Culture",
                              "Tibia",
                              "Lung"))

F420_umap <- DimPlot(F420,
        reduction = "umap",
        group.by = "src",
        pt.size = 1,
        label = F) +
  coord_equal() +
  ggtitle("F420 by Source") +
  scale_color_npg()
F420_umap

pdf("Figures/supplement/requests/k7m2_f420/F420_umap.pdf",
        width = 5,
        height = 5)
F420_umap
dev.off()
```

### S4

#### Repeat Figure 2C

C) Pathway enrichment analysis with adjusted p values for hallmark gene sets associated with these shared genes

```{r Fig2c_repeat, fig.height=6, fig.width=10}
### Extract genes that are shared in these datasets
load("Data/OS_merged_postCCR.RData")
load("Data/Supplement/k7m2_f420/k7m2_tumor.RData")
load("Data/Supplement/k7m2_f420/F420_tumor.RData")
# Set seed before subsetting and DE analysis to ensure consistent results
set.seed(100)

k7m2$tissue <- k7m2$src 
F420$tissue <- F420$src 
data <- list(
  os17 = subset(OS, cells = WhichCells(OS, expression = model == "OS-17")),
  t143B = subset(OS, cells = WhichCells(OS, expression = model == "143B")),
  NCHOS2 = subset(OS, cells = WhichCells(OS, expression = model == "NCH-OS-2")),
  NCHOS7 = subset(OS, cells = WhichCells(OS, expression = model == "NCH-OS-7")),
  K7M2 = k7m2,
  F420 = F420)

# Rename "Flank" to "Culture" for NCHOS samples
data$NCHOS2$src[grep("Flank", data$NCHOS2$src)] <- "Culture"
data$NCHOS7$src[grep("Flank", data$NCHOS7$src)] <- "Culture"

repfig_2c_data <- tibble()

for (tissue_name in c("Tibia", "Lung")) {
  for (direction in c("up", "down")) {
    de_genes <- list()
    for (i in seq_len(length(data))) {
      Idents(data[[i]]) <- data[[i]]$src
      de_genes[[i]] <-
        FindMarkers(data[[i]],
                    ident.1 = tissue_name,
                    ident.2 = "Culture",
                    only.pos = FALSE,
                    min.pct = 0.1)
      
      # Select genes that change in the direction of interest
      if (direction == "up") {
        de_genes[[i]] <-
          de_genes[[i]][de_genes[[i]]$p_val_adj < 0.05 &
                          de_genes[[i]]$avg_log2FC > 0, ]
      } else {
        de_genes[[i]] <-
          de_genes[[i]][de_genes[[i]]$p_val_adj < 0.05 &
                          de_genes[[i]]$avg_log2FC < 0, ]
      }
    }
    
    # Convert mouse genes to human genes
    library(nichenetr)
    rowname_df <- tibble(human = convert_mouse_to_human_symbols(rownames(de_genes[[5]])) %>%
                           as.vector(), mouse = rownames(de_genes[[5]]))
    #table(rowname_df$human) %>% as.data.frame() %>% arrange(Freq * -1) %>% head()
    #summary(as.factor(rowname_df$human))
    
    rowname_df2 <- tibble(human = convert_mouse_to_human_symbols(rownames(de_genes[[6]])) %>%
                            as.vector(), mouse = rownames(de_genes[[6]]))
    #table(rowname_df2$human) %>% as.data.frame() %>% arrange(Freq * -1) %>% head()
    #summary(as.factor(rowname_df2$human))
    
    if (tissue_name == "Tibia") {
      if (direction == "up") {
        de_genes[[5]] <- de_genes[[5]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) & human != "HLA-A") %>%
          column_to_rownames("human")
        message(paste0(tissue_name, " ",  direction, " for de_genes_5 DONE"))
        
        de_genes[[6]] <- de_genes[[6]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df2) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) &
                   human != "BCL2A1" & 
                   human != "CSF2RB" & 
                   human != "FTL" &
                   human != "IFI16" & 
                   human != "LILRB3" &
                   human != "LILRB4" &
                   human != "HLA-A" & 
                   human != "HBA2" &
                   human != "H3F3A") %>% 
          column_to_rownames("human")
        message(paste0(tissue_name, " ",  direction, " for de_genes_6 DONE"))

      } else {
        de_genes[[5]] <- de_genes[[5]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) & human != "EIF3J") %>%
          column_to_rownames("human")
	message(paste0(tissue_name, " ",  direction, " for de_genes_5 DONE"))
        
        de_genes[[6]] <- de_genes[[6]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df2) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human)) %>%
          column_to_rownames("human")
	message(paste0(tissue_name, " ",  direction, " for de_genes_6 DONE"))

      }
    }
    
    if (tissue_name == "Lung") {
      if (direction == "up") {
        de_genes[[5]] <- de_genes[[5]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) & human != "H3F3A" & human != "HLA-A") %>%
          column_to_rownames("human")
	message(paste0(tissue_name, " ", direction, " for de_genes_5 DONE"))
        
        de_genes[[6]] <- de_genes[[6]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df2) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) &
                   human != "CSF2RB" & 
                   human != "FTL" & 
                   human != "IFI16" &
                   human != "LILRB3" & 
                   human != "LILRB4" & 
                   human != "LYZ" &
                   human != "MTHFS" &
                   human != "BCL2A1" & 
                   human != "BCL2A1" & 
                   human != "HLA-A" & 
                   human != "HBA2") %>% 
          column_to_rownames("human")
        message(paste0(tissue_name, " ", direction, " for de_genes_6 DONE"))

      } else {
        de_genes[[5]] <- de_genes[[5]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) & human != "EIF3J") %>%
          column_to_rownames("human")
	message(paste0(tissue_name, " ", direction, " for de_genes_5 DONE"))
        
        de_genes[[6]] <- de_genes[[6]] %>%
          rownames_to_column("mouse") %>%
          full_join(rowname_df2) %>%
          select(-mouse) %>%
          # Remove NAs and repeated genes (from conversion)
          filter(!is.na(human) & human != "AMD1") %>%
          column_to_rownames("human")
	message(paste0(tissue_name, " ", direction, " for de_genes_6 DONE"))
      }
    }
    
    A <- rownames(de_genes[[1]])
    B <- rownames(de_genes[[2]])
    C <- rownames(de_genes[[3]])
    D <- rownames(de_genes[[4]])
    E <- rownames(de_genes[[5]]) #interprets F as FALSE
    G <- rownames(de_genes[[6]])
    
    # Preparing clusterProfiler to perform hypergeometric test on msigdb signatures 
    msig.gene.set <-
      msigdbr(species = "Homo sapiens", category = "H") %>%
      dplyr::select(gs_name, human_gene_symbol)
    msig.name <-
      msigdbr(species = "Homo sapiens", category = "H") %>%
      dplyr::select(gs_id, gs_name)
    
    # Getting the middle "flower" of the Venn diagram
    intersect <- c((intersect(intersect(intersect(intersect(intersect(A, B), C), D), E), G)),
                   (intersect(intersect(intersect(A, B), C), E)),
                   (intersect(intersect(intersect(A, B), D), E)),
                   (intersect(intersect(intersect(A, D), C), E)),
                   (intersect(intersect(intersect(B, D), C), E)),
                   (intersect(intersect(intersect(A, B), C), G)),
                   (intersect(intersect(intersect(A, B), D), G)),
                   (intersect(intersect(intersect(A, D), C), G)),
                   (intersect(intersect(intersect(B, D), C), G)),
                   (intersect(intersect(intersect(A, B), E), G)),
                   (intersect(intersect(intersect(A, C), E), G)),
                   (intersect(intersect(intersect(A, D), E), G))) %>%
      unique()
    temp <-
      enricher(intersect,
               TERM2GENE = msig.gene.set,
               TERM2NAME = msig.name)@result %>%
      as_tibble() %>%
      select(ID, p.adjust) %>%
      dplyr::rename(pathway = ID) %>%
      arrange(p.adjust) %>%
      slice_head(n = 5) %>%
      mutate(tissue = tissue_name,
             label_y = if_else(direction == "up", 1, -1),
             p.adjust = -log10(p.adjust),
             pathway = str_remove(pathway, "HALLMARK_") %>%
               str_replace_all("_", " "),
             order = seq_len(n()))
    
    if (direction == "down") {
      temp$p.adjust <- temp$p.adjust * -1
    }
    
    repfig_2c_data <- bind_rows(repfig_2c_data, temp)
  }
}

################ Two sided Barplot
# Properly name and order tissue types
repfig_2c_data <- repfig_2c_data %>%
  mutate(tissue = stringr::str_replace(tissue,
                                       "Tibia",
                                       "Tibia Gene Sets") %>%
           stringr::str_replace("Lung",
                                "Lung Gene Sets"),
         pathway = stringr::str_wrap(pathway, width = 25))

repfig_2c_data$tissue <- factor(as.factor(repfig_2c_data$tissue),
                                levels = c("Tibia Gene Sets",
                                           "Lung Gene Sets"))

# Create tibble that is properly ordered to use in proper plotting of tibia and lung data
lab4plot <- tibble(y = c(-4, 4, -4, 4),
                   x = c(0.5, 0.5, 0.5, 0.5),
                   label = as.factor(c("Downregulated\nduring tibia colonization",
                                       "Upregulated\nduring tibia colonization",
                                       "Downregulated\nduring lung colonization",
                                       "Upregulated\nduring lung colonization")),
                   tissue = factor(as.factor(c("Tibia Gene Sets",
                                               "Tibia Gene Sets",
                                               "Lung Gene Sets",
                                               "Lung Gene Sets")),
                                   levels = c("Tibia Gene Sets",
                                              "Lung Gene Sets")))

repfigure_2c <-
  ggplot() +
  geom_bar(data = repfig_2c_data,
           aes(x = -1 * order,
               y = p.adjust,
               fill = p.adjust > 0),
           stat = "identity",
           alpha = 0.8) +
  coord_flip() +
  facet_wrap(~tissue,
             ncol = 1,
             scales = "free") +
  geom_hline(yintercept = c(-1.5, 1.5),
             color = "gray",
             linetype = 2,
             size = 1) +
  geom_hline(yintercept = 0,
             color = "black",
             size = 1) +
  geom_text(data = repfig_2c_data,
            aes(x = -1 * order,
                y = label_y * 4,
                label = pathway)) +
  geom_text(data = lab4plot,
            aes(x = x,
                y = y,
                label = label),
            fontface = "bold") +
  scale_fill_manual(values = c("#377eb8", "#e41a1c")) +
  theme_bw() +
  theme(strip.background = element_rect(color = "white", fill = "white"),
        strip.text.x = element_text(size = 20, face = "bold"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  labs(title = "",
       y = "",
       x = "") +
  ylim(-12, 12) +
  xlim(-6, 1)

repfigure_2c

# Save as PDF
pdf("Figures/supplement/requests/k7m2_f420/repeat_figure_2c_k7m2-f420.pdf",
    width = 8,
    height = 9)
repfigure_2c
dev.off()

save(repfigure_2c, file = "Figures/supplement/requests/k7m2_f420/repfigure_2c.RData")
save(repfigure_2c, file = "Figures/supplement/requests/k7m2_f420/repfig_2c_data.RData")
```

### S6: Repeat Figure 3

Figure 3. Osteosarcoma cells retain phenotypic heterogeneity despite adaptive changes in response to changing microenvironments. A) Schematic outlining scRNA-seq bioinformatics workflow. In B-D, samples subset to equal number of cells to allow inter-model comparison (n=1500 per condition).

##### S6A (Repeat Figure 3B)

*statistical method subject to change*

B) Osteosarcoma cells exhibited higher ITH scores than that of osteoblasts grown in cell culture (####p <0.001). Within model comparison, identified increase in ITH scores as cells colonize tibia and lung microenvironments, with the exception of NCH-OS-7 (****p <0.001). P values were adjusted using Šidák multiple comparisons test.

```{r repfig3b, eval=T}
# Load OS data that was subsetted to be equal cell numbers across 
load("Data/Supplement/k7m2_f420/F420_tumor.RData")
load("Data/Supplement/k7m2_f420/k7m2_tumor.RData")
load("Data/OSlist_1500_CCR.RData")

if (!dir.exists("Figures/supplement/requests/k7m2_f420")) {
  dir.create("Figures/supplement/requests/k7m2_f420")
}

# Subset down to 1500 where we can, but some of the groups have much less.
# Since we're calculating the average distance, the results will still be valid,
# just less accurate than if we had greater numbers.
OS.list_1500 <- c(OS.list_1500,
                  k7m2.Cx = subset(k7m2, subset = src == "Culture") %>%
                      subset(cells = sample(Cells(.), 1500)),
                  k7m2.Tibia = subset(k7m2, subset = src == "Tibia"),
                  k7m2.Lung = subset(k7m2, subset = src == "Lung"),
                  F420.Cx = subset(F420, subset = src == "Culture") %>%
                      subset(cells = sample(Cells(.), 1500)),
                  F420.Tibia = subset(F420, subset = src == "Tibia") %>%
                      subset(cells = sample(Cells(.), 1500)),
                  F420.Lung = subset(F420, subset = src == "Lung"))

# k7m2d[, .N, by = "src"]
#       src    N
#1: Culture 2271
#2:   Tibia  605
#3:    Lung  393
# F420.d[, .N, by = "src"]
#       src    N
#1: Culture 4067
#2:   Tibia 1600
#3:    Lung   15

merged_sobjs <-
    SeuratObject:::merge.Seurat(x = OS.list_1500[[1]],
                                y = OS.list_1500[2:length(OS.list_1500)],
                                add.cell.ids = names(OS.list_1500))

merged_sobjs$group <-
    colnames(merged_sobjs) %>%
    str_remove("_.+")

table(merged_sobjs$group)

group_names <-
    merged_sobjs$group %>%
    unique()

better_names <-
    group_names %>%
    str_replace("\\.", " ") %>%
    str_replace("cx", "Culture") %>%
    str_replace("Cx", "Culture") %>%
    str_replace("OB", "Osteoblast") %>%
    str_replace("t143", "143")

# Find the genes that are variable within each sample, and combine into a vector
var_genes <- c()
for (i in seq_along(OS.list_1500)) {
  OS.list_1500[[i]] <- OS.list_1500[[i]] %>%
    NormalizeData(verbose = FALSE) %>%
    FindVariableFeatures(verbose = FALSE)
  var_genes <- c(var_genes, VariableFeatures(OS.list_1500[[i]]))
}
var_genes <- unique(var_genes)

raw_ith <-
    tibble(Sample = character(),
           ITHScore = numeric())

for (i in seq_len(length(group_names))) {
  better_name <-
    better_names[i]

  sub_obj <- subset(merged_sobjs, group == group_names[i]) %>%
    NormalizeData(verbose = FALSE) %>%
    FindVariableFeatures(verbose = FALSE) %>%
    ScaleData(
      features = var_genes,
      vars.to.regress = c("S.Score", "G2M.Score"),
      block.size = 100000)

  z <-
    GetAssayData(sub_obj, slot = "scale.data") %>%
    as.data.frame() %>%
    rownames_to_column("gene") %>%
    as_tibble() %>%
    filter(gene %in% VariableFeatures(sub_obj)) %>%
    select(-gene) %>%
    t()

  mat <- dist(z, diag = TRUE, upper = FALSE)
  mat2 <- as.matrix(mat)
  mat2[upper.tri(mat2, diag = TRUE)] <- NA
  raw_ith <- bind_rows(raw_ith,
                       tibble(ITHScore = colMeans(mat2, na.rm = TRUE),
                              Sample = better_name))
}

# Some rows will be NAs because we removed the diagonal from the matrix
raw_ith <- na.omit(raw_ith)

ith <-
    raw_ith %>%
    mutate(Sample = as.factor(Sample) %>%
                fct_relevel(better_names),
           logITH = log10(ITHScore))

ridge_cols <- c("#B8B8B8FF",
                "#D43F3AFF", "#D43F3AFF", "#D43F3AFF",
                "#EEA236FF", "#EEA236FF", "#EEA236FF",
                "#357EBDFF", "#357EBDFF", "#357EBDFF",
                "#5CB85CFF", "#5CB85CFF", "#5CB85CFF",
                "#E495A5", "#E495A5", "#E495A5",
                "#ABB065", "#ABB065", "#ABB065")
rplot <-
  ggplot(ith, aes(x = logITH, y = Sample)) +
  geom_density_ridges(scale = 2, aes(fill = Sample)) +
  scale_y_discrete(expand = c(0, 0), limits = rev) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_fill_manual(values = ridge_cols) +
  theme_ridges() +
  xlab("log[ITH Score]") +
  theme(legend.position = "none",
        axis.text.y = element_text(size = 6),
        axis.text.x = element_text(size = 6),
        axis.title.x = element_text(hjust = 0.5, size = 6),
        axis.title.y = element_blank())

rplot

# Perform a stat test for each group
comps <- tribble(
  ~ctl, ~comp,
  1, 2,
  1, 5,
  1, 8,
  1, 11,
  1, 14,
  1, 17,
  2, 3,
  2, 4,
  5, 6,
  5, 7,
  8, 9,
  8, 10,
  11, 12,
  11, 13,
  14, 15,
  14, 16,
  17, 18,
  17, 19
)

comps$ctl_name <- unique(ith$Sample)[comps$ctl]
comps$comp_name <- unique(ith$Sample)[comps$comp]
comps$mwu_p <- NA
comps$overlap <- NA
comps$eff_size <- NA

# Calculate Mann-Whitney test (in figure, p-value)
for (i in seq_along(comps$ctl)) {
  message("compare ",
          group_names[comps$ctl[i]],
          " to ",
          group_names[comps$comp[i]])

  comps$mwu_p[i] <-
    wilcox.test(filter(ith, Sample == better_names[comps$ctl[i]]) %>%
                    pull(ITHScore),
                filter(ith, Sample == better_names[comps$comp[i]]) %>%
                    pull(ITHScore))$p.value

  comps$overlap[i] <-
    overlapping::overlap(list(filter(ith, Sample == better_names[comps$ctl[i]]) %>%
                    pull(ITHScore),
                              filter(ith, Sample == better_names[comps$comp[i]]) %>%
                    pull(ITHScore)))$OV

  comps$eff_size[i] <-
    effsize::cohen.d(filter(ith, Sample == better_names[comps$ctl[i]]) %>%
                        pull(ITHScore),
                     filter(ith, Sample == better_names[comps$comp[i]]) %>%
                        pull(ITHScore))$estimate
}

# Convert overlap to percentage
comps$overlap <- paste0(format(comps$overlap * 100, digits = 2, trim = T), "%")

# Plot with p-values superimposed for each sample
repfigure_3b_15 <-
    rplot +
    geom_text(data = comps,
              aes(x = 1.95, y = comp_name),
              label = paste0("overlap ",
                             comps$overlap,
                             ", eff. size ",
                             format(comps$eff_size, digits = 1, trim = T),
                             "\n"),
              size = 1.8,
              nudge_y = 0.2,
              color = "#2b2a2a")
repfigure_3b_15

save(repfigure_3b_15,
     file = "Figures/supplement/requests/k7m2_f420/repfigure_3b_15.RData")

# Save as PDF
pdf("Figures/supplement/requests/k7m2_f420/repfigure_3b_15.pdf",
       width = 8,
       height = 8)
repfigure_3b_15
dev.off()
```

### S6B (Repeat Figure 3C)

C) UMAP analysis for merged primary and metastatic samples in each of the four models. Cells in grey represent remaining cells in merged sample. Cluster enrichment analysis shows distribution of cells in each cluster in the two microenvironment conditions (tibia, lung). While some cells in the primary and metastatic lesions adopted shared phenotypes, others adopted distinct phenotypes.

```{r repfig3c, eval=T}
load("Data/OS_list_CCR.RData")

OS_list$OS17_TL <- OS_list$OS17_TL %>%
    FindClusters(resolution = 0.2)
OS_list$`143B_TL` <- OS_list$`143B_TL` %>%
    FindClusters(resolution = 0.15)
OS_list$NCHOS2_TL <- OS_list$NCHOS2_TL %>%
    FindClusters(resolution = 0.15)
OS_list$NCHOS7_TL <- OS_list$NCHOS7_TL %>%
    FindClusters(resolution = 0.2)
OS_list$K7M2 <-
    k7m2 %>%
    subset(src == "Lung" | src == "Tibia") %>%
    FindVariableFeatures(selection.method = "vst",
                         verbose = FALSE) %>%
    RunPCA(pc.genes = k7m2_all@var.genes,
           npcs = 20,
           verbose = FALSE) %>%
    RunUMAP(reduction = "pca",
            dims = 1:20,
            verbose = FALSE) %>%
    FindNeighbors(reduction = "pca",
                  dims = 1:20,
                  verbose = FALSE) %>%
    FindClusters(resolution = 0.15,
                 verbose = FALSE)

OS_list$F420 <-
    F420 %>%
    subset(src == "Lung" | src == "Tibia") %>%
    FindVariableFeatures(selection.method = "vst",
                         verbose = FALSE) %>%
    RunPCA(pc.genes = k7m2_all@var.genes,
           npcs = 20,
           verbose = FALSE) %>%
    RunUMAP(reduction = "pca",
            dims = 1:20,
            verbose = FALSE) %>%
    FindNeighbors(reduction = "pca",
                  dims = 1:20,
                  verbose = FALSE) %>%
    FindClusters(resolution = 0.15,
                 verbose = FALSE)

# Cluster distribution in each sample
# proportion of cells in lung in each cluster
lung_prop_tbl <- tibble()
for (s_obj in c(OS_list$OS17_TL,
                OS_list$`143B_TL`,
                OS_list$NCHOS2_TL,
                OS_list$NCHOS7_TL,
                OS_list$K7M2,
                OS_list$F420)) {
    s_obj$src <- factor(s_obj$src, levels = c("Tibia", "Lung"))
    lung_prop_tbl <-
        s_obj@meta.data %>%
        as_tibble() %>%
        select(model, src, seurat_clusters) %>%
        group_by(model, seurat_clusters) %>%
        summarize(lung_perc = sum(grepl("Lung", src)) / n() * 100,
                  .groups = "drop") %>%
        dplyr::rename(sample = model,
                      cluster = seurat_clusters) %>%
        bind_rows(lung_prop_tbl)
}

write_tsv(lung_prop_tbl, "Data/Supplement/k7m2_f420/repfig3d_lung_prop_tbl.tsv")

# UMAP Plots
sample_names <- c("OS17_TL",
                  "143B_TL",
                  "NCHOS2_TL",
                  "NCHOS7_TL",
                  "K7M2",
                  "F420")

suppfigure_3c_umaps <- list()
for (s_obj in sample_names) {
  # Tibia plot with lung-derived cells greyed out
  temp_obj <- OS_list[[s_obj]]
  lung <- vector(mode = "character")
  temp <- WhichCells(temp_obj,
                     expression = src == "Lung")
  lung <- append(lung, temp)
  Idents(temp_obj, cells = lung) <- " "

  cust_dim_plot <- function(x_obj, title, subtitle) {
      DimPlot(x_obj,
              reduction = "umap",
              pt.size = 0.1,
              label = T,
              label.size = 1,
              cols = c("light grey",
                       "#E64B35FF",
                       "#4DBBD5FF",
                       "#00A087FF",
                       "#3C5488FF"),
              order = T) +
    coord_fixed() +
    NoLegend() +
    NoAxes() +
    ggtitle(title, subtitle = subtitle) +
    theme(plot.title = element_text(hjust = 0.5, size = 7),
          plot.subtitle = element_text(hjust = 0.5, size = 6),
          plot.margin = unit(c(0, 0, 0, 0), "pt"))
  }

  tibia_plot <-
    cust_dim_plot(temp_obj,
                  title = s_obj,
                  subtitle = "Tibia")

  # Lung plot with lung-derived cells greyed out
  temp_obj2 <- OS_list[[s_obj]]
  tib <- vector(mode = "character")
  temp <- WhichCells(temp_obj2,
                     expression = src == "Tibia")
  tib <- append(tib, temp)
  Idents(temp_obj2, cells = tib) <- " "

  lung_plot <-
    cust_dim_plot(temp_obj2,
                  title = "",
                  subtitle = "Lung")

  suppfigure_3c_umaps[[s_obj]] <- tibia_plot + lung_plot

  pdf(paste0("Figures/supplement/requests/k7m2_f420/umap_",
             s_obj,
             ".pdf"),
      width = 5,
      height = 5)
  print(suppfigure_3c_umaps[[s_obj]])
  dev.off()
}

save(suppfigure_3c_umaps,
     file = "Figures/supplement/requests/k7m2_f420/suppFigure_3c_umaps.RData")

# Bar plots

suppFigure_3c_barplots <- list()

for (sobj in c(OS_list$OS17_TL,
               OS_list$`143B_TL`,
               OS_list$NCHOS2_TL,
               OS_list$NCHOS7_TL,
               OS_list$K7M2,
               OS_list$F420)) {
  tident <- table(Idents(sobj),
                  sobj$src,
                  dnn = c("cluster",
                          "tissue")) %>%
    as.data.frame() %>%
    group_by(tissue) %>%
    mutate(perc = Freq / sum(Freq) * 100)

    tident$cluster <- as.character(tident$cluster)
    tident$tissue <- factor(tident$tissue, levels = c("Tibia", "Lung"))
    suppFigure_3c_barplots[[sobj$model[1]]] <-
        ggplot(tident,
               aes(x = tissue,
                   y = perc,
                   fill = cluster)) +
        theme_bw(base_size = 5) +
        geom_col(position = "fill", width = 0.5) +
        xlab("Tissue") +
        ylab("Proportion") +
        scale_fill_npg(alpha = 1) +
        theme_minimal() +
        NoLegend() +
        theme(panel.grid = element_blank(),
              axis.text.y = element_text(size = 5),
              axis.text.x = element_text(size = 5,
                                         angle = 45,
                                         hjust = 1),
              axis.title.x = element_text(size = 6),
              axis.title.y = element_text(size = 6))

  # Save plot as PDF
  pdf(paste0("Figures/supplement/requests/k7m2_f420/suppFigure_3c2_",
             head(sobj$model, 1),
             ".pdf"),
      width = 1,
      height = 2)
  print(suppFigure_3c_barplots[[sobj$model[1]]])
  dev.off()
}

save(suppFigure_3c_barplots,
     file = "Figures/supplement/requests/k7m2_f420/suppFigure_3c_barplots.RData")
```

### S6C (Repeat Figure 3D)

D) Metabolic heterogeneity in glycolysis activation. We used a pathway enrichment
analysis for hallmark gene sets using genes differentially upregulated in each
cluster relative to every other cluster within the same model. P values were
adjusted for multiple comparisons. Boxes in grey identify non-significant pathway
enrichments, whereas boxes in red identify statistically significant enrichments.
Bar plot shows percentage of cells identified per cluster in lung metastases.
In B-D, samples subset to equal number of cells to allow inter-model comparison
(n=1500 per condition).

```{r repfig3d_process, dependson='repfig3c', eval=T}
source("Downstream.v3.R")
library(data.table)
#library(pheatmap)
#Create list of objects used in repfigure 3 heatmap
B.list <- list(OS17 = OS_list$`OS17_TL`,
               t143b = OS_list$`143B_TL`,
               OS2 = OS_list$`NCHOS2_TL`,
               OS7 = OS_list$`NCHOS7_TL`,
               K7M2 = OS_list$K7M2,
               F420 = OS_list$F420)

species_vect <- c("human", "human", "human", "human", "mouse", "mouse")

# Use DGEA function to find hallmark pathway expression
em.hm.list <- list()
set.seed(1337)
for (i in seq_len(length(B.list))) {
  em.hm.list[[i]] <- DGEA(B.list[[i]], spec = species_vect[i]) %>%
    setDT(keep.rownames = TRUE) %>%
    as_tibble() %>%
    rename_with(.cols = where(is.numeric), ~ str_c(names(B.list[i]), "_", .x))
}

cx.pd <-
    purrr::reduce(em.hm.list, inner_join, by = "rn")

# Make the "rn" column the rownames -Emily
# cx.pd <- cx.pd %>%
#     column_to_rownames("rn")

head(cx.pd)
##Write table to edit rownames to easily convert to dataframe
write_tsv(cx.pd, file = "Data/Supplement/k7m2_f420/cxpd.tsv")
rm(OS_list)
```

```{r repfig3d, eval=T}
lung_percent <-
    read_tsv("Data/Supplement/k7m2_f420/repfig3d_lung_prop_tbl.tsv") %>%
    mutate(sample_order = rank(sample) * 1000 - lung_perc,
           x_label = paste0(sample, " c", cluster) %>%
                reorder(sample_order))

cx.pd <-
    read_tsv("Data/Supplement/k7m2_f420/cxpd.tsv",
             show_col_types = FALSE) %>%
    pivot_longer(c(-rn),
                 names_to = "sample",
                 values_to = "pval") %>%
    mutate(cluster = str_remove(sample, ".+_") %>%
                as.numeric(),
           sample = str_replace(sample, "OS2", "NCH-OS-2") %>%
                str_replace("OS7", "NCH-OS-7") %>%
                str_replace("t143b", "143B") %>%
                str_replace("OS17", "OS-17") %>%
                str_remove("_.+"),
            rn = str_remove(rn, "HALLMARK_") %>%
                str_replace("_", " ")) %>%
    left_join(lung_percent) %>%
    mutate(pval = -log10(pval),
           sample_order = rank(sample) * 1000 - lung_perc,
           signif_bool = pval >= (-1 * log10(0.05)),
           signif = if_else(signif_bool == TRUE, "p-value < 0.05", "p-value > 0.05"))

P1 <-
  cx.pd %>%
  ggplot(., aes(x = x_label, y = rn, fill = signif)) +
    geom_tile() +
    geom_text(aes(label = sprintf("%0.2f", pval)), size = 2.2) +
    scale_fill_manual(values = c("red", "gray"), name = "") +
    scale_y_discrete(limits = rev, position = "right") +
    theme(panel.background = element_blank(),
          strip.background = element_blank(),
          strip.text.x = element_blank(),
          plot.margin = unit(c(0, 0, 0, 0), "cm"),
          panel.border = element_rect(colour = "black", fill = NA),
          axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust = 0.4,
                                     size = 6),
          axis.text.y = element_text(size = 6),
          legend.background = element_blank(),
          legend.text = element_text(size = 8)) +
    labs(x = "",
         y = "") +
    facet_grid(. ~ sample,
               scales = "free_x",
               space = "free")

P1_legend <- cowplot::get_legend(P1)
P1 <- P1 + theme(legend.position = "none")

P2 <-
  lung_percent %>%
  ggplot(., aes(x = x_label, y = lung_perc, fill = sample)) +
    geom_bar(stat = "identity") +
    ylab("Percent\nin lung") +
    xlab("") +
    labs(fill = "Model",
         y = "Percent in\nmetastases(%)") +
    scale_fill_brewer(palette = "Set2") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(size = 6),
          axis.text.y = element_text(size = 6),
          panel.background = element_blank(),
          strip.background = element_blank(),
          strip.text.x = element_blank(),
          plot.margin = unit(c(0, 0, -0.2, 0), "cm"),
          panel.border = element_rect(colour = "black", fill = NA),
          legend.text = element_text(size = 8)) +
    facet_grid(. ~ sample,
               scales = "free_x",
               space = "free")

P2_legend <- cowplot::get_legend(P2)
P2 <- P2 + theme(legend.position = "none")

repfigure_3d <-
    cowplot::plot_grid(cowplot::plot_grid(P2,
                                          P1,
                                          align = "v",
                                          rel_heights = c(1, 6),
                                          ncol = 1),
                       cowplot::plot_grid(P1_legend,
                                          P2_legend,
                                          ncol = 1,
                                          axis = "t"),
                        nrow = 1,
                        rel_widths = c(0.95, 0.05))

repfigure_3d

save(repfigure_3d, file = "Figures/supplement/requests/k7m2_f420/repfigure_3d.RData")

pdf("Figures/supplement/requests/k7m2_f420/repfigure_3d.pdf",
       width = 15,
       height = 8)
repfigure_3d
dev.off()
```

```{r assemblingFig3, eval=T}
load("Figures/supplement/requests/k7m2_f420/repfigure_3b_15.RData")
load("Figures/supplement/requests/k7m2_f420/suppFigure_3c_umaps.RData")
load("Figures/supplement/requests/k7m2_f420/suppFigure_3c_barplots.RData")
load("Figures/supplement/requests/k7m2_f420/repfigure_3d.RData")

layout_matrix <-
    matrix(c(1,  1,  1,  1,  2,  2,  2,  8,  4,  4,  4,  10,
             1,  1,  1,  1,  2,  2,  2,  8,  4,  4,  4,  10,
             1,  1,  1,  1,  3,  3,  3,  9,  5,  5,  5,  11,
             1,  1,  1,  1,  3,  3,  3,  9,  5,  5,  5,  11,
             1,  1,  1,  1,  6,  6,  6,  12, 7,  7,  7,  13,
             1,  1,  1,  1,  6,  6,  6,  12, 7,  7,  7,  13,
             14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
             14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
             14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
             14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
             14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14),
         nrow = 11,
         byrow = TRUE)

pdf("Figures/supplement/requests/k7m2_f420/suppFigure_3_with_mouse.pdf",
    width = 8,
    height = 9)
gridExtra::grid.arrange(grobs = list(repfigure_3b_15,
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$OS17_TL),
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$`143B_TL`),
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$NCHOS2_TL),
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$NCHOS7_TL),
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$K7M2),
                                     patchwork::patchworkGrob(suppfigure_3c_umaps$F420),
                                     suppFigure_3c_barplots$`OS-17`,
                                     suppFigure_3c_barplots$`143B`,
                                     suppFigure_3c_barplots$`NCH-OS-2`,
                                     suppFigure_3c_barplots$`NCH-OS-7`,
                                     suppFigure_3c_barplots$K7M2,
                                     suppFigure_3c_barplots$F420,
                                     repfigure_3d),
                        layout_matrix = layout_matrix)
dev.off()
```

## Patient tumors and cell lines

*Note: this code chunk depends on objects created in OSHetero2021_Supplemental_Figures.Rmd*

```{r patientvpdx}
library(harmony)
# Start from this stopping point from OSHetero2021_Supplemental_Figures.Rmd
comb <- qs::qread("Data/Supplement/PatientData/comb.qs")

DimPlot(comb, reduction = "umap", label = T, repel = T) +
  coord_fixed() +
  ggtitle("GSE152048 composite") +
  theme(legend.position = "none")

DimPlot(comb, reduction = "umap", group.by = "src") +
  coord_fixed() +
  ggtitle("GSE152048 composite") 

DimPlot(comb, reduction = "umap", group.by = "type") +
  coord_fixed() +
  ggtitle("GSE152048 composite") 

# Isolate the tumor cells from mets
# Calculate and show module scores from the paper
ms <- list()
ms$Osteoblastic <- c("RUNX2", "COL1A1", "CDH11", "IBSP")
ms$Chondroblastic <- c("SOX9", "ACAN", "PTH1R")
ms$Osteoclast <- c("ACP5", "CTSK", "MMP9")
ms$Myeloid <- c("CD74", "CD14", "FCGR3A")
ms$TCell <- c("CD3E", "IL7R", "CD8A", "CD4", "NKG7")
ms$NKCell <- c("NKG7", "GNLY")
ms$NKTCell <- c("NKG7", "GNLY", "CD3E")
ms$DCCell <- c("CD1C", "FCER1A", "CLEC9A", "CCR7", "CD14", "CD163")
ms$Fibroblast <- c("DCN", "COL1A1")
ms$Pericyte <- c("RGS5", "ACTA2")
ms$MSC <- c("MME", "THY1", "CXCL12", "SFRP2")
ms$Endothelial <- c("PECAM1", "VWF")
ms$Myoblast <- c("MYL1", "MYLPF")
ms$BCell <- c("MS4A1", "CD19", "JCHAIN")

mod_names <- names(ms)
plots <- list()
#Added min.cutoff to isolate for cells most module like
for(i in 1:length(mod_names)) {
  comb <- AddModuleScore(comb,
                         ms[i],
                         name = names(ms[i]))
  plots[[i]] <- FeaturePlot(comb,
                          features = str_c(names(ms[i]),
                                           "1"),
                          min.cutoff = 2,
                          pt.size = 1,
                          order = T,
                          label = T,
                          cols = c("lightgoldenrod",
                                   "darkred"))
  print(plots[[i]])
}

#look for proliferation
comb<-kill_cc(comb)
DimPlot(comb, reduction = "umap", label = T, repel = T) +
  coord_fixed() +
  ggtitle("G1/S/G2M cell cycle analysis - comb tumor") +
  theme(legend.position = "none")

#Select Clusters which are tumors. Notice how one group is G2M only. 
#Maybe possible that this group is tumor and due to 0inflation there are transcripts which are not active. 
#Could be cycling back to other tumor states
Idents(comb) <- comb$seurat_clusters
DimPlot(comb, reduction = "umap", label = T, repel = T) +
coord_fixed() +
  ggtitle("GSE152048 composite clusters - primaries") +
  theme(legend.position = "none")
#saveRDS(comb,'Data/Supplement/PatientData/comb.rds')
#comb<-readRDS('Data/Supplement/PatientData/comb.rds')

#### Tumor - Met ####

# Separate out the met tumor lesions
met <- subset(comb, subset = type == "Lung Met")

# Re-cluster after subsetting to tumor cells
# Since have not re-processed, the cluster names are the same as those used in primary where clusters are shared between the two
Idents(met) <- met$seurat_clusters
met <- subset(met,
              idents=c(2,
                       3,
                       27,
                       16,
                       15,
                       6,
                       #28,
                       5,
                       20,
                       #29,
                       0,
                       24,
                       14,
                       21,
                       33,
                       #35,
                       13,
                       26,
                       7,
                       9
                       #23
                       ))
DimPlot(met, reduction = "umap", label = T, repel = T) +
coord_fixed() +
  ggtitle("GSE152048 composite clusters - primaries") +
  theme(legend.position = "none")

met <- met %>%
  NormalizeData() %>%
  ScaleData() %>%
  FindVariableFeatures()
met <- met %>%
  RunPCA(features = VariableFeatures(met)) %>%
  RunUMAP(dims=1:20) %>%
  FindNeighbors() %>%
  FindClusters(resolution=0.3)

#Plot results
DimPlot(met,
        pt.size=1,
        reduction="umap",
        label=T)+
  coord_fixed()+
  ggtitle("Reclustering of Met Cells")

#kill_cc adds metadata
met <- kill_cc(met)

DimPlot(met, pt.size=1, reduction="umap", label=T)+
  coord_fixed()+
  ggtitle("Cell Cycle of met Cells")

metplots2 <- list()
#Added min.cutoff to isolate for cells most module like
Idents(met) <- met$seurat_clusters
for(i in 1:length(mod_names)) {
  met <- AddModuleScore(met,
                        ms[i],
                        name = names(ms[i]))
  metplots2[[i]] <- FeaturePlot(met,
                                features = str_c(names(ms[i]),
                                                 "1"),
                                pt.size = 1,
                                order = T,
                                label = T,
                                cols = c("lightgoldenrod",
                                         "darkred"))
  print(metplots2[[i]])
}

#### Met_comb ####

#Now that sure of tumor cells. Take the tumor subsets from "met" and combine across data with harmony
#cluster 12 has markers of epithelial nature(top 10 DEgs: SFN,S100A14,ELF3,LRP2,MLPH,LAMP3)
met_comb<-subset(met,
                 idents = c(0,
                            1,
                            2,
                            6))
met_comb <- met_comb %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F)

met_comb <- RunHarmony(met_comb,
                       group.by.vars = "src")
met_comb <- RunUMAP(met_comb,
                    reduction = "harmony",
                    dims = 1:30)
met_comb <- met_comb %>%
  FindNeighbors(reduction = "harmony",
                dims = 1:30) %>%
  FindClusters(resolution=0.3)

# Save a stopping point - merged and harmony aligned
qs::qsave(met_comb, "Data/Supplement/PatientData/met_comb.qs")
# Start from this stopping point
met_comb <- qs::qread("Data/Supplement/PatientData/met_comb.qs")

DimPlot(met_comb, reduction = "umap", label = T, repel = T) +
  coord_fixed() +
  ggtitle("met Clusters post-integration") +
  theme(legend.position = "none")

DimPlot(met_comb, reduction = "umap", group.by = "src") +
  coord_fixed() +
  ggtitle("met by sample post-integration") 

met_comb<-kill_cc(met_comb)

#Run cell type module as before
metplots3 <- list()
Idents(met_comb)<-met_comb$seurat_clusters
for(i in 1:length(mod_names)) {
  met_comb <- AddModuleScore(met_comb,
                             ms[i], name = names(ms[i]))
  metplots3[[i]] <- FeaturePlot(met_comb,
                                features = str_c(names(ms[i]), "1"),
                                pt.size = 1,
                                order = T,
                                label = T,
                                cols = c("lightgoldenrod",
                                         "darkred"))
  print(metplots3[[i]])
}
```

```{r patientdeg}
tumor_comb <- qs::qread("Data/Supplement/PatientData/tumor_comb.qs")
met_comb <- qs::qread("Data/Supplement/PatientData/met_comb.qs")

all_tumor <- merge(x = tumor_comb,
                   y = met_comb,
                   project = "Patient Tumor Comparison")

all_tumor <- all_tumor %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F)

all_tumor <- RunHarmony(all_tumor,
                       group.by.vars = "src")
all_tumor <- RunUMAP(all_tumor,
                    reduction = "harmony",
                    dims = 1:30)
all_tumor <- all_tumor %>%
  FindNeighbors(reduction = "harmony",
                dims = 1:30) %>%
  FindClusters(resolution=0.3)

# Save a stopping point - merged and harmony aligned
qs::qsave(all_tumor, "Data/Supplement/PatientData/all_tumor.qs")
# Start from this stopping point
all_tumor <- qs::qread("Data/Supplement/PatientData/all_tumor.qs")

set.seed(100)
# DEG list between mets and primary patient tumors
patient_de_genes <- list()
patient_de_genes_dir <- list()
for (tissue_name in c("Primary", "Lung Met")) {
  for (direction in c("up", "down")) {
    Idents(all_tumor) <- all_tumor$type
    patient_de_genes <-
      FindMarkers(all_tumor,
                  ident.1 = "Lung Met",
                  ident.2 = "Primary",
                  ionly.pos = FALSE,
                  min.pct = 0.1)
    
    # Select genes that change in the direction of interest
    if (direction == "up") {
      patient_de_genes_dir[["up"]] <-
        patient_de_genes[patient_de_genes$p_val_adj < 0.05 &
                           patient_de_genes$avg_log2FC > 0, ]
    } else {
      patient_de_genes_dir[["down"]] <-
        patient_de_genes[patient_de_genes$p_val_adj < 0.05 &
                           patient_de_genes$avg_log2FC < 0, ]
    }
  }
}

save(patient_de_genes_dir, file="Data/Supplement/patient_de_genes_dir.RData")
```

## Fig4C, S19B-28B: Met Module scores and Scatter Plots

### Fig4C, S19B-22B Module Score Feature plots

```{r METs5-s6-s7-s8_modules, dependson="s5-s6-s7_processing"}
# Load met samples
met_comb <- qs::qread("Data/Supplement/PatientData/met_comb.qs")

# Set up directory
if(!dir.exists("Figures/supplement/met")) {
  dir.create("Figures/supplement/met")
}

Idents(met_comb) <- met_comb$src
(table(Idents(met_comb)))

#Subset to each sample
s <- c("BC10", "BC17")

met_tumor <- c()
met_tumor[["BC10"]] <- subset(met_comb, subset = src == "BC10")
met_tumor[["BC17"]] <- subset(met_comb, subset = src == "BC17")

met_tumor[["BC10"]] <- met_tumor[["BC10"]] %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F) %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters(resolution = 0.1) %>%
  RunUMAP(dims = 1:20)

met_tumor[["BC17"]] <- met_tumor[["BC17"]] %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA(verbose = F) %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters(resolution = 0.1) %>%
  RunUMAP(dims = 1:20)

for (i in 1:length(met_tumor)){
  print(DimPlot(met_tumor[[i]], pt.size = 1, label = T) +
          coord_fixed() +
          theme(legend.position = "none") +
          ggtitle(str_c(s[i], " basic clustering")))
}

for(i in 1:length(met_tumor)){
  met_tumor[[i]] <- met_tumor[[i]] %>% 
    FindClusters(resolution = 0.7)
  print(DimPlot(met_tumor[[i]], pt.size = 1, label = T) +
          coord_fixed() +
          theme(legend.position = "none") +
          ggtitle(str_c(s[i], " basic clustering")))
}

if (!dir.exists("Figures/supplement/met")) {
  dir.create("Figures/supplement/met")
}
# For Glycolysis
# Apply Glycolysis module score to each of the primary tumor cells
gl <- list()
Glycolysis <- read.table(file = "pathways/Glycolysis.txt",
                         header = FALSE,
                         sep = ",")
Glycolysis <- Glycolysis[, 1]

# for s6
Hypoxia <- read.table(file = "pathways/Hypoxia.txt",
                      header = FALSE,
                      sep = ",")
Hypoxia <- Hypoxia[, 1]

# For EMT
EMT <- read.table(file = "pathways/EMT.txt",
                  header = FALSE,
                  sep = ",")
EMT <- EMT[, 1]

# For TNFA
TNFA <- read.table(file = "pathways/TNFA.txt",
                   header = FALSE,
                   sep = ",")
TNFA <- TNFA[, 1]

gl$Glycolysis <- Glycolysis
gl$Hypoxia <- Hypoxia
gl$EMT <- EMT
gl$TNFA <- TNFA

mod_names <- names(gl)

#### Glycolysis ####
# Glycolysis - Added min.cutoff to isolate for cells most module like
gly <- list()
for (i in 1:length(met_tumor)) {
  met_tumor[[i]] <- AddModuleScore(met_tumor[[i]],
                             gl[1],
                             name = names(gl[1]))
  p <- FeaturePlot(met_tumor[[i]],
                   features = str_c(names(gl[1]),
                                    "1"),
                   min.cutoff = 1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste0(s[i])) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/met/",
                    names(gl[1]),
                    "_FeaturePlot_",
                    s[i],
                    ".pdf"))
  print(p)
  dev.off()
  
  gly[[i]] <- p
}

# Eliminate outlier cell skewing feature plot
df <- as.data.frame(met_tumor[["BC17"]]$Glycolysis1)
colnames(df) <- "gly"
df$cid <- rownames(df)
df[order(-df$gly), ] %>% head()
#       gly                     cid
# 0.2999779 BC17_TTTAGTCAGAGAACCC-1
# 0.2001013 BC17_CTCATGCGTTCAACGT-1
# 0.1995038 BC17_CACATGACATCCTTCG-1
# 0.1954288 BC17_TGAGGTTAGTTGAAGT-1
# 0.1880591 BC17_TTGCATTTCTCCATAT-1
# 0.1797822 BC17_GGAGCAATCACAAGGG-1
odd_cell <- "BC17_TTTAGTCAGAGAACCC-1"
met_tumor[["BC17"]] <-
  met_tumor[["BC17"]][, !colnames(met_tumor[["BC17"]]) %in% odd_cell]

# Re-plot
for (i in 1:length(met_tumor)) {
  met_tumor[[i]] <- AddModuleScore(met_tumor[[i]],
                             gl[1],
                             name = names(gl[1]))
  p <- FeaturePlot(met_tumor[[i]],
                   features = str_c(names(gl[1]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste0(s[i])) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/met/",
                    names(gl[1]),
                    "_FeaturePlot_",
                    s[i],
                    ".pdf"))
  print(p)
  dev.off()
  
  gly[[i]] <- p
  
}

gly_msplots <- gly[[1]] +
  gly[[2]] 

pdf(file = "Figures/supplement/met/gly_msplots.pdf",
    height=10,
    width=10)
gly_msplots
dev.off()

#### Hypoxia ####
# Hypoxia - Added min.cutoff to isolate for cells most module like
hypx <- list()
for (i in 1:length(met_tumor)) {
  met_tumor[[i]] <- AddModuleScore(met_tumor[[i]],
                             gl[2],
                             name = names(gl[2]))
  p <- FeaturePlot(met_tumor[[i]],
                   features = str_c(names(gl[2]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste0(s[i])) +
    theme(aspect.ratio = 1)

  pdf(file = paste0("Figures/supplement/met/",
                    names(gl[2]),
                    "_FeaturePlot_",
                    s[i],
                    ".pdf"))
  print(p)
  dev.off()
  
  hypx[[i]] <- p
}
hypx_msplots <- hypx[[1]] +
  hypx[[2]] 

pdf(file = "Figures/supplement/met/hypx_msplots.pdf",
    height=10,
    width=10)
hypx_msplots
dev.off()

#### EMT ####
# EMT - Added min.cutoff to isolate for cells most module like
emt <- list()
for (i in 1:length(met_tumor)) {
  met_tumor[[i]] <- AddModuleScore(met_tumor[[i]],
                             gl[3],
                             name = names(gl[3]))
  p <- FeaturePlot(met_tumor[[i]],
                   features = str_c(names(gl[3]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste0(s[i])) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/met/",
                    names(gl[3]),
                    "_FeaturePlot_",
                    s[i],
                    ".pdf"))
  print(p)
  dev.off()
  
  emt[[i]] <- p
}

emt_msplots <- emt[[1]] +
  emt[[2]]

pdf(file = "Figures/supplement/met/emt_msplots.pdf",
    height=10,
    width=10)
emt_msplots
dev.off()

#### TNFA ####
# TNFA - Added min.cutoff to isolate for cells most module like
tnf <- list()
for (i in 1:length(met_tumor)) {
  met_tumor[[i]] <- AddModuleScore(met_tumor[[i]],
                             gl[4],
                             name = names(gl[4]))
  p <- FeaturePlot(met_tumor[[i]],
                   features = str_c(names(gl[4]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste0(s[i])) +
    theme(aspect.ratio = 1)

  pdf(file = paste0("Figures/supplement/met/",
                    names(gl[4]),
                    "_FeaturePlot_",
                    s[i],
                    ".pdf"))
  print(p)
  dev.off()
  
  tnf[[i]] <- p
}

tnf_msplots <- tnf[[1]] +
  tnf[[2]]

pdf(file = "Figures/supplement/met/tnf_msplots.pdf",
    height=10,
    width=10)
tnf_msplots
dev.off()

# Note that this object has one outlier cell removed
save(met_tumor, file = "Data/Supplement/PatientData/met_tumor.RData")
```

### S23B-28B: Met Module score scatter plots

```{r METs5-s6-s7-s8_scatter, dependson=c("s5-s6-s7_processing","s5-s6-s7-s8_modules")}
# Define colors
colors <- c("#F8766D",
            "#C49A00",
            "#53B400",
            "#00C094",
            "#00B6EB",
            "#A58AFF",
            "#FB61D7")

s <- c("BC10", "BC17")

#### Glycolysis vs all others ####
combo <- list()
for (i in 1:3) {
  for (j in 1:length(met_tumor)) {
    p <- FeatureScatter(met_tumor[[j]],
                        feature1 = str_c(names(gl[1]),
                                         "1"),
                        
                        feature2 = str_c(names(gl[i+1]),
                                         "1"),
                        group.by = "src",
                        cols = colors[[j]],
                        pt.size = 0.5) +
      theme(aspect.ratio = 1)
    print(p)
    

    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[1]),
                      names(gl[i+1]),
                      "_ScatterPlot_",
                      s[j],
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()
    
    combo[[j]] <- p
  }
  sfig <- combo[[1]] +
    combo[[2]] 

    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[1]),
                      names(gl[i+1]),
                      "_combined",
                      ".pdf"),
        height=10,
        width=10)
    print(sfig)
    dev.off()
}

#### Hypoxia vs Remaining ####
combo <- list()
for (i in 1:2) {
  for (j in 1:length(met_tumor)) {
    p <- FeatureScatter(met_tumor[[j]],
                        feature1 = str_c(names(gl[2]),
                                         "1"),
                        
                        feature2 = str_c(names(gl[i+2]),
                                         "1"),
                        group.by = "src",
                        cols = colors[[j]],
                        pt.size = 0.5) +
      theme(aspect.ratio = 1)
    print(p)

    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[2]),
                      names(gl[i+2]),
                      "_ScatterPlot_",
                      s[j],
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()

    combo[[j]] <- p
  }
    sfig <- combo[[1]] +
      combo[[2]]

    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[2]),
                      names(gl[i+2]),
                      "_combined",
                      ".pdf"),
        height=10,
        width=10)
    print(sfig)
    dev.off()
}

#### Last Two vs Each other ####
combo <- list()

  for (i in 1:length(met_tumor)) {
    p <- FeatureScatter(met_tumor[[i]],
                        feature1 = str_c(names(gl[3]),
                                         "1"),
                        
                        feature2 = str_c(names(gl[4]),
                                         "1"),
                        group.by = "src",
                        cols = colors[[i]],
                        pt.size = 0.5) +
      theme(aspect.ratio = 1)
    print(p)
    
    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[3]),
                      names(gl[4]),
                      "_ScatterPlot_",
                      s[i],
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()
    
    combo[[i]] <- p
  }
    sfig <- combo[[1]] +
      combo[[2]] 
    
    pdf(file = paste0("Figures/supplement/met/",
                      names(gl[3]),
                      names(gl[4]),
                      "_combined",
                      ".pdf"),
        height=10,
        width=10)
    print(sfig)
    dev.off()
```

## S9-18: Module Scores and Scatter Plots for Cell lines and PDXs

### S9-12 Module Scores

```{r ms}
library(rrrSingleCellUtils)
library(patchwork)
library(Seurat)
library(ggplot2)
library(tidyverse)
library(stringr)
# library(future)
library(harmony)
library(devtools)
source("scSeurat.R")
set.seed(888)

# Equal number since feature plots can appear misleading if have uneven numbers
load("Data/OSlist_1500_CCR.RData")

# Apply Glycolysis module score to each of the primary tumor cells
gl <- list()
Glycolysis <- read.table(file = "pathways/Glycolysis.txt",
                         header = FALSE,
                         sep = ",")
Glycolysis <- Glycolysis[, 1]


Hypoxia <- read.table(file = "pathways/Hypoxia.txt",
                      header = FALSE,
                      sep = ",")
Hypoxia <- Hypoxia[, 1]


EMT <- read.table(file = "pathways/EMT.txt",
                  header = FALSE,
                  sep = ",")
EMT <- EMT[, 1]


TNFA <- read.table(file = "pathways/TNFA.txt",
                   header = FALSE,
                   sep = ",")
TNFA <- TNFA[, 1]

gl$Glycolysis <- Glycolysis
gl$Hypoxia <- Hypoxia
gl$EMT <- EMT
gl$TNFA <- TNFA

mod_names <- names(gl)

#### Glycolysis ####
# Glycolysis - Added min.cutoff to isolate for cells most module like
gly <- list()
for (i in 1:length(OS.list_1500)) {
  OS.list_1500[[i]] <- AddModuleScore(OS.list_1500[[i]],
                             gl[1],
                             name = names(gl[1]))
  p <- FeaturePlot(OS.list_1500[[i]],
                   features = str_c(names(gl[1]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste(head(OS.list_1500[[i]]@meta.data$model,1),
                   head(OS.list_1500[[i]]@meta.data$src,1),
                   "-",
                   names(gl[1]))) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/requests/",
                    names(gl[1]),
                    "_FeaturePlot_",
                    paste0(head(OS.list_1500[[i]]@meta.data$model,1),
                           head(OS.list_1500[[i]]@meta.data$src,1),
                           ".pdf")))
      print(p)
      dev.off()
  
  gly[[i]] <- p
  
}

gly_plot <- gly[[1]] +
  gly[[2]] +
  gly[[3]] +
  gly[[4]] +
  gly[[5]] +
  gly[[6]] +
  gly[[7]] +
  gly[[8]] +
  gly[[9]] +
  gly[[10]] +
  gly[[11]] +
  gly[[12]] +
  gly[[13]] +
  plot_layout(ncol = 3)

pdf(file = "Figures/supplement/requests/gly_plot.pdf",
    height=26,
    width=15)
gly_plot
dev.off()

#### Hypoxia ####
# Hypoxia - Added min.cutoff to isolate for cells most module like
hypx <- list()
for (i in 1:length(OS.list_1500)) {
  OS.list_1500[[i]] <- AddModuleScore(OS.list_1500[[i]],
                             gl[2],
                             name = names(gl[2]))
  p <- FeaturePlot(OS.list_1500[[i]],
                   features = str_c(names(gl[2]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste(head(OS.list_1500[[i]]@meta.data$model,1),
                   head(OS.list_1500[[i]]@meta.data$src,1),
                   "-",
                   names(gl[2]))) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/requests/",
                    names(gl[2]),
                    "_FeaturePlot_",
                    paste0(head(OS.list_1500[[i]]@meta.data$model,1),
                           head(OS.list_1500[[i]]@meta.data$src,1),
                           ".pdf")))
  print(p)
  dev.off()
  
  hypx[[i]] <- p
}
hypx_plot <- hypx[[1]] +
  hypx[[2]] +
  hypx[[3]] +
  hypx[[4]] +
  hypx[[5]] +
  hypx[[6]] +
  hypx[[7]] +
  hypx[[8]] +
  hypx[[9]] +
  hypx[[10]] +
  hypx[[11]] +
  hypx[[12]] +
  hypx[[13]] +
  plot_layout(ncol = 3)

pdf(file = "Figures/supplement/requests/hypx_plot.pdf",
    height=26,
    width=15)
hypx_plot
dev.off()

#### EMT ####
# EMT - Added min.cutoff to isolate for cells most module like
emt <- list()
for (i in 1:length(OS.list_1500)) {
  OS.list_1500[[i]] <- AddModuleScore(OS.list_1500[[i]],
                             gl[3],
                             name = names(gl[3]))
  p <- FeaturePlot(OS.list_1500[[i]],
                   features = str_c(names(gl[3]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste(head(OS.list_1500[[i]]@meta.data$model,1),
                   head(OS.list_1500[[i]]@meta.data$src,1),
                   "-",
                   names(gl[3]))) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/requests/",
                    names(gl[3]),
                    "_FeaturePlot_",
                    paste0(head(OS.list_1500[[i]]@meta.data$model,1),
                           head(OS.list_1500[[i]]@meta.data$src,1)),
                    ".pdf"))
  print(p)
  dev.off()
  
  emt[[i]] <- p
}

emt_plot <- emt[[1]] +
  emt[[2]] +
  emt[[3]] +
  emt[[4]] +
  emt[[5]] +
  emt[[6]] +
  emt[[7]] +
  emt[[8]] +
  emt[[9]] +
  emt[[10]] +
  emt[[11]] +
  emt[[12]] +
  emt[[13]] +
  plot_layout(ncol = 3)

pdf(file = "Figures/supplement/requests/emt_plot.pdf",
    height=26,
    width=15)
emt_plot
dev.off()

#### TNFA ####
# TNFA - Added min.cutoff to isolate for cells most module like
tnf <- list()
for (i in 1:length(OS.list_1500)) {
  OS.list_1500[[i]] <- AddModuleScore(OS.list_1500[[i]],
                             gl[4],
                             name = names(gl[4]))
  p <- FeaturePlot(OS.list_1500[[i]],
                   features = str_c(names(gl[4]),
                                    "1"),
                   min.cutoff = 0.1,
                   pt.size = 1,
                   order = T,
                   cols = c("lightgoldenrod",
                            "darkred")) +
    ggtitle(paste(head(OS.list_1500[[i]]@meta.data$model,1),
                   head(OS.list_1500[[i]]@meta.data$src,1),
                   "-",
                   names(gl[4]))) +
    theme(aspect.ratio = 1)
  
  pdf(file = paste0("Figures/supplement/requests/",
                    names(gl[4]),
                    "_FeaturePlot_",
                    paste0(head(OS.list_1500[[i]]@meta.data$model,1),
                           head(OS.list_1500[[i]]@meta.data$src,1)),
                    ".pdf"))
  print(p)
  dev.off()
  
  tnf[[i]] <- p
}

tnf_plot <- tnf[[1]] +
  tnf[[2]] +
  tnf[[3]] +
  tnf[[4]] +
  tnf[[5]] +
  tnf[[6]] +
  tnf[[7]] +
  tnf[[8]] +
  tnf[[9]] +
  tnf[[10]] +
  tnf[[11]] +
  tnf[[12]] +
  tnf[[13]] +
  plot_layout(ncol = 3)

pdf(file = "Figures/supplement/requests/tnf_plot.pdf",
    height=26,
    width=15)
tnf_plot
dev.off()

save(OS.list_1500, file = "Data/Supplement/PatientData/OS.list_1500_ms.RData")
```

## S13-18 Scatter Plots

```{r cellpdx_scatter}
load("Data/Supplement/PatientData/OS.list_1500_ms.RData")

# Define colors
# library(scales)
# print(hue_pal()(13))
colors <- c("#F8766D",
            "#DE8C00",
            "#B79F00",
            "#7CAE00",
            "#00BA38",
            "#00C08B",
            "#00BFC4",
            "#00B4F0",
            "#619CFF",
            "#C77CFF",
            "#F564E3",
            "#FF64B0",
            "#FF65AC")


if (!dir.exists("Figures/supplement/requests/cellms")) {
  dir.create("Figures/supplement/requests/cellms")
}
#### Glycolysis vs all others ####
combo <- list()
for (j in 1:3) {
  for (i in 1:length(OS.list_1500)) {
    p <- FeatureScatter(OS.list_1500[[i]],
                        feature1 = str_c(names(gl[1]),
                                         "1"),
                        
                        feature2 = str_c(names(gl[j+1]),
                                         "1"),
                        group.by = "src",
                        cols = colors[[i]],
                        pt.size = 0.5) +
      theme(aspect.ratio = 1)
    print(p)
    
    pdf(file = paste0("Figures/supplement/requests/cellms/",
                      names(gl[1]),
                      names(gl[j+1]),
                      "_ScatterPlot_",
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()
    
    combo[[i]] <- p
  }
  sfig <- combo[[1]] +
    combo[[2]] +
    combo[[3]] +
    combo[[4]] +
    combo[[5]] +
    combo[[6]] +
    combo[[7]] +
    combo[[8]] +
    combo[[9]] +
    combo[[10]] +
    combo[[11]] +
    combo[[12]] +
    combo[[13]] +
    plot_layout(ncol = 3)
  
  pdf(file = paste0("Figures/supplement/requests/cellms/cellscatter_",
                    j+5,
                    ".pdf"),
      height=13,
      width=15)
  print(sfig)
  dev.off()
}

#### Hypoxia vs Remaining ####
combo <- list()
for (j in 1:2) {
  for (i in 1:length(OS.list_1500)) {
  p <- FeatureScatter(OS.list_1500[[i]],
                      feature1 = str_c(names(gl[2]),
                                       "1"),
                      
                      feature2 = str_c(names(gl[j+2]),
                                       "1"),
                      group.by = "src",
                      cols = colors[[i]],
                      pt.size = 0.5) +
    theme(aspect.ratio = 1)
  print(p)
    
    pdf(file = paste0("Figures/supplement/requests/cellms/",
                      names(gl[2]),
                      names(gl[j+2]),
                      "_ScatterPlot_",
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()
    
    combo[[i]] <- p
  }
  sfig <- combo[[1]] +
    combo[[2]] +
    combo[[3]] +
    combo[[4]] +
    combo[[5]] +
    combo[[6]] +
    combo[[7]] +
    combo[[8]] +
    combo[[9]] +
    combo[[10]] +
    combo[[11]] +
    combo[[12]] +
    combo[[13]] +
    plot_layout(ncol = 3)
  
  pdf(file = paste0("Figures/supplement/requests/cellms/cellscatter_",
                    j+6,
                    ".pdf"),
      height=13,
      width=15)
  print(sfig)
  dev.off()
}

#### Last Two vs Each other ####
combo <- list()

  for (i in 1:length(OS.list_1500)) {
    p <- FeatureScatter(OS.list_1500[[i]],
                        feature1 = str_c(names(gl[3]),
                                         "1"),
                        
                        feature2 = str_c(names(gl[4]),
                                         "1"),
                        group.by = "src",
                        cols = colors[[i]],
                        pt.size = 0.5) +
      theme(aspect.ratio = 1)
    print(p)
    
    pdf(file = paste0("Figures/supplement/requests/cellms/",
                      names(gl[3]),
                      names(gl[4]),
                      "_ScatterPlot_",
                      ".pdf"),
        height = 5,
        width = 5)
    print(p)
    dev.off()
    
    combo[[i]] <- p
  }
sfig <- sfig <- combo[[1]] +
  combo[[2]] +
  combo[[3]] +
  combo[[4]] +
  combo[[5]] +
  combo[[6]] +
  combo[[7]] +
  combo[[8]] +
  combo[[9]] +
  combo[[10]] +
  combo[[11]] +
  combo[[12]] +
  combo[[13]] +
  plot_layout(ncol = 3)

pdf(file = paste0("Figures/supplement/requests/cellms/cellscatter_9.pdf"),
    height=13,
    width=15)
print(sfig)
dev.off()
```

# S31 - Cell Line Imaging 
(no coding needed)

```{r}
sessionInfo()
```