Early_vs_Late_Analysis.Rmd

---
title: "Data Analysis for Early and Late Lung-colonizing Tumors"
author: "Emily Franz, Ryan Roberts, Yogesh Budhathoki and Matt Cannon"
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
)
```

```{r lib}
# Load packages
library(Seurat)
library(tidyverse)
library(ggplot2)
library(rrrSingleCellUtils)
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(patchwork)
library(msigdbr)

# Set random generator seed to facilitate reproducibility
set.seed(888)
```

## Helper functions
```{r helper_functions}
theme_robert <- function(axis_font_size = 5,
                         axis_title_font_size = 8,
                         facet_font_size = 5,
                         font_name = "Arial",
                         legend_key_size = "0.2",
                         legend_text_font_size = 5,
                         legend_title_font_size = 6,
                         subtitle_font_size = 6,
                         title_font_size = 10) {
    if (font_name == "Arial" &
        !"Arial" %in% names(grDevices::pdfFonts())) {
        #    !"Arial" %in% sysfonts::font_families()) {
        message("Adding Arial font to system")
        tryCatch(
            {
                extrafont::font_import(paths = "/home/gdrobertslab/lab/Tools/fonts/Arial/",
                                       prompt = FALSE)
                extrafont::loadfonts(quiet = TRUE)
            },
                error = function(e) {
                    message("Arial font files not found in /home/gdrobertslab/lab/Tools/fonts/Arial/")
                    message("Either use a different font or install Arial")
                    print(e)
            }
        )
    }
    list(ggplot2::theme(
            axis.text = ggplot2::element_text(size = axis_font_size),
            axis.title = ggplot2::element_text(size = axis_title_font_size),
            plot.title = ggplot2::element_text(size = title_font_size,
                                               face = "bold",
                                               hjust = 0.5),
            strip.text = ggplot2::element_text(size = facet_font_size,
                                               face = "bold"),
            plot.subtitle = ggplot2::element_text(hjust = 0.5,
                                                  size = subtitle_font_size),
            legend.text = ggplot2::element_text(size = legend_text_font_size),
            legend.title = ggplot2::element_text(size = legend_title_font_size,
                                                 face = "bold"),
            strip.background = ggplot2::element_blank(),
            legend.key.size = ggplot2::unit(0.2, "cm"),
            legend.position = "top"),
    ggplot2::scale_color_manual(values = plot_cols))
}
```

## Make up directory structure
```{bash mkdirs, eval=FALSE}
for directoryName in \
  output \
  output/figures \
  output/rdata \
  output/de \
  output/gsea
do
  if [ ! -d ${directoryName} ]
  then
    mkdir ${directoryName}
  fi
done
```

# Load and Process Datasets
Datasets are from tumors collected from mouse lungs at early (day 14 post injection of tumor cells) and late (day 35 post injection of tumor cells) timepoints and at the time of injection. mCherry was used to label the tumor cells as well as the surrounding niche cells.
```{r load}
sample_list <- read_tsv("misc/sample_cutoffs.txt", show_col_types = FALSE)

sobj_list <- parallel::mclapply(seq_len(nrow(sample_list)), function(i) {
    message(paste("Starting sample", i))

    species_pattern <- sample_list$species_pattern[i]
    obj_name <- sample_list$obj_name[i]

    sobj <-
        tenx_load_qc(paste0("/home/gdrobertslab/lab/Counts/",
                            sample_list$sample_name[i],
                            "/filtered_feature_bc_matrix"),
                     violin_plot = FALSE,
                     species_pattern = species_pattern)

    # Add metadata to dataset
    for (colname in colnames(sample_list)) {
        sobj[[colname]] <- sample_list[[colname]][i]
    }

    cutoff_tibble <-
        tribble(~feature, ~min_val, ~max_val,
                "nCount_RNA",
                sample_list$ncount_rna_min[i],
                sample_list$ncount_rna_max[i],

                "percent.mt",
                0,
                sample_list$mt_max[i])

    plotted <-
        feature_hist(sobj,
                     features = c("nCount_RNA", "percent.mt"),
                     cutoff_table = cutoff_tibble)

    ggsave(paste0("output/figures/feature_hist_",
                  obj_name,
                  ".png"),
           width = 10,
           height = 10,
           plot = plotted)

    sobj <-
        sobj %>%
        subset(nCount_RNA >= sample_list$ncount_rna_min[i] &
               nCount_RNA <= sample_list$ncount_rna_max[i] &
               percent.mt <= sample_list$mt_max[i]) %>%
        process_seurat()

    message(paste("Sample", i, "completed"))
    return(sobj)
}, mc.cores = parallelly::availableCores())

names(sobj_list) <- sample_list$obj_name

qs::qsave(sobj_list,
          file = "output/rdata/sobj_list.qs")
```

## UMAP with all samples
```{r umap_alldependson='load', cache.vars=''}
sobj_list <- qs::qread("output/rdata/sobj_list.qs")
all_samples <-
    merge(sobj_list[[1]],
          y = sobj_list[2:length(sobj_list)],
          add.cell.ids = names(sobj_list)) %>%
    process_seurat(resolution = 0.2)

plot_name <-
    DimPlot(all_samples,
            group.by = "sample_description",
            split.by = "sample_description",
            shuffle = TRUE,
            ncol = 2,
            cols = plot_cols[seq_along(sobj_list)],
            pt.size = 0.01,
            order = FALSE) +
    theme_robert() +
    labs(title = NULL)

plot_name

ggsave("output/figures/umap_all_samples.pdf",
       width = 4,
       height = 6,
       plot = plot_name)
```

## Human Cell Only Processing
```{r human-processing, dependson="load"}
# Human
# Normalize the number of cells per sample
n_cells <- c()

for (i in seq_along(sobj_list)) {
    n_cells[i] <- length(colnames(sobj_list[[i]]))
}

names(n_cells) <- names(sobj_list)

# Merge day 14 and day 28 Seurat objects together
tumor <-
    merge(sobj_list[["d14_pos_h"]],
          y = sobj_list[c("d28_pos_h", "tibia_tumor_h")],
          add.cell.ids = c("d14", "d28", "tibia"),
          project = "Early vs Late OS Mets") %>%
    NormalizeData() %>%
    ScaleData() %>%
    FindVariableFeatures() %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    FindClusters(resolution = 0.15, k.param = 30) %>%
    RunUMAP(dims = 1:30, n.neighbors = 30L, metric = "euclidean")

tumor$group <-
    factor(str_replace_all(tumor$group, "_", " "),
           levels = c("Primary tumor",
                      "Early metastatic tumor",
                      "Midpoint metastatic tumor",
                      "Late metastatic tumor"))

# Save human mCherry+ merged object
qs::qsave(tumor, file = "output/rdata/human_tumor.qs")

# Use timepoint naming (Day 14/Day 28) for UMAP plotting
umap_timepoint <-
    r_dim_plot(tumor,
        title = "Tumor at metastatic timepoints",
        group.by = "timepoint",
        label = TRUE)
umap_timepoint

pdf("output/figures/human_umap_timepoint.pdf",
    width = 5,
    height = 5)
umap_timepoint
dev.off()

# Use basic naming (early/late) for UMAP plotting
umap_group <-
    r_dim_plot(tumor,
        title = "Tumor at metastatic timepoints",
        group.by = "group",
        label = TRUE,
        repel = TRUE)
umap_group

pdf("output/figures/human_umap_group.pdf",
    width = 5,
    height = 5)
umap_group
dev.off()

# Reorder the group metadata variable
tumor$group[tumor$group == "Early metastatic tumor"] <- "Early"
tumor$group[tumor$group == "Midpoint metastatic tumor"] <- "Established"
tumor$group[tumor$group == "Primary tumor"] <- "Primary"
tumor$group <- factor(tumor$group, levels = c("Primary", "Early", "Established"))

umap_group_split <-
    r_dim_plot(tumor,
        title = "Tumor at metastatic timepoints",
        split.by = "group",
        label = TRUE,
        repel = TRUE)
umap_group_split

```

## Human Cells Only: Pre-cell cycle regression
```{r human_geneplots, dependson=c("load","human-processing")}
# Dot plot of expression of genes of interest by group (early/late)
gene_list <-
    c("BCL2",
      "BCL2L1",
      "BCL2L2",
      "MCL1")

dot_plot <-
    DotPlot(tumor,
            features = gene_list,
            group.by = "group",
            scale = FALSE,
            cols = c("lightgoldenrod", "darkred")) +
  RotatedAxis() +
  coord_fixed() +
  scale_y_discrete(limits = rev) +
  labs(y = NULL,
       x = NULL)

# Save as pdf
pdf("output/figures/early_late_dot_plot.pdf",
    width = 7,
    height = 4)
dot_plot
dev.off()

expr_data <-
    GetAssayData(tumor) %>%
    as.matrix() %>%
    as.data.frame() %>%
    rownames_to_column("gene") %>%
    filter(gene %in% gene_list) %>%
    as_tibble() %>%
    pivot_longer(-gene,
                 names_to = "cell",
                 values_to = "expression") %>%
    mutate(sample = str_remove(cell, "_.*")) %>%
    left_join(tumor@meta.data %>%
                rownames_to_column("cell") %>%
                mutate(sample = str_remove(cell, "_.*")) %>%
                select(sample, group) %>%
                distinct()) %>%
    group_by(group, gene) %>%
    summarize(mean_expr = mean(expression),
              pct_expr = sum(expression > 0) / n(),
              .groups = "drop") %>%
    mutate(gene_group = if_else(gene == "MCL1",
                           "MCL1",
                           "BCL2 family"),
           group = factor(group,
                          levels = c("Primary tumor",
                                     "Early metastatic tumor",
                                     "Midpoint metastatic tumor",
                                     "Late metastatic tumor")))

bcl_dotplot <- function(data) {
    plot_name <-
        ggplot(data,
               aes(x = gene,
                   y = group,
                   size = pct_expr,
                   color = mean_expr)) +
        geom_point() +
        theme_robert() +
        scale_color_gradient(low = "lightgoldenrod",
                            high = "darkred",
                            name = "Mean expr") +
        scale_size_continuous(limits = c(0, 1),
                              name = "% cells expr") +
        theme(panel.grid.major = element_blank(),
              legend.box = "vertical",
              legend.key.width = unit(15, "points")) +
        labs(x = "",
             y = "") +
        scale_y_discrete(limits = rev)

    return(plot_name)
}

split_dotplot <-
    cowplot::plot_grid(bcl_dotplot(expr_data),
                       bcl_dotplot(subset(expr_data, gene_group == "MCL1")) +
                            theme(axis.text.y = element_blank()),
                    align = "h",
                    rel_widths = c(1, 0.5))

# Save as pdf
pdf("output/figures/early_late_split_dot_plot.pdf",
    width = 6,
    height = 4)
split_dotplot
dev.off()

expr_data <-
    GetAssayData(tumor) %>%
    as.matrix() %>%
    as.data.frame() %>%
    rownames_to_column("gene") %>%
    filter(gene %in% gene_list) %>%
    as_tibble() %>%
    pivot_longer(-gene,
                 names_to = "cell",
                 values_to = "expression") %>%
    mutate(sample = str_remove(cell, "_.*")) %>%
    left_join(tumor@meta.data %>%
                rownames_to_column("cell") %>%
                mutate(sample = str_remove(cell, "_.*")) %>%
                select(sample, group) %>%
                distinct()) %>%
    group_by(group, gene) %>%
    summarize(mean_expr = mean(expression),
              pct_expr = sum(expression > 0) / n(),
              .groups = "drop") %>%
    mutate(gene_group = if_else(gene == "MCL1",
                           "MCL1",
                           "BCL2 family"),
           group = factor(group,
                          levels = c("Primary tumor",
                                     "Early metastatic tumor",
                                     "Midpoint metastatic tumor",
                                     "Late metastatic tumor")))

bcl_dotplot <- function(data) {
    plot_name <-
        ggplot(data,
               aes(x = gene,
                   y = group,
                   size = pct_expr,
                   color = mean_expr)) +
        geom_point() +
        theme_robert() +
        scale_color_gradient(low = "lightgoldenrod",
                            high = "darkred",
                            name = "Mean expr") +
        scale_size_continuous(limits = c(0, 1),
                              name = "% cells expr") +
        theme(panel.grid.major = element_blank(),
              legend.box = "vertical",
              legend.key.width = unit(15, "points")) +
        labs(x = "",
             y = "") +
        scale_y_discrete(limits = rev)

    return(plot_name)
}

split_dotplot <-
    cowplot::plot_grid(bcl_dotplot(expr_data),
                       bcl_dotplot(subset(expr_data, gene_group == "MCL1")) +
                            theme(axis.text.y = element_blank()),
                    align = "h",
                    rel_widths = c(1, 0.5))

# Save as pdf
pdf("output/figures/early_late_split_dot_plot.pdf",
    width = 6,
    height = 4)
split_dotplot
dev.off()

dot_plot <-
    DotPlot(tumor,
            features = gene_list[1:3],
            group.by = "group",
            scale = FALSE,
            cols = c("lightgoldenrod", "darkred")) +
  RotatedAxis() +
  coord_fixed() +
  scale_y_discrete(limits = rev) +
  labs(y = NULL,
       x = NULL) +
  DotPlot(tumor,
            features = gene_list[4],
            group.by = "group",
            scale = FALSE,
            cols = c("lightgoldenrod", "darkred")) +
  RotatedAxis() +
  coord_fixed() +
  scale_y_discrete(limits = rev) +
  labs(y = NULL,
       x = NULL)

# Save as pdf
pdf("output/figures/early_late_dot_plot.pdf",
    width = 7,
    height = 4)
dot_plot
dev.off()


# Violin plot of expression of genes of interest by group (early/late)
vln_plot <-
    VlnPlot(tumor,
            features = gene_list,
            group.by = "group") +
    RotatedAxis()

# Save as pdf
pdf("output/figures/early_late_vln_plot.pdf",
    width = 6,
    height = 5)
vln_plot
dev.off()

# Feature plot (UMAP) of expression of genes of interest by group (early/late)
feature_plot <-
    FeaturePlot(tumor,
                features = gene_list,
                split.by = "group",
                order = TRUE)

# Save as pdf
pdf("output/figures/early_late_feature_plot.pdf",
    width = 15,
    height = 18)
feature_plot
dev.off()
```

## Human Cells Only: Post-cell cycle regression

Regress out cell cycle genes so that the dataset does not cluster on the UMAP
by cell cycle gene expression.

```{r cc, dependson=c("load","human-processing")}
# Calculate cell cycle phase and view its effects on clustering
# Note: default for kill_cc() is not performing cell cycle regression

# Save as pdf
pdf("output/figures/early_late_tumor_pre-kill_cc.pdf",
    width = 5,
    height = 5)
kill_cc(tumor)
dev.off()

# Calculate cell cycle phase and remove its effects on clustering
tumor_cc <-
    kill_cc(tumor,
            cc_regress = "Y") %>%
    FindNeighbors(dims = 1:20) %>%
    FindClusters(resolution = 0.15) %>%
    RunUMAP(dims = 1:20)

# Plot tumor UMAP after cell cycle genes are regressed out
DimPlot(tumor_cc,
        group.by = "Phase",
        label = TRUE,
        pt.size = 1) +
    theme_robert()

# Save as pdf
pdf("output/figures/early_late_post-kill_cc.pdf",
    width = 5,
    height = 5)
DimPlot(tumor_cc,
        group.by = "Phase",
        label = TRUE,
        pt.size = 1) +
    theme_robert()
dev.off()

# Plot tumor UMAP by timepoint post cell cycle regression
umap_cc <-
    DimPlot(tumor_cc,
            pt.size = 1,
            label = TRUE,
            group.by = "timepoint") +
    theme_robert() +
    theme(aspect.ratio = 1)
umap_cc

# Save as pdf
pdf("output/figures/early_late_umap_cc_bytrtment.pdf",
    width = 5,
    height = 5)
umap_cc
dev.off()

# Save the cell cycle regressed tumor object
save(tumor_cc, file = "output/rdata/early_late_tumor_cc.RData")
```

Repeat plots of gene expression post cell cycle regression.
Add S0131 to dotplot
```{r human_geneplots_postcc, dependson=c("load","human-processing","cc")}
# Dot plot of expression of genes of interest by group (early/late)
dot_plot <-
    DotPlot(tumor_cc,
            features = gene_list,
            group.by = "group",
            scale = FALSE,
            cols = c("lightgoldenrod", "darkred")) +
  RotatedAxis() +
  coord_fixed()

# Save as pdf
pdf("output/figures/early_late_dot_plot_postccr.pdf",
    width = 7,
    height = 4)
dot_plot
dev.off()

# Violin plot of expression of genes of interest by group (early/late)
vln_plot <-
    VlnPlot(tumor_cc,
            features = gene_list,
            group.by = "group") +
  RotatedAxis()

# Save as pdf
pdf("output/figures/early_late_vln_plot_postccr.pdf",
    width = 6,
    height = 5)
vln_plot
dev.off()

# Feature plot (UMAP) of expression of genes of interest by group (early/late)
feature_plot <-
    FeaturePlot(tumor_cc,
                features = gene_list,
                split.by = "group",
                order = TRUE) +
  theme(aspect.ratio = 1)

# Save as pdf
pdf("output/figures/early_late_feature_plot_postccr.pdf",
    width = 6,
    height = 12)
feature_plot
dev.off()
```

## Growth factor expression in normal vs lungs with mets
Also plot each time point separately as a time-course
```{r deGrowthFactorMets, dependson='load'}
sobj_list <- qs::qread("output/rdata/sobj_list.qs")

sobj_list$S0006_mouse$group_name <- "Mouse control lung"

for (sample_name in c("d14_pos_m", "d28_pos_m", "d49_pos_m")) {
    sobj_list[[sample_name]]$group_name <- "Stromal combined timepoints"
}

combined_data <-
    merge(x = sobj_list$S0006_mouse,
          y = list(sobj_list[["d14_pos_m"]],
                   sobj_list[["d28_pos_m"]],
                   sobj_list[["d49_pos_m"]])) %>%
    process_seurat(umap_metric = "cosine",
                   run_umap_dims = 1:10,
                   umap_n_neighbors = 30,
                   neighbor_k_param = 20,
                   resolution = 1)

combined_data$group_name <-
    factor(combined_data$group_name,
           levels = c("Mouse control lung",
                      "Stromal combined timepoints"))

combined_data$sample_description <-
    factor(combined_data$sample_description,
           levels = c("SCID control lung",
                      "OS17 mets - 2wk",
                      "OS17 mets - 4wk",
                      "OS17 mets - 7wk"))

combined_data$sample_name <-
    combined_data$sample_name %>%
    str_remove("/outs") %>%
    factor(., levels = c("S0006-WTLung",
                         "S0082",
                         "S0084",
                         "S0093"))

combined_data$group <-
    factor(combined_data$group,
           levels = c("Mouse_control_lung",
                      "Early_mouse_stroma",
                      "Midpoint_mouse_stroma",
                      "Late_mouse_stroma"))
```

## Code to try to assess stromal cell type assignment
To assess the stromal cell type assignment, we will use the SingleR package to assign cell types from a couple of different cell references. We then use marker genes taken from several different papers to assess the cell type assignments. We then use a combination of these to assign cell types for each cluster.
```{r stromal_cell_type_assignment, eval=FALSE}
# Load the already-processed mouse lung reference
ref_path <- "/home/gdrobertslab/lab/GenRef/sc_ref_datasets/mouse"
load(paste0(ref_path, "/GSE151974/mouse_lung_ref.RData"))
mouse_lung_ref$label.fine <- mouse_lung_ref$cell_type

# Add it to immune and whole-mouse references
ref <- list(
    immune = celldex::ImmGenData(),
    whole_mouse = celldex::MouseRNAseqData(),
    lung = as.SingleCellExperiment(mouse_lung_ref)
)

# ref_slots <- c("label.fine", "label.main", "cell_type")
# names(ref_slots) <- names(ref)

# Cluster the combined data to a finer resolution
combined_data <- FindClusters(combined_data, resolution = 1)

# Make cell assignments by cluster for each of the references separately
cell_assign <- lapply(ref, function(x) {
    a <- SingleR::SingleR(as.SingleCellExperiment(combined_data),
        ref = x,
        labels = x$label.fine,
        clusters = combined_data$seurat_clusters)
    a <- tibble(cluster = rownames(a), type = a$labels)
    return(a)
})
# Combine the assignments into a table
cell_assign <- left_join(cell_assign[[1]], cell_assign[[2]], by = "cluster") %>%
    left_join(., cell_assign[[3]], by = "cluster") %>%
    as_tibble()
colnames(cell_assign) <- c("cluster", "immune", "whole_mouse", "lung")
cell_assign
# Add the tables that Yogesh assembled from PMID 35690521
mac_subsets <-
    readxl::read_excel("misc/Finalmacrophage_and_monocyte_subsets.xlsx") %>%
    filter(Species == "Mouse" & Assay == "RNA")
mac_subsets$Major[mac_subsets$Major == 0] <- "Minor"
mac_subsets$Major[mac_subsets$Major == 1] <- "Major"
mac_subsets$Annotation <- paste0(mac_subsets$Annotation, "_", mac_subsets$Major)
mac_subsets <- select(mac_subsets, Annotation, ID)
mac_subset_list <- lapply(unique(mac_subsets$Annotation), function(s) {
    return(mac_subsets$ID[mac_subsets$Annotation == s])
})
names(mac_subset_list) <- unique(mac_subsets$Annotation)
# Examine the macrophage clusters more closely
# Look at markers published in PMID 35732118, 35812401, 32302573, 35690521
murine_mac_marks <- c(mac_subset_list, list(
    inflam_mono = c("Ms4a4c", "Irf7", "Plac8", "Ly6c2"),
    Cd36_Spp1_mac = c("Clec4d", "Hmox1", "Cd36", "Spp1", "Arg1"),
    C1qa_mac = c("C1qa", "Fcgr3", "Ccr2", "Tgfbi"),
    MhcHigh = c("Cd74", "H2-Aa", "Plbd1", "Ciita", "Il1r2", "Ifitm1", "Napsa",
        "H2-Oa", "H2-DMb2", "Cd209a", "Plet1"),
    Batf3_Dc = c("Flt3", "Batf3", "Zbtb46", "Ccr7", "Fscn1", "Ccl22", "Nfkb2"),
    Cd14Mono = c("S100a8", "S100a9", "S100a12", "Vcan", "Cd36", "Cd14"),
    Cd16Mono = c("Tcf7l2", "Mtss1", "Serpina1", "Rhoc", "Ms4a7", "Lilra1",
        "Ifitm2", "Fcgr3a", "Siglec10", "Cx3cr1", "Lilrb1"),
    Nlrp3Mac = c("Nlrp3", "Ereg", "Itgax", "Il1b", "Vegfa", "Ccl3"),
    C1qcTAM = c("Apoe", "C1qa", "C1qb", "C1qc", "Cd81", "Trem2", "Acp5",
        "Apoc1", "Prdm1", "Csf1r", "Slamf8"),
    Spp1TAM = c("Cstb", "Spp1", "Fn1", "Abl2", "Pparg", "Sdc4", "Rgcc",
        "Adm", "Marco"),
    Mast = c("Ms4a2", "Cpa3", "Ccl3", "Ccl6", "Tpsb2", "Tpsd1", "Hdc",
        "Slc18a2", "Gata2", "Kit", "Il1r1", "Hpgds", "Vwa5a"),
    pDC = c("Gzmb", "Tcf4", "Jchain", "Mzb1", "Irf7", "Il3ra", "Slc15a4"),
    Cd1c_cDC2 = c("Axl", "Fcer1a", "Cd1c", "Fcgr2b", "Cd1e", "Cd1a"),
    Batf3_cDC1 = c("Cst7", "Dapp1", "Flt3", "C1orf54", "Ccr7", "Ido1"),
    C1qaMac = c("C1qa", "Marco", "Apoe", "Cxcl10", "Cd163", "Cd68"),
    Thbs1Mac = c("Thbs1", "Fcn1", "Vcan", "S100a12", "Nlrp3", "Il1b"),
    MafbMac = c("Klf2", "Klf4", "Nr4a1", "Nr4a2", "Hsph1", "Lyve1", "Vegfa"),
    Mki67Mac = c("Stmn1", "Top2a", "Mki67", "Mcm4", "Mcm5", "Mcm7"),
    S100a8Mac = c("S100a8", "Ctsa", "Ctsd", "S100a2", "S100a6"),
    Spp1Mac = c("Spp1", "Tmsb4x", "Tmsb10", "Ctsb"),
    Neutro = c("S011a9", "S100a8", "Csfr", "Clec4d"),
    cDC2_Cd209a = c("Cd209a"),
    cDC2_Itgax = c("Itgax"),
    cDC1_Clec9a = c("Cd207", "Clec9a", "Batf3"),
    cDC1_Ccl22 = c("Ccl22", "Ccr7", "Fscn1", "Ccl5"),
    Mono_Ly6c2 = c("Cxcl10", "Ly6c1", "Fcgr1", "Ly6c2"),
    Mono_Nr4a1 = c("Fn1", "Nr4a1"),
    Mono_Itga1 = c("Itga1", "Cd300a", "Clec4a1"),
    Macro_Maf = c("Maf", "Cx3cr1", "Ccl2"),
    Macro_Mgl2 = c("Mgl2", "Clec4b1"),
    Macro_Vegfa = c("Vegfa", "Hilpda", "Hmox1", "Spp1"))
)
pdf("output/figures/mac_subset_marks.pdf",
    width = 12,
    height = 12)
gridExtra::grid.table(cell_assign)
DimPlot(combined_data,
    pt.size = 1,
    label = TRUE) +
    coord_fixed()
DimPlot(combined_data,
    pt.size = 1,
    split.by = "timepoint") +
    coord_fixed()
for (i in names(murine_mac_marks)) {
    p <- r_feature_plot(combined_data,
        title = i,
        features = murine_mac_marks[[i]],
        ncol = 3)
    print(p)
}
dev.off()
r_feature_plot(combined_data,
        title = "Irf4",
        features = "Irf4")
```

### Final cell type assignment

The following was run previously to establish the cell types within the stromal populations based on the preceding analysis:
cell_assignments <-
    read_csv("misc/macrophage_cell_type.csv",
             show_col_types = FALSE) %>%
    pull(label.main, name = cluster)

combined_data <-
    RenameIdents(combined_data, cell_assignments)
combined_data$cell_type_final <- Idents(combined_data)

cell_type_assignments <- combined_data$cell_type_final %>%
    as.data.frame() %>%
    rownames_to_column(var = "cell_id") %>%
    as_tibble()

saveRDS(cell_type_assignments, file = "misc/cell_type_assignments")

In the final step, assignments were saved by cell ID to avoid erors related to future changes in clustering.

```{r assignCellTypes}
cell_type_assignments <- readRDS(file = "misc/cell_type_assignments")

cell_type_assignments <- left_join(
    tibble(cell_id = names(combined_data$cell_type_final)),
    cell_type_assignments)

combined_data$cell_type_final <- cell_type_assignments$assignment

cell_type_markers <- c(
    "Sftpc", "Cldn18", "Ager", # Epithelial
    "Pdpn", "Upk3b", "Cdh2", # Mesothelial
    "Nkg7", "Klrd1", "Gzma", # NK
    "Zbtb46", "Il12b", "Serpinb9", # Dendritic
    "Gsn", "Mfap4", "Lum", # Fibroblast
    "Plac8", "S100a4", "Ly6c2", # Monocyte
    "Cstb", "Spp1", "Fn1", # Interstitial Mac
    "Ms4a4c", "Irf7", # TIM
    "Cd36", "Pecam1", "Cd93", # Endothelial
    "Ccl8", "C1qa", "C1qb", "TAM", # TAM
    "Tcf7l2", "Axl", "Siglecf" # Alveolar Mac
)

DotPlot(combined_data, features = cell_type_markers,
    cols = c("lightgoldenrod", "darkred")) +
    theme(axis.text.x = element_text(
        angle = 90, vjust = 0))

qs::qsave(combined_data,
          file = "output/rdata/combined_data.qs")
```

### Run DE on the mouse samples
```{r deGrowthFactorMetsDE, dependson='deGrowthFactorMets'}
combined_data <- qs::qread("output/rdata/combined_data.qs")
Idents(combined_data) <- combined_data$group_name

overall_de <-
    FindMarkers(combined_data,
                ident.1 = "Mouse control lung",
                logfc.threshold = 0)
qs::qsave(overall_de,
          file = "output/rdata/overall_de.qs")
```

### Count cell types
```{r deGrowthFactorMetsCellType, dependson='deGrowthFactorMetsDE'}
cell_type_list <-
    table(combined_data$cell_type_final, combined_data$group) %>%
    as.data.frame() %>%
    group_by(Var1) %>%
    summarize(sum = sum(Freq), .groups = "drop") %>%
    filter(sum > 10)
```

### Make supplemental figures with all cell types
```{r deGrowthFactorMetsLoop, dependson='deGrowthFactorMetsCellType'}
overall_de <-
    qs::qread("output/rdata/overall_de.qs")

target_genes <-
    read_tsv("misc/growth_factor_genes.txt",
             col_names = "human_gene",
             show_col_types = FALSE) %>%
    left_join(read_tsv("misc/HMD_HumanPhenotype.rpt",
                       col_names = c("human_gene",
                                     "h_entrez",
                                     "mouse_gene",
                                     "mgi",
                                     "pheno",
                                     "junk"),
                       show_col_types = FALSE) %>%
        select(human_gene, mouse_gene),
        by = "human_gene") %>%
    select(mouse_gene)

target_genes <-
    left_join(target_genes,
              overall_de %>%
                rownames_to_column("mouse_gene"))

################################################################################
# DE between lung and primary mets by cell type
combined_subset <-
    combined_data %>%
    subset(cell_type_final %in% cell_type_list$Var1)

combined_subset <-
    combined_subset[rownames(combined_subset) %in% target_genes$mouse_gene, ]

combined_subset$four_groups <-
    combined_subset$group %>%
    str_replace_all("_", " ") %>%
    factor(levels = c("Mouse control lung",
                      "Early mouse stroma",
                      "Midpoint mouse stroma",
                      "Late mouse stroma"))

Idents(combined_subset) <- combined_subset$group_name

qs::qsave(combined_subset,
          file = "output/rdata/combined_subset.qs")

table(combined_subset$cell_type_final, combined_subset$group_name) %>%
    as.data.frame() %>%
    pivot_wider(names_from = Var2,
                values_from = Freq)# %>%
    # gt::gt() %>%
    # gt::gtsave(filename = "output/figures/cell_type_table.pdf")
    ########### Clean this up

cell_type_de <- tibble(gene = character())

for (cell_type in unique(combined_subset$cell_type_final)) {
    cell_type_de <-
        subset(combined_subset,
               subset = cell_type_final == {{ cell_type }}) %>%
        FindMarkers(ident.1 = "Mouse control lung",
                    logfc.threshold = 0,
                    min.pct = 0) %>%
        as.data.frame() %>%
        rownames_to_column("gene") %>%
        mutate(cell_type = cell_type) %>%
        select(gene, avg_log2FC, p_val_adj, cell_type) %>%
        arrange(p_val_adj) %>%
        bind_rows(cell_type_de)
}

write_tsv(cell_type_de,
          file = "output/de/cell_type_de.txt")

for (factor_name in c("group_name", "four_groups")) {
    ############################################################################
    # Plot DE between lung and primary mets by cell type
    temp_data <-
        GetAssayData(combined_subset) %>%
        as.matrix() %>%
        t() %>%
        as.data.frame() %>%
        select(target_genes$mouse_gene[target_genes$mouse_gene %in%
                                          rownames(combined_subset)]) %>%
        rownames_to_column("cell") %>%
        mutate(cell_type = combined_subset$cell_type_final,
               group_col = combined_subset[[factor_name]] %>% pull()) %>%
        pivot_longer(cols = c(-group_col,
                             -cell_type,
                             -cell),
                     names_to = "gene",
                     values_to = "expression") %>%
        group_by(cell_type, group_col, gene) %>%
        summarize(expr_mean = mean(expression),
                  num_cells = n(),
                  pct_expr = sum(expression > 0) / n(),
                  .groups = "drop") %>%
        mutate(gene_group = paste(gene, group_col, sep = "\n"))

    plot_name <-
        ggplot(temp_data,
           aes(x = group_col,
               y = cell_type,
               group = group_col,
               color = expr_mean,
               size = pct_expr)) +
        geom_point() +
        facet_wrap(~ gene, nrow = 1) +
        theme_robert() +
        scale_color_gradient(low = "white",
                             high = "darkred") +
        theme(axis.text.x = element_text(angle = 90,
                                         hjust = 1,
                                         vjust = 0.5),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank()) +
        scale_y_discrete(limits = rev) +
        labs(x = "",
             y = "",
             title = paste0("Gene Expression Split by ",
                            str_replace(factor_name, "_", " ") %>%
                                str_to_title(),
                            " and Cell Type"))

    ggsave(paste0("output/figures/dotplot_cell_type_by_",
                  factor_name,
                  ".pdf"),
           plot = plot_name,
           width = 15,
           height = 8)

    plot_name <-
        temp_data %>%
        left_join(cell_type_de) %>%
        mutate(fct_order = as.numeric(group_col),
               group = paste0(group_col,
                              " (n=",
                              num_cells,
                              ")") %>%
                   factor() %>%
                   fct_reorder(fct_order),
               p_val_adj = replace_na(p_val_adj, 1),
               sig = ifelse(p_val_adj <= 0.05,
                            "Significant",
                            "Non-significant")) %>%
        ggplot(aes(x = group,
                   y = gene,
                   fill = pct_expr)) +
        geom_tile(aes(color = sig),
                  size = 0.5,
                  width = 0.9,
                  height = 0.9) +
        geom_text(aes(label = sprintf("%.1f%%", pct_expr * 100)),
                  size = 1,
                  color = "lightgray") +
        facet_wrap(~ cell_type, nrow = 1, scales = "free_x") +
        labs(x = "",
             y = "",
             title = paste0("Percent of Cells Expressing Select Gene")) +
        theme_robert() +
        scale_fill_continuous(labels = scales::label_percent(accuracy = 2),
                              name = "Percent of\nCells Expressing Gene") +
        scale_color_manual(values = c(plot_cols[1], alpha("black", 0)),
                           name = NULL,
                           breaks = c("Significant", "Non-significant")) +
        theme(legend.position = "right",
              axis.text.x = element_text(angle = 90, hjust = 1)) +
        guides(color = guide_legend(override.aes = list(size = 2, alpha = 1)))

    ggsave(file = paste0("output/figures/growth_factor_perc_expressing_",
                         factor_name,
                         ".pdf"),
        plot = plot_name,
        width = 10,
        height = 8)
}
```

## Make main figures for growth factors

############################################ Maybe use scaled gene expression values?
I'd like the following cell types to be in the main figure: dendritic cells, fibroblasts, macrophages, and monocytes And the following genes: AREG, EREG, FGF2, FGF10, FGF18, HBEGF, HGF, IFG1, NRG1, PDGFa, PDGFb, PDGFc, TGFa, TGFb1, and VEGFa
```{r main_figs, dependson='deGrowthFactorMets'}
combined_subset <- qs::qread("output/rdata/combined_subset.qs")

cell_type_de <- read_tsv("output/de/cell_type_de.txt",
                         show_col_types = FALSE)

main_figure_genes <-
    c("Ereg", "Fgf2", "Hbegf", "Nrg1", "Vegfa")

main_figure_cell_types <-
    c("TAM", "TIM", "Fibroblast", "Alveolar Mac", "Monocyte", "Mesothelial")

combined_subset_main <-
    subset(combined_subset,
           subset = cell_type_final %in% main_figure_cell_types)

combined_subset_main <-
    combined_subset_main[main_figure_genes, ]

for (factor_name in c("group_name", "four_groups")) {
    ############################################################################
    # Plot DE between lung and primary mets by cell type
    temp_data <-
        GetAssayData(combined_subset_main) %>%
        as.matrix() %>%
        t() %>%
        as.data.frame() %>%
        rownames_to_column("cell") %>%
        mutate(cell_type = combined_subset_main$cell_type_final,
               group_col = combined_subset_main[[factor_name]] %>% pull()) %>%
        pivot_longer(cols = c(-group_col,
                              -cell_type,
                              -cell),
                     names_to = "gene",
                     values_to = "expression") %>%
        filter(cell_type %in% main_figure_cell_types) %>%
        group_by(cell_type, group_col, gene) %>%
        summarize(expr_mean = mean(expression),
                  pct_expr = sum(expression > 0) / n(),
                  num_cells = n(),
                  .groups = "drop") %>%
        mutate(gene_group = paste(gene, group_col, sep = "\n")) %>%
        group_by(gene) %>%
        mutate(norm_expr = scale(expr_mean, center = FALSE))

    plot_name <-
        ggplot(temp_data,
           aes(x = group_col,
               y = cell_type,
               color = norm_expr,
               size = pct_expr)) +
        geom_point() +
        facet_wrap(~ gene, nrow = 1) +
        theme_robert() +
        scale_color_gradient(low = "white",
                             high = "darkred",
                             name = "Normalized\nexpression") +
        scale_size_continuous(name = "% of cells\nexpressing gene") +
        theme(axis.text.x = element_text(angle = 90,
                                         hjust = 1,
                                         vjust = 0.5),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank()) +
        scale_y_discrete(limits = rev) +
        labs(x = "",
             y = "Expression",
             title = paste0("Gene Expression Split by ",
                            str_replace(factor_name, "_", " ") %>%
                                str_to_title(),
                            " and Cell Type"))

    ggsave(paste0("output/figures/dotplot_cell_type_by_",
                  factor_name,
                  "_main.pdf"),
           plot = plot_name,
           width = 8,
           height = 4)

    plot_name <-
        temp_data %>%
        left_join(cell_type_de) %>%
        mutate(fct_order = as.numeric(group_col),
               group = paste0(group_col,
                              " (n=",
                              num_cells,
                              ")") %>%
                   factor() %>%
                   fct_reorder(fct_order),
               p_val_adj = replace_na(p_val_adj, 1),
               sig = ifelse(p_val_adj <= 0.05,
                            "Significant",
                            "Non-significant")) %>%
        ggplot(aes(x = group,
                   y = gene,
                   fill = pct_expr)) +
        geom_tile(aes(color = sig),
                  size = 1,
                  width = 0.9,
                  height = 0.9) +
        geom_text(aes(label = sprintf("%.1f%%", pct_expr * 100)),
                  size = 1.5,
                  color = "lightgray") +
        facet_wrap(~ cell_type, nrow = 1, scales = "free_x") +
        labs(x = "",
             y = "",
             title = paste0("Percent of Cells Expressing Select Gene")) +
        theme_robert() +
        scale_fill_continuous(labels = scales::label_percent(accuracy = 2),
                              name = "Percent of\nCells Expressing Gene") +
        scale_color_manual(values = c(plot_cols[1], alpha("black", 0)),
                           name = NULL,
                           breaks = c("Significant", "Non-significant")) +
        theme(legend.position = "right",
              axis.text.x = element_text(angle = 90, hjust = 1)) +
        guides(color = guide_legend(override.aes = list(size = 2, alpha = 1)))

    ggsave(file = paste0("output/figures/growth_factor_perc_expressing_",
                         factor_name,
                         "_main.pdf"),
        plot = plot_name,
        width = 8,
        height = 6)
}
```

## Human Down Sampling for human
```{r human_downsample}
#load the sobj
sobj_list <- qs::qread("output/rdata/sobj_list.qs")

#subset human only
human_objects_list <-
    sobj_list[c("tibia_tumor_h",
                "d14_pos_h",
                "d28_pos_h",
                "d49_pos_h")]

set.seed(3333)
#Downsample
for (i in seq_len(length(human_objects_list))) {
    temp_obj <- human_objects_list[[i]]
    current_n_cells <- length(Cells(temp_obj))
    sub_n_cells <- min(1500, current_n_cells)
    temp_obj <- temp_obj[, sample(1:current_n_cells, sub_n_cells)]
    human_objects_list[[i]] <- temp_obj
}

#Dimplot with down sample
downsample_umap_human_only <-
    merge(human_objects_list[["tibia_tumor_h"]],
            y = human_objects_list[c("d14_pos_h",
                                     "d28_pos_h",
                                     "d49_pos_h")],
            add.cell.ids = c("Primary_tumor",
                             "Early_metastatic_tumor",
                             "Midpoint_metastatic_tumor",
                             "Late_metastatic_tumor"),
            project = "downsample_umap_human_only") %>%
            process_seurat()

downsample_umap_human_only$group <-
    factor(str_replace_all(downsample_umap_human_only$group, "_", " "),
           levels = c("Primary tumor",
                      "Early metastatic tumor",
                      "Midpoint metastatic tumor",
                      "Late metastatic tumor"))

#Dimplot
downsample_human_combined <-
    DimPlot(downsample_umap_human_only) +
    ggtitle(label = "Primary and metastatic tumor") +
    theme_robert() +
    theme(plot.title = element_text(hjust = 0.5)) +
    coord_fixed()

ggsave("output/figures/downsample_human_combined.pdf",
       width = 6,
       height = 6,
       plot = downsample_human_combined)

#split by
downsample_human_split <-
    DimPlot(downsample_umap_human_only, split.by = "group") +
            theme_robert() +
            theme(plot.title = element_text(hjust = 0.5)) +
            NoLegend() +
            coord_fixed()

ggsave("output/figures/downsample_human_split.pdf",
       width = 12,
       height = 4,
       plot = downsample_human_split)
```

## Mouse Down Sampling and Dimplot
```{r mouse_downsample}
#subset human only
sobj_list <- qs::qread("output/rdata/sobj_list.qs")
combined_data <- qs::qread("output/rdata/combined_data.qs")

Idents(combined_data) <- combined_data$cell_type_final

#subset the combined mouse object so that you can downsample
set.seed(4444)
mouse_annotated_list <- list()
for (stuff in unique(combined_data$obj_name)) {
    print(stuff)
    temp_obj <- subset(combined_data,
                       subset = obj_name == stuff)
    mouse_annotated_list[[stuff]] <- temp_obj
}

#Downsample
for (i in seq_len(length(mouse_annotated_list))) {
    temp_obj <- mouse_annotated_list[[i]]
    current_n_cells <- length(Cells(temp_obj))
    sub_n_cells <- min(4500, current_n_cells)
    temp_obj <- temp_obj[, sample(1:current_n_cells, sub_n_cells)]
    mouse_annotated_list[[i]] <- temp_obj
}

#downsampled objects are merged and seurat process for dimplot
downsample_umap_mouse_only <-
    merge(mouse_annotated_list[["S0006_mouse"]],
            y = mouse_annotated_list[c("d14_pos_m",
                                       "d28_pos_m",
                                       "d49_pos_m")],
            add.cell.ids = c("Mouse_control_lung",
                             "Early_mouse_stroma",
                             "Midpoint_mouse_stroma",
                             "Late_mouse_stroma"),
            project = "downsample_umap_mouse_only") %>%
            process_seurat()

downsample_umap_mouse_only$group <-
    factor(str_replace_all(downsample_umap_mouse_only$group,
                           "_",
                           " "),
           levels = c("Mouse control lung",
                      "Early mouse stroma",
                      "Midpoint mouse stroma",
                      "Late mouse stroma"))

#combined
downsample_mouse_combined <-
    DimPlot(downsample_umap_mouse_only, group.by = "cell_type_final") +
            theme_robert() +
            labs(title = "Mouse control and stroma") +
            theme(plot.title = element_text(hjust = 0.5)) +
            coord_fixed()

#split by
downsample_mouse_split <-
            DimPlot(downsample_umap_mouse_only, split.by = "four_groups") +
            theme_robert() +
            theme(plot.title = element_text(hjust = 0.5)) +
            NoLegend() +
            coord_fixed()

#save the plot
ggsave("output/figures/downsample_mouse_combined.pdf",
       width = 6,
       height = 6,
       plot = downsample_mouse_combined)

ggsave("output/figures/downsample_mouse_split.pdf",
       width = 12,
       height = 4,
       plot = downsample_mouse_split)

#all combined
plots_all <-
    downsample_human_combined + downsample_human_split +
    downsample_mouse_combined + downsample_mouse_split +
    plot_layout(widths = c(1, 4))

#save the plot
ggsave("output/figures/plots_all.pdf",
       width = 10,
       height = 6,
       plot = plots_all)
```

## Non-downsample UMAP plots
```{r}
#load the sobj
sobj_list <- qs::qread("output/rdata/sobj_list.qs")

#subset the sobj of interest for human and mouse and merge
#for human
umap_human_only <-
    merge(sobj_list[["tibia_tumor_h"]],
          y = sobj_list[c("d14_pos_h",
                            "d28_pos_h",
                            "d49_pos_h")],
          add.cell.ids = c("Primary_tumor",
                            "Early_metastatic_tumor",
                            "Midpoint_metastatic_tumor",
                            "Late_metastatic_tumor"),
           project = "umap_human_only") %>%
    process_seurat()

#Human only combined umap: this puts them in the same space/comparable
combined_human_only <-
    DimPlot(umap_human_only) +
    ggtitle(label = "Primary and metastatic tumor") +
    theme_robert() +
    theme(plot.title = element_text(hjust = 0.5)) +
    coord_fixed()

combined_human_only

#Dimplot
human_umaps <- DimPlot(umap_human_only, split.by = "group") +
            theme_robert() +
            theme(plot.title = element_text(hjust = 0.5)) +
            NoLegend() + coord_fixed()

#rearrange the panels for umap space
human_umaps$data$group <-
    factor(x = human_umaps$data$group,
           levels = c("Primary_tumor",
                      "Early_metastatic_tumor",
                      "Midpoint_metastatic_tumor",
                      "Late_metastatic_tumor"))

#for mouse
umap_mouse_only <-
    merge(sobj_list[["S0006_mouse"]],
            y = sobj_list[c("d14_pos_m",
                            "d28_pos_m",
                            "d49_pos_m")],
            add.cell.ids = c("Mouse_control_lung",
                            "Early_mouse_stroma",
                            "Midpoint_mouse_stroma",
                            "Late_mouse_stroma"),
            project = "umap_mouse_only") %>%
            process_seurat()

#Dimplot
mouse_umaps_combined <- DimPlot(umap_mouse_only) +
                        theme_robert() +
                        theme(plot.title = element_text(hjust = 0.5)) +
                        NoLegend() +
                        coord_fixed()

mouse_umaps <- DimPlot(umap_mouse_only, split.by = "group") +
               theme_robert() +
               theme(plot.title = element_text(hjust = 0.5)) +
               NoLegend() +
               coord_fixed()

#rearrange the panels for umap space
mouse_umaps$data$group <- factor(x = mouse_umaps$data$group,
                                levels = c("Mouse_control_lung",
                                        "Early_mouse_stroma",
                                        "Midpoint_mouse_stroma",
                                        "Late_mouse_stroma"))

combined_plots <- (human_umaps / mouse_umaps)

combined_plots_all <- (combined_human_only + (human_umaps / mouse_umaps))
                         + plot_layout(widths = c(1, 4))

#save the plot
ggsave("output/figures/combined_human_only.pdf",
       width = 6,
       height = 6,
       plot = combined_plots)

ggsave("output/figures/combined_plots.pdf",
       width = 8,
       height = 6,
       plot = combined_plots)

ggsave("output/figures/combined_plots_all.pdf",
       width = 12,
       height = 6,
       plot = combined_plots)


```


## GSEA_KEGG
GSEA-KEGG figure showing upregulated genes/pathways as lesions progress from:
                  S0037 (primary) redo human cells to S0082(early) human cells
                  S0082 (early) human cell to S0084(midpoint) human cells
                  S0084(midpoint) human cells to S0093 (late) human cells
                  S0082(early) human cells to S0093(late) human cells
### Run DE for KEGG
```{r}
sobj_list <- qs::qread("output/rdata/sobj_list.qs")
#we want all of the item in the same space, do combine and run process seurat
human_only <-
    merge(sobj_list[["tibia_tumor_h"]],
          y = sobj_list[c("d14_pos_h",
                          "d28_pos_h",
                          "d49_pos_h")],
          add.cell.ids = c("Primary_tumor",
                           "Early_metastatic_tumor",
                           "Midpoint_metastatic_tumor",
                           "Late_metastatic_tumor"),
          project      =   "human_only") %>%
          process_seurat()

#set the idents
Idents(human_only) <- "group"

#make a table to do the find markers in loop
marker_table <-
    tribble(~item1,                      ~item2,                        ~cmp,
            "Primary_tumor",             "Early_metastatic_tumor",      "Primary_tumor_vs_Early_metastatic_tumor",
            "Early_metastatic_tumor",    "Midpoint_metastatic_tumor",   "Early_metastatic_tumor_vs_Midpoint_metastatic_tumor",
            "Early_metastatic_tumor",    "Late_metastatic_tumor",       "Early_metastatic_tumor_vs_Late_metastatic_tumor",
            "Midpoint_metastatic_tumor", "Late_metastatic_tumor",       "Midpoint_metastatic_tumor_vs_Late_metastatic_tumor",
            "Late_metastatic_tumor",     "Primary_tumor",               "Late_metastatic_tumor_vs_Primary_tumor",
            "Early_metastatic_tumor",    c("Midpoint_metastatic_tumor",
                                           "Late_metastatic_tumor"),    "Early_vs_Mid_Late",
            "Early_metastatic_tumor",    c("Midpoint_metastatic_tumor",
                                           "Late_metastatic_tumor",
                                           "Primary_tumor"),            "Early_vs_Mid_Late_Primary")


#make a list that you can append to
#FindMarkers since we comaring across only two groups
gsea_lists <-
    parallel::mclapply(1:nrow(marker_table), function(i) {
        df <- FindMarkers(human_only,
                          ident.1 = marker_table$item1[[i]],
                          ident.2 = marker_table$item2[[i]],
                          min.pct = 0.1,
                          logfc.threshold = 0) %>%
            rownames_to_column("gene")
        filename1 <-
            paste0("output/de/",
                   marker_table$cmp[[i]],
                   ".tsv")
        write_tsv(as.data.frame(df), file = filename1)
        return(df)
    }, mc.cores = parallelly::availableCores())

names(gsea_lists) <- marker_table$cmp

# Drop Early_vs_Mid_Late since we're just using this for IPA
gsea_lists <-
    gsea_lists[!names(gsea_lists) %in% c("Early_vs_Mid_Late",
                                         "Early_vs_Mid_Late_Primary")]

qs::qsave(gsea_lists, file = "output/rdata/gsea_gene_lists.qs")
```


## Ran IPA on the output of FindMarkers above
Used a logFC cutoff of -1/1 and p_val_adj cutoff of 0.05
```{r plotIPA}
path_top_5 <-
    read_tsv(list.files("output/pathway_analysis",
                        pattern = "_paths.txt",
                        full.names = TRUE),
         id = "sample",
         col_names = c("pathway",
                       "neg_log_pval",
                       "ratio",
                       "z_score",
                       "genes"),
         show_col_types = FALSE,
         skip = 1) %>%
    filter(neg_log_pval > log10(0.05) &
           !is.nan(z_score) &
           z_score != 0) %>%
    arrange(sample, desc(abs(z_score))) %>%
    mutate(sample = str_remove(sample, ".+/") %>%
               str_remove("_paths.txt")) %>%
    group_by(sample, z_score > 0) %>%
    slice_head(n = 5) %>%
    group_by(sample) %>%
    arrange(desc(z_score), .by_group = TRUE) %>%
    select(sample, pathway, z_score) %>%
    mutate(order = seq_len(n()),
           justify_y = if_else(z_score > 0, 1, 0),
           y_pos = if_else(z_score > 0, -0.1, 0.1),
           sample = str_remove(sample, "_tumor$") %>%
                str_replace_all("metastatic", "metastasis") %>%
                str_replace("tumor", "vs") %>%
                str_replace_all("_", " ") %>%
                str_to_sentence() %>%
                factor(., levels = c("Primary vs early metastasis",
                                     "Early metastasis vs midpoint metastasis",
                                     "Midpoint metastasis vs late metastasis",
                                     "Late metastasis vs primary"))) %>%
    ungroup()

two_way_plot <- function(data,
                         x_col = "z_score",
                         ylim = c(-5, 5),
                         lab_y = 3) {
    lab4plot <-
        tibble(y = c(-lab_y, lab_y),
            x = c(0.2, 0.2),
            label = as.factor(c(paste("Upregulated in",
                stringr::str_remove(data$sample[1], ".+vs") %>%
                stringr::str_to_lower()),
            paste("Upregulated in",
                str_remove(data$sample[1], "vs.+") %>%
                str_to_lower()))))

    plot_name <-
        ggplot() +
        geom_bar(data = data,
                 aes(x = -1 * order,
                     y = get(x_col),
                     fill = get(x_col) > 0),
                 stat = "identity",
                 alpha = 0.8) +
            coord_flip() +
            geom_hline(yintercept = 0,
                       color = "black",
                       linewidth = 0.5) +
            geom_text(data = data,
                      aes(x = -1 * order,
                          y = y_pos,
                          hjust = justify_y,
                          label = pathway),
                      size = 1.5) +
            geom_text(data = lab4plot,
                      aes(x = x,
                          y = y,
                          label = label),
                      fontface = "bold",
                      size = 1.5) +
            scale_fill_manual(values = plot_cols,
                              name = paste0(x_col, " > 0")) +
            theme_robert() +
            theme(panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.border = element_blank(),
                  axis.ticks.y = element_blank(),
                  axis.text.y = element_blank(),
                  axis.line.x = element_line(color = "black")) +
            labs(title = data$sample[1],
                 y = "",
                 x = "") +
            ylim(ylim)

    return(plot_name)
}

plot_list <- lapply(levels(path_top_5$sample),
                    function(x) {
    two_way_plot(filter(path_top_5, sample == x))
})

pdf("output/figures/Top_5_IPA_pathways.pdf",
    width = 4,
    height = 6)
# Adjusted heights to get the bars to be the "same" height
wrap_plots(plot_list, ncol = 1, heights = c(5, 5, 0.75, 1.22))
dev.off()
```

```{r plotIPA2}
dis_functs <-
    read_tsv("output/pathway_analysis/Early_vs_Mid_Late_dis_functs.txt",
         col_names = c("category",
                       "dis_funct",
                       "pval",
                       "predict",
                       "act_score",
                       "genes",
                       "num_genes"),
         show_col_types = FALSE,
         skip = 1) %>%
    mutate(sample = "Early vs Established Metastasis") %>%
    select(-category, -num_genes) %>%
    filter(pval <= 0.05 &
           !is.nan(act_score) &
           abs(act_score) >= 0 &
           !is.na(predict)) %>%
    arrange(sample, pval) %>%
    mutate(sample = str_remove(sample, ".+/") %>%
               str_remove("_dis_functs.txt")) %>%
    group_by(sample, act_score > 0) %>%
    slice_head(n = 5) %>%
    group_by(sample) %>%
    arrange(desc(act_score), .by_group = TRUE) %>%
    mutate(order = seq_len(n()),
           justify_y = if_else(act_score > 0, 1, 0),
           y_pos = if_else(act_score > 0, -0.1, 0.1)) %>%
    ungroup() %>%
    rename(pathway = dis_funct)

pdf("output/figures/IPA_dis_funct_early_vs_midlate.pdf",
    width = 4,
    height = 4)
two_way_plot(dis_functs,
             x_col = "act_score",
             ylim = c(-6, 6),
             lab_y = 3) +
    theme(legend.position = "none")
dev.off()
```

### Running the GSEA_KEGG for each pathway category
select dataset and then aggregate the genes into a list
```{r}
gsea_lists <- qs::qread("output/rdata/gsea_gene_lists.qs")
cat_tib <-
    tribble(~category,          ~subcat,
            "C2",               "CP:KEGG",
            "C4",               NULL,
            "C6",               NULL)

for (i in seq_len(nrow(cat_tib))) {
    message(paste("Starting analysis for", cat_tib$category[i]))
    kegg_ref <-
        msigdbr::msigdbr(species = "Homo sapiens",
                         category = cat_tib$category[i]) %>%
                         split(x = .$gene_symbol, f = .$gs_name)

    #Run the GSEA in loop for all the things
    kegg_results_list <-
        parallel::mclapply(names(gsea_lists),
                           function(item) {

            gsea_input <-
                gsea_lists[[item]] %>%
                rownames_to_column("gene") %>%
                arrange(desc(avg_log2FC)) %>%
                pull(avg_log2FC, name = gene)

            result <-
                fgsea::fgseaMultilevel(kegg_ref,
                                       gsea_input,
                                       minSize = 15,
                                       maxSize = 500,
                                       nPerm = 1000) %>%
                arrange(desc(NES)) %>%
                    mutate(pathway =
                        str_replace_all(pathway, "_", " ") %>%
                        str_remove("KEGG_")) %>%
                    filter(padj < 0.05)

            write_tsv(result,
                      file = paste0("output/gsea/",
                                    item,
                                    "category_",
                                    cat_tib$category[i],
                                    ".tsv"))
            return(result)
    }, mc.cores = parallelly::availableCores())

    names(kegg_results_list) <- names(gsea_lists)

    for (item in names(gsea_lists)) {
        top5up_down <-
            rbind(kegg_results_list[[item]] %>%
                    filter(NES > 0) %>%
                    slice_head(n = 5),
                kegg_results_list[[item]] %>%
                    filter(NES < 0) %>%
                    slice_tail(n = 5)) %>%
                dplyr::select(pathway,
                              NES,
                              size,
                              padj) %>%
                arrange(desc(NES)) %>%
                mutate(pathway = as.factor(pathway) %>%
                            str_wrap(80) %>%
                            fct_reorder(NES),
                       sample = item,
                       order = 1:n(),
                       justify_y = if_else(NES > 0, 1, 0),
                       y_pos = if_else(NES > 0, -0.1, 0.1),
                       sample = str_remove(sample, "_tumor$") %>%
                            str_replace_all("metastatic", "metastasis") %>%
                            str_replace_all("_", " ") %>%
                            str_replace("tumor vs", "vs") %>%
                            str_to_sentence() %>%
                            factor(.,
                                   levels = c("Primary vs early metastasis",
                                              "Early metastasis vs midpoint metastasis",
                                              "Midpoint metastasis vs late metastasis",
                                              "Late metastasis vs primary")))

        pdf(paste0("output/figures/gsea_",
                   cat_tib$category[i],
                   "_",
                   item,
                   "barplots.pdf"),
            width = 4,
            height = 3)
        print(two_way_plot(top5up_down, x_col = "NES"))
        dev.off()
    }
}
```

## Ran IPA on the output of FindMarkers above
Used a logFC cutoff of -1/1 and p_val_adj cutoff of 0.05
```{r plotIPA}
path_top_5 <-
    read_tsv(list.files("output/pathway_analysis",
                        pattern = "_paths.txt",
                        full.names = TRUE),
         id = "sample",
         col_names = c("pathway",
                       "neg_log_pval",
                       "ratio",
                       "z_score",
                       "genes"),
         show_col_types = FALSE,
         skip = 1) %>%
    filter(neg_log_pval > log10(0.05) &
           !is.nan(z_score) &
           z_score != 0) %>%
    mutate(sample = str_remove(sample, ".+/") %>%
               str_remove("_paths.txt")) %>%
    group_by(sample, z_score > 0) %>%
    arrange(desc(abs(z_score)), .by_group = TRUE) %>%
    slice_head(n = 5) %>%
    select(sample, pathway, z_score) %>%
    mutate(order = seq_len(n()),
           label_y = z_score / abs(z_score),
           sample = str_remove(sample, "_tumor$") %>%
                str_replace_all("metastatic", "metastasis") %>%
                str_replace("tumor", "vs") %>%
                str_replace_all("_", " ") %>%
                str_to_sentence() %>%
                factor(., levels = c("Primary vs early metastasis",
                                     "Early metastasis vs midpoint metastasis",
                                     "Midpoint metastasis vs late metastasis",
                                     "Late metastasis vs primary"))) %>%
    ungroup()

two_way_plot <- function(data) {
    lab4plot <-
        tibble(y = c(-3, 3),
            x = c(0.2, 0.2),
            label = as.factor(c(paste("Upregulated in",
                stringr::str_remove(data$sample[1], ".+vs") %>%
                stringr::str_to_lower()),
            paste("Upregulated in",
                str_remove(data$sample[1], "vs.+") %>%
                str_to_lower()))))

    plot_name <-
        ggplot() +
        geom_bar(data = data,
                 aes(x = -1 * order,
                     y = z_score,
                     fill = z_score > 0),
                 stat = "identity",
                 alpha = 0.8) +
            coord_flip() +
            geom_hline(yintercept = 0,
                       color = "black",
                       linewidth = 0.5) +
            geom_text(data = data,
                      aes(x = -1 * order,
                          y = label_y * 3,
                          label = pathway),
                      size = 1.5) +
            geom_text(data = lab4plot,
                      aes(x = x,
                          y = y,
                          label = label),
                      fontface = "bold",
                      size = 1.5) +
            scale_fill_manual(values = plot_cols) +
            # theme_robert(title_font_size = 5,
            #               axis_font_size = 4) +
            theme(strip.background = element_rect(color = "white",
                                                  fill = "white"),
                  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.line.y = element_blank(),
                  axis.text.y = element_blank(),
                  axis.ticks.y = element_blank(),
                  plot.title = element_text(size = 5, face = "bold"),
                  axis.text.x = element_text(size = 3)) +
            labs(title = data$sample[1],
                 y = "",
                 x = "") +
            ylim(-5, 5)

    return(plot_name)
}

plot_list <- lapply(levels(path_top_5$sample),
                    function(x) {
    two_way_plot(filter(path_top_5, sample == x))
})

pdf("output/figures/Top_5_IPA_pathways.pdf",
    width = 4,
    height = 4)
wrap_plots(plot_list, ncol = 1, heights = c(5, 5, 1.8, 3))
dev.off()
```

# Nichenetr
sender : all stroma/mouse cells from the early+mid
receiver : Human cells/tumor from the early+mid
Genes/markers tested : differential expression from ident1=early and ident2=mid
```{r}
combined_data <- qs::qread("output/rdata/combined_data.qs")

# Seurat object: sender of signal
sender_sobj <- combined_data %>%
    subset(cell_type_final != "Epithelial" &
        cell_type_final != "Endothelial")
Idents(sender_sobj) <- "cell_type_final"
sender_sobj$nichenet_id <- sender_sobj$cell_type_final

# Seurat object: receiver of signal
receiver_sobj <- human_only %>%
    subset(group != "Late_metastatic_tumor" &
        group != "Midpoint_metastatic_tumor")
receiver_sobj$group[receiver_sobj$group == "Early_metastatic_tumor"] <-
    "Early_metastasis"
receiver_sobj$nichenet_id <- receiver_sobj$group

# Seurat object: combined objects
sender_receiver_sobj <- merge(sender_sobj,
          y = receiver_sobj,
          add.cell.ids = c("sender", "receiver"),
          project      =   "sender_receiver_sobj") %>%
          process_seurat()
```

## Get the pre-defined model data
```{r}
ligand_target_matrix <-
    readRDS("nichenetrfiles/ligand_target_matrix.rds")

# Use a strict definition of the network
lr_network <-
    readRDS("nichenetrfiles/lr_network.rds") %>%
    filter(database != "ppi_prediction_go" & database != "ppi_prediction")

weighted_networks <-
    readRDS("nichenetrfiles/weighted_networks.rds")

weighted_networks_lr <-
    weighted_networks$lr_sig %>%
    inner_join(lr_network %>% distinct(from, to),
               by = c("from", "to"))
```


## Convert the gene names in model to mouse
```{r}
#makes capital lettered human genes (CXCL) to mouse format (Cxcl)
lr_network <-
    lr_network %>%
    mutate(from = convert_human_to_mouse_symbols(from),
           to = convert_human_to_mouse_symbols(to)) %>%
    drop_na()

#converts the colnames human genes name format to mouse format
colnames(ligand_target_matrix) <-
    ligand_target_matrix %>%
    colnames() %>%
    convert_human_to_mouse_symbols()

#converts the rownames human genes name format to mouse format
rownames(ligand_target_matrix) <-
    ligand_target_matrix %>%
    rownames() %>%
    convert_human_to_mouse_symbols()

#removes the NAs or the genes that were umable to convert
ligand_target_matrix <-
    ligand_target_matrix %>%
    .[!is.na(rownames(ligand_target_matrix)),
        !is.na(colnames(ligand_target_matrix))]

#Convert the human format genes to mouse format in the file and drops NAs
weighted_networks_lr <-
    weighted_networks_lr %>%
    mutate(from = convert_human_to_mouse_symbols(from),
           to = convert_human_to_mouse_symbols(to)) %>%
    drop_na()
```

## Nichenet Analysis

If it has not been run previously, run:
rrrSingleCellUtils:::gen_lig_receptor_ref()
to download the models to your cache.

```{r}
# Set the same idents for all the objects
Idents(sender_sobj) <- "nichenet_id"
Idents(receiver_sobj) <- "nichenet_id"
Idents(sender_receiver_sobj) <- "nichenet_id"

# Extract the genesets identified in the GSEA
pathways <- c("PDGF_UP.V1_DN", "PDGF_UP.V1_UP", "VEGF_A_UP.V1_DN",
    "VEGF_A_UP.V1_DN", "ERBB2_UP.V1_UP", "ERBB2_UP.V1_DN",
    "PGF_UP.V1_UP", "PGF_UP.V1_DN")

pathway_genes <- msigdbr(species = "human", category = "C6") %>%
    filter(gs_name %in% pathways) %>%
    pull(gene_symbol)

# Gives the list of target genes from the receiver dataset
gset_stuff <- find_tar_genes(receiver_sobj,
        id1 = c("Early_metastasis"),
        id2 = c("Primary_tumor"),
        pval = 0.05,
        logfc = 0.25,
        spec = "human") %>%
    intersect(pathway_genes)

# Finds the ligands and genes target relationship in the sender and receiver:
# Basically runs the nichenetr analysis
fl_object <- find_ligands(sender_receiver_sobj,
    gset = gset_stuff,
    receiver = unique(Idents(receiver_sobj)),
    senders = unique(Idents(sender_sobj)),
    gset_spec = "human",
    rec_spec = "human",
    send_spec = "mouse")

# Makes the plot, plots (fl_object[[7]]) object. Check (fl_object[[7]])
complex_heatmap <- plot_complex_heatmap(fl_object,
    grid_color_low = "gray95",
    grid_color_high = "steelblue4",
    point_color_low = "gray95",
    point_color_high = "darkred")

complex_heatmap
```

```{r, cache=FALSE}
sessionInfo()
```