HSVref_A673.Rmd

---
title: "Repeat A673 tumor only scRNAseq with HSV reference"
author: "Emily Franz"
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, cache = FALSE}
library(Seurat)
library(tidyverse)
library(ggsci)
library(rrrSingleCellUtils)
library(ggplot2)
library(patchwork)
source("scSeurat-addhsv.R")
# Set the random generator seed so that results are reproducible.
set.seed(888)

# Create files for saving if do not exist
if(!file.exists("Data")){
  dir.create("Data")
}
if(!file.exists("Plots")){
  dir.create("Plots")
}
if(!file.exists("Plots/A673_1")){
  dir.create("Plots/A673_1")
}
```

## Tumor Samples - Human
Load all the different A673 datasets, note that HSV genes included but filter out mouse

```{r load, fig.height=5, fig.width=5}
sample_names <- c('S0136_HSV_exon_edits',
                  'S0137_HSV_exon_edits',
                  'S0138_HSV_exon_edits',
                  'S0139_HSV_exon_edits')

sample_labels <- c("1 - Control",
                   "2 - oHSV",
                   "3 - Trabectedin",
                   "4 - oHSV+Trabectedin")
# correlate sample names and sample labels
names(sample_labels) <- sample_names

sample_titles <- c("Control",
                   "oHSV",
                   "Trabectedin",
                   "Combination")
names(sample_titles) <- sample_names

# set max ncounts (will come back to adjust after vlnplot made, start at 1000000)
nCount_RNA_max <- list(1000000,
                       1000000,
                       1000000,
                       1000000)
names(nCount_RNA_max) <- sample_names

# set min for ncounts (will come back to adjust after vlnplot made, start at 0)
nCount_RNA_min <- c(0, 0, 0, 0)
names(nCount_RNA_min) <- sample_names

# set % mt genes (will come back to adjust after vlnplot made, start at 50)
mt_max <- c(50, 50, 50, 50)
names(mt_max) <- sample_names

for (item in sample_names) {
  obj <- tenx_load_qc(paste("/gpfs0/home2/gdrobertslab/lab/Counts/",
                            item,
                            "/filtered_feature_bc_matrix/",
                            sep = ""),
                      human_hsv = TRUE)
  plot1 <- FeatureScatter(obj,
                          feature1 = "nCount_RNA",
                          feature2 = "percent.mt") +
  theme(legend.position="none")
  
  plot2 <- FeatureScatter(obj,
                          feature1 = "nCount_RNA",
                          feature2 = "nFeature_RNA")  +
  theme(legend.position="none")
  
  print(plot1 + plot2 +
    plot_annotation(title = sample_titles[[item]]))

  obj <- subset(obj,
                subset = 
                  nCount_RNA < nCount_RNA_max[[item]] &
                  nCount_RNA > nCount_RNA_min[[item]] &
                  percent.mt < mt_max[[item]])
  
  assign(sample_titles[[item]], obj) #makes each sample an individual Seurat obj without overwriting
}
```

Look at HSV transcripts vs mt transcripts and HSV transcripts vs nCounts.

oHSV-treated sample:

```{r scatter_hsv, fig.height=5, fig.width=10}
plot1_hsv <- FeatureScatter(oHSV,
                            feature1 = "nCount_RNA",
                            feature2 = "percent.hsv") +
  theme(legend.position="none")

plot2_hsv <- FeatureScatter(oHSV,
                            feature1 = "percent.mt",
                            feature2 = "percent.hsv") +
  theme(legend.position="none")

plot1_hsv +
  plot2_hsv +
  plot_annotation(title = 'oHSV-treated Tumor')
```

Combination-treated sample:

```{r scatter_combo, fig.height=5, fig.width=10}
plot1_hsv <- FeatureScatter(Combination,
                        feature1 = "nCount_RNA",
                        feature2 = "percent.hsv") +
  theme(legend.position="none")

plot2_hsv <- FeatureScatter(Combination,
                        feature1 = "percent.mt",
                        feature2 = "percent.hsv") +
  theme(legend.position="none")

plot1_hsv +
  plot2_hsv +
  plot_annotation(title = 'Combination-treated Tumor')
```

Using the previous violin plots, choose cutoffs to proceed with for processing.

```{r load_final}
# set max ncounts 
nCount_RNA_max <- list(45000,
                       25000,
                       20000,
                       20000)
names(nCount_RNA_max) <- sample_names

# set min for ncounts 
nCount_RNA_min <- c(3000,
                    3000,
                    3000,
                    3000)
names(nCount_RNA_min) <- sample_names

# set % mt genes (will come back to adjust after vlnplot made, start at 50)
# very high percent mt in tumor samples - caution
mt_max <- c(20,
            25,
            20,
            25)
names(mt_max) <- sample_names

# do not need to plot violin plot since completed this previously to determine the cutoffs used here.
for (item in sample_names) {
  obj <- tenx_load_qc(paste("/gpfs0/home2/gdrobertslab/lab/Counts/",
                            item,
                            "/filtered_feature_bc_matrix/",
                            sep = ""),
                      human_hsv = TRUE,
                      violin_plot = FALSE)
  
  obj <- subset(obj,
                  subset = 
                    nCount_RNA < nCount_RNA_max[[item]] &
                    nCount_RNA > nCount_RNA_min[[item]] &
                    percent.mt < mt_max[[item]])
  
  obj$treatment <- sample_labels[[item]]
  
  # Make each sample an individual Seurat obj without overwriting
  assign(sample_titles[[item]], obj)
}

# Create a list of seurat objects
#seurat_samples <- list(Control,
#                       oHSV,
#                       Trabectedin,
#                       Combination)
# Correlate sample names and sample labels
#names(seurat_samples) <- sample_names

#save(seurat_samples,
#     file = "raw_seurat_samples.RData")
```

# Normalize and Process

```{r process, fig.height=3, fig.width=3}
A673_hsv <- merge(Control,
                   y = c(oHSV,
                         Trabectedin,
                         Combination),
                   add.cell.ids = c("Control",
                                    "oHSV",
                                    "Trabectedin",
                                    "Combination"),
                   project = "Trabectedin plus HSV1716 in A673 - repeat")

#### Normalize, scale, process, cluster ####
A673_hsv <- NormalizeData(A673_hsv) %>%
  ScaleData() %>%
  FindVariableFeatures()

#### Variable Features ####
# Use Elbow Plots to determine the PCs to use

A673_hsv <- RunPCA(A673_hsv,
                     features = VariableFeatures(object = A673_hsv))
ElbowPlot(A673_hsv,
          ndims = 30)
```

### Clustering Resolution {.tabset}

#### Res0.1
```{r umap, fig.height=5, fig.width=10}
A673_hsv <- FindNeighbors(A673_hsv, dims = 1:13) %>%
FindClusters(resolution = 0.1) %>%
  RunUMAP(dims = 1:13)

DimPlot(A673_hsv,
        reduction = "umap",
        split.by = "treatment",
        pt.size = 1,
        label = T) +
  coord_fixed()
```

#### Res0.2
```{r, umap2, fig.height=5, fig.width=10}
A673_hsv <- FindClusters(A673_hsv,
                           resolution = 0.2) %>%
  RunUMAP(dims = 1:13)

DimPlot(A673_hsv,
        reduction = "umap",
        split.by = "treatment",
        pt.size = 1,
        label = T) +
  coord_fixed()
```

#### Res0.3
```{r, umap3, fig.height=5, fig.width=10}
A673_hsv <- FindClusters(A673_hsv,
                           resolution = 0.3) %>%
  RunUMAP(dims = 1:13)

DimPlot(A673_hsv,
        reduction = "umap",
        split.by = "treatment",
        pt.size = 1,
        label = T) +
  coord_fixed()
```

### DGEA by Cluster

Differentially Expressed Genes by Cluster (res0.3):
```{r deg}
A673_hsv_markers <- FindAllMarkers(A673_hsv,
                                     min.pct = 0.25,
                                     logfc.threshold = 0.25,
                                     only.pos = T)
```

```{r degt, dependson="deg"}

# Hard to work with this function/datatable so this step removes genes not in the top 10 DEG for all clusters 
# Print top 10 in each cluster
A673_hsv_markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC) -> t

# Get dataframe into proper format 
t <- as.data.frame(t)
t$cluster <- as.numeric(as.character(t$cluster))
t <- t[order(t$cluster),]
t <- t[,c(7, 1:6)]
t <- t[,c(7, 1:6)]
#could also use arrange() to sort by column

# Create interactive datatable showing top 10 DEG in each cluster 
DT::datatable(t, 
              rownames = FALSE,
              extensions = c('FixedColumns',
                             'Buttons'),
              options = list(
                pageLength = 10,
                scrollX = TRUE,
                scrollCollapse = TRUE,
                dom = 'Bfrtip',
                buttons = c('copy',
                            'csv',
                            'excel')
              ))
```

Downsample so there is the same number of cells in each treatment group.

```{r downsample, fig.height=5, fig.width=10}
# Subset the highest number that makes all treatment groups equal in cell number
# subset # of cells from each ident of the seurat object 
Idents(A673_hsv) <- A673_hsv$treatment
A673_hsv_even <- subset(A673_hsv,
                          downsample = table(A673_hsv$treatment) %>%
                            min())

# Check that worked (should be 348)
## extract meta data for group numbers
library(data.table)
A673_hsv_even.d <- A673_hsv_even@meta.data %>%
  as.data.table

# the resulting md object has one "row" per cell
A673_hsv_even.d[, .N, by = "treatment"]

DimPlot(A673_hsv_even,
        split.by = "treatment") +
  coord_fixed()

```

# oHSV Analysis

```{r matthelp_hsvgenes, dependson = c("umap","downsample")}
# Find head/prefix ^HSV1 in all genes (rownames) of tumor object
# Output is a vector of genes beginning with HSV1 that were found in the data:

hsv <- str_subset(row.names(A673_hsv_even),
                  pattern = "^HSV1")
```

**All oHSV Genes Present:**

```{r hsv_genes, dependson = c("umap","downsample","matthelp_hsvgenes")}
# Create interactive datatable showing hsv genes
hsv_dt <- as.data.table(hsv)
DT::datatable(hsv_dt, 
              rownames = FALSE,
              extensions = c('FixedColumns',
                             'Buttons'),
              options = list(
                pageLength = 10,
                scrollX = TRUE,
                scrollCollapse = TRUE,
                dom = 'Bfrtip',
                buttons = c('copy',
                            'csv',
                            'excel')
              ))

```

A closer look at HSV transcript presence in the data.

```{r matthelp_assaydata, dependson = c("umap","downsample","matthelp_hsvgenes")}
#Make a compressed matrix of HSV1 genes as rownames and the corresponding single cell data for each gene (getassaydata is scaled data)
assay.data <- GetAssayData(A673_hsv_even[hsv])

# Rownames of compressed matrix (assay.data) should be the genes from the subset of genes above (hsv <- str_subset)
#row.names(assay.data)

# Make a numeric vector of the column sums 
assay_col <- colSums(assay.data)

assay_col_unscaled <- GetAssayData(A673_hsv_even[hsv], ) %>%
  colSums()
(assay_col > 0) %>%
  summary()
```

Distribution of HSV1 transcripts:

```{r tibble_hist, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata")}

#Use tibble and ggplot to look at histogram distribution of HSV genes 
tibble(cell_sum = assay_col) %>%
  ggplot(aes(cell_sum)) +
  geom_histogram(bins = 200)
```

Log normalized distribution of HSV1 transcripts: 
```{r tibble_hist_log, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata")}

#log scale helps to separate out low expression values and zero. We add 0.000001 to make sure we don't have log(0) (equals infinity so algorithm will remove these values) but can still view all the data points. 
tibble(cell_sum = assay_col) %>%
  ggplot(aes(cell_sum + 0.000001)) +
  geom_histogram(bins = 200) +
  scale_x_log10()
```

Given that this new analysis better captures viral presence, it is added to the tumor data as metadata. The UMAP distribution of oHSV transcripts is shown. 
```{r new_hsv_meta, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata")}
A673_hsv_even <- AddMetaData(A673_hsv_even,
                               assay_col,
                               col.name = "oHSV_abundance")

A673_hsv_even <- AddMetaData(A673_hsv_even,
                          assay_col > 0,
                          col.name = "oHSV_presence")
```

### DimPlots: oHSV presence

```{r hsvdimplot, fig.height=5, fig.width=10, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}

Idents(A673_hsv_even) <- A673_hsv_even$oHSV_presence

# Change identities from TRUE/FALSE to infected/uninfected 
A673_hsv_even <- RenameIdents(A673_hsv_even, 
                                c("TRUE" = "oHSV+", 
                                  "FALSE" = "oHSV-"))

DimPlot(A673_hsv_even,
        pt.size = 0.5,
        order = T) +
  ggtitle("oHSV Transcript Presence") +
  coord_fixed()

ohsvpresence <- DimPlot(A673_hsv_even,
                        pt.size = 0.5,
                        split.by = "treatment",
                        order = TRUE) +
  ggtitle("oHSV Transcript Presence") +
  coord_fixed()
ohsvpresence

pdf("Plots/A673_1/ohsv_presence_dimplot.pdf",
    height = 5, width = 10)
ohsvpresence
dev.off()
```

### FeaturePlots: oHSV abundance 

```{r featurehsv, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
FeaturePlot(A673_hsv_even,
            features = "oHSV_abundance",
            order = TRUE) +
  theme(aspect.ratio = 1)
```

**Split by treatment group**
```{r feature_trt_hsv, fig.height = 3, fig.width = 12, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
FeaturePlot(A673_hsv_even,
            split.by="treatment",
            features = "oHSV_abundance",
            order = TRUE) +
  theme(aspect.ratio = 1)
```


**Max Cutoff at 100**

```{r feature_trt_hsv_100, fig.height = 3, fig.width = 12, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
# Max cutoff at 100
FeaturePlot(A673_hsv_even,
            split.by="treatment",
            features = "oHSV_abundance",
            order = TRUE,
            max.cutoff = 100) +
  theme(aspect.ratio = 1)
```


**Max Cutoff at 50**

```{r feature_trt_hsv_50, fig.height = 3, fig.width = 12, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
# Max cutoff at 50
FeaturePlot(A673_hsv_even,
            split.by="treatment",
            features = "oHSV_abundance",
            order = TRUE,
            max.cutoff = 50) +
  theme(aspect.ratio = 1)
```


**Max Cutoff at 20**

```{r feature_trt_hsv_20, fig.height = 3, fig.width = 12, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
# Max cutoff at 20
FeaturePlot(A673_hsv_even,
            split.by="treatment",
            features = "oHSV_abundance",
            order = TRUE,
            max.cutoff = 20) +
  theme(aspect.ratio = 1)
```

### Violin Plots: oHSV abundance

oHSV abundance is shown below across treatment groups, as well as by cluster and treatment groups. 

```{r hsv_ab, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
Idents(A673_hsv_even) <- A673_hsv_even$treatment

# Violin plots of HSV abundance
VlnPlot(A673_hsv_even,
        features = "oHSV_abundance",
        split.by = "treatment")

obj4vln <- RenameIdents(A673_hsv_even,
                        c("1 - Control" = "1",
                          "2 - oHSV" = "2",
                          "3 - Trabectedin" = "3",
                          "4 - oHSV+Trabectedin" = "4"))
log_ab <- VlnPlot(obj4vln,
                  features = "oHSV_abundance",
                  split.by = "treatment",
                  log = TRUE) +
  NoLegend()

log_ab
# save as pdf
pdf("Plots/A673_1/vln_log_hsvabundance.pdf",
    height = 5, width = 5)
print(log_ab)
dev.off()

pdf("Plots/A673_1/vln_log_hsvabundance2.pdf",
    height = 2, width = 8)
print(log_ab)
dev.off()

VlnPlot(A673_hsv_even,
        features = "oHSV_abundance",
        split.by = "treatment",
        group.by = "seurat_clusters")

# Histogram of HSV abundance per treatment group
for (trtmnt in levels(A673_hsv_even)) {
  ab_hist <- subset(A673_hsv_even, idents = trtmnt) %>%
    feature_hist(features = "oHSV_abundance",
                 cutoff_table = NULL) +
    ggtitle("oHSV Abundance",
            subtitle = paste(trtmnt))
  # Print in Rmd
  print(ab_hist)
  # Save as pdf
  pdf(paste0("Plots/A673_1/abundancehist_", trtmnt, ".pdf"),
      height = 5, width = 5)
  print(ab_hist)
  dev.off()
}
```

### HSV Time Point Module Scores 

Look at immediate early, early, and late gene expression and add to the data as a module score. Below the baseline plot,
cutoffs are used to better visualize the spread of oHSV transcripts for each stage of gene expression. Ensuring that the expression patterns 
are not exactly the same supports that presence of oHSV transcripts is not due to random spread during 10x sample preparation. 

Immediate early: "HSV1-RL2",
                 "HSV1-RS1",
                 "HSV1-UL54",
                 "HSV1-US1",
                 "HSV1-US12"

Early: "HSV1-UL23",
       "HSV1-UL29",
       "HSV1-UL50",
       "HSV1-UL2"

Late: "HSV1-UL48",
      "HSV1-UL19",
      "HSV1-US6" (note: not present in our data), 
      "HSV1-UL27",
      "HSV1-UL53",
      "HSV1-UL44",
      "HSV1-UL41"


```{r hsv_timepoints, fig.height = 4, fig.width = 16, dependson = c("umap","downsample","matthelp_hsvgenes","matthelp_assaydata","new_hsv_meta")}
# Create lists of time points for HSV gene expression 
imm_early <- list(c("HSV1-RL2",
                    "HSV1-RS1",
                    "HSV1-UL54",
                    "HSV1-US1",
                    "HSV1-US12"))

early <- list(c("HSV1-UL23",
                "HSV1-UL29",
                "HSV1-UL50",
                "HSV1-UL2"))

late <- list(c("HSV1-UL48",
               "HSV1-UL19",
               "HSV1-US6",
               "HSV1-UL27",
               "HSV1-UL53",
               "HSV1-UL44",
               "HSV1-UL41"))

# Add HSV timepoint module scores 
timepoints <- c("Immediate_Early_Genes",
                "Early_Genes",
                "Late_Genes")

hsv_timepoints <- c(imm_early,
                    early,
                    late)

names(hsv_timepoints) <- timepoints

plots <- list()
for (hsv_time in names(hsv_timepoints)) {
  # Add timepoint data as module score
  A673_hsv_even <- AddModuleScore(A673_hsv_even,
                                  features = list(hsv_timepoints[[hsv_time]]),
                                  ctrl = 5,
                                  name = hsv_time)
  
  # Rename AddModuleScore naming syntax
  colnames(A673_hsv_even@meta.data)[which(names(A673_hsv_even@meta.data) == paste0(hsv_time, "1"))] <- hsv_time
  
  # Create FeaturePlots based on the timepoint module scores
  plots[[paste0(hsv_time, "_plot")]] <- 
    FeaturePlot(A673_hsv_even,
                features = hsv_time,
                split.by = "treatment",
                order = TRUE,
                min.cutoff = 0) +
    theme(aspect.ratio = 1)

  print(plots[[paste0(hsv_time, "_plot")]])

# Create FeaturePlots based on the timepoint module scores
  plots[[paste0(hsv_time, "_plot_cutoff15")]] <- 
    FeaturePlot(A673_hsv_even,
                features = hsv_time,
                split.by = "treatment",
                order = TRUE,
                max.cutoff = 15,
                min.cutoff = 0) +
    theme(aspect.ratio = 1)
  
  plots[[paste0(hsv_time, "_plot_cutoff1")]] <- 
    FeaturePlot(A673_hsv_even,
                features = hsv_time,
                split.by = "treatment",
                order = TRUE,
                max.cutoff = 1,
                min.cutoff = 0) +
    theme(aspect.ratio = 1)
  
  plots[[paste0(hsv_time, "rfeatureplot_plot_cutoff1")]] <- 
    r_feature_plot(A673_hsv_even,
                   features = hsv_time,
                   split.by = "treatment",
                   order = TRUE,
                   max.cutoff = 1,
                   min.cutoff = 0) +
    theme(aspect.ratio = 1)
  
  print(plots[[paste0(hsv_time, "rfeatureplot_plot_cutoff1")]])
  
  pdf(paste0("Plots/A673_1/rfeaturecutoff1_", hsv_time, trtmnt),
      height = 5, width = 10)
  print(plots[[paste0(hsv_time, "rfeatureplot_plot_cutoff1")]])
  dev.off()
}

save(A673_hsv_even, file  = "Data/tumor-only-A673_hsv_even-labeled.RData")
```

```{r hsv_patchwork_cutoff15, fig.height = 12, fig.width = 16}
library(patchwork)

hsv_patchwork_hi <- 
  (plots[["Immediate_Early_Genes_plot_cutoff15"]] / 
     plots[["Early_Genes_plot_cutoff15"]] / 
     plots[["Late_Genes_plot_cutoff15"]]) + 
  plot_annotation(title = 'Max Cutoff of 15',
                  theme = theme(plot.title = element_text(size = 18)))

hsv_patchwork_hi

```

```{r hsv_patchwork_cutoff1, fig.height = 12, fig.width = 16}
library(patchwork)

hsv_patchwork_lo <- 
  (plots[["Immediate_Early_Genes_plot_cutoff1"]] / 
     plots[["Early_Genes_plot_cutoff1"]] / 
     plots[["Late_Genes_plot_cutoff1"]]) +
  plot_annotation(title = 'Max Cutoff of 1',
                  theme = theme(plot.title = element_text(size = 18)))

hsv_patchwork_lo

```

```{r hsv_vln, eval = F}
for (hsv_time in names(hsv_timepoints)) {
# Create Violin Plots based on the timepoint module scores
  hsv_timepoints[["plots"]][[paste0(hsv_time, "_vln")]] <- 
    VlnPlot(A673_hsv_even,
            features = hsv_time,
            group.by = "treatment",
            log = TRUE) +
    theme(aspect.ratio = 1)
  
print(hsv_timepoints[["plots"]][[paste0(hsv_time, "_vln")]])
}

hsv_patchwork_vln <- 
  (hsv_timepoints[["plots"]][["Immediate_Early_Genes_vln"]] / 
     hsv_timepoints[["plots"]][["Early_Genes_vln"]] / 
     hsv_timepoints[["plots"]][["Late_Genes_vln"]]) + 
  plot_annotation(title = 'Violin Plot of HSV Timepoint Expression',
                  theme = theme(plot.title = element_text(size = 18)))

hsv_patchwork_vln
```


```{r hsv-vln, dependson="hsv_timepoints"}
Idents(A673_hsv_even) <- A673_hsv_even$treatment

VlnPlot(A673_hsv_even,
        features = "Immediate_Early_Genes",
        split.by = "treatment")

VlnPlot(A673_hsv_even,
        features = "Early_Genes",
        split.by = "treatment")

VlnPlot(A673_hsv_even,
        features = "Late_Genes",
        split.by = "treatment")

```

```{r hsv_hist, dependson="hsv_timepoints"}
for (trtmnt in c("2 - oHSV", "4 - oHSV+Trabectedin")) {
  for (hsv_time in names(hsv_timepoints)) {
    # Create histogram Plots based on the timepoint module scores
    plots[[paste0(hsv_time, "_hist")]][[trtmnt]] <-
      subset(A673_hsv_even,
             idents = trtmnt) %>%
      feature_hist(features = hsv_time,
                   cutoff_table = NULL) +
      ggtitle(paste(trtmnt),
              subtitle = hsv_time)
    # Print in Rmd
    print(plots[[paste0(hsv_time, "_hist")]][[trtmnt]])
    # Save as pdf
    pdf(paste0("Plots/A673_1/hist_", hsv_time, "_", trtmnt, ".pdf"),
        height = 5, width = 5)
    print(plots[[paste0(hsv_time, "_hist")]][[trtmnt]])
    dev.off()
  }
}

```

```{r ridge, dependson="hsv_timepoints"}
for (hsv_time in names(hsv_timepoints)) {
  Idents(A673_hsv_even) <- A673_hsv_even$oHSV_presence 
  sobj <- subset(A673_hsv_even,
                 idents = "TRUE")
  Idents(sobj) <- sobj$treatment
  sobj <- subset(sobj,
                 idents = c("2 - oHSV", "4 - oHSV+Trabectedin"))
  ridge <- RidgePlot(sobj,
                     features = hsv_time,
                     y.max = 10,
                     same.y.lims = T,
                     log = F,
                     sort = T)
  print(ridge)
  
  pdf(paste0("Plots/A673_1/ridge", hsv_time, ".pdf"),
      height = 5, width = 10)
  print(ridge)
  dev.off()
}
```

```{r ridge2}
for (trtmnt in c("2 - oHSV", "4 - oHSV+Trabectedin")) {
     for (hsv_time in names(hsv_timepoints)) {
  Idents(A673_hsv_even) <- A673_hsv_even$oHSV_presence 
  sobj <- subset(A673_hsv_even,
                 idents = "TRUE")
  Idents(sobj) <- sobj$treatment
  sobj <- subset(sobj,
                 idents = trtmnt)
  sobj <- RenameIdents(sobj,
                       "2 - oHSV" = "2",
                       "4 - oHSV+Trabectedin" = "4")
  ridge <- RidgePlot(sobj,
                     features = hsv_time,
                     y.max = 10,
                     log = T,
                     sort = T) +
    NoLegend()
  print(ridge)
  
  pdf(paste0("Plots/A673_1/ridge_p2", hsv_time, trtmnt, ".pdf"),
      height = 5, width = 10)
  print(ridge)
  dev.off()
     }
}
```


```{r height = 10, width = 10}
Idents(A673_hsv_even) <- A673_hsv_even$treatment
dot <- subset(A673_hsv_even,
              idents = c("2 - oHSV", "4 - oHSV+Trabectedin")) %>%
  DotPlot(features = c("Immediate_Early_Genes",
                       "Early_Genes",
                       "Late_Genes"),
          scale = F,
          group.by = "treatment",
          cols = c("lightgoldenrod", "darkred"),
          dot.min = 0)
print(dot)

dot2 <- subset(A673_hsv_even,
              idents = c("2 - oHSV", "4 - oHSV+Trabectedin")) %>%
  DotPlot(features = c("Immediate_Early_Genes",
                       "Early_Genes",
                       "Late_Genes"),
          scale = F,
          group.by = "treatment",
          cols = c("lightgoldenrod", "darkred"),
          dot.scale = 20,
          dot.min = -10)
print(dot2)


feature <- subset(A673_hsv_even,
                  idents = c("2 - oHSV", "4 - oHSV+Trabectedin")) %>%
  r_feature_plot(features = c("Immediate_Early_Genes",
                              "Early_Genes",
                              "Late_Genes"),
                 split.by = "treatment",
                 order = TRUE,
                 max.cutoff = 1,
                 min.cutoff = 0) +
  theme(aspect.ratio = 1)

print(feature)
pdf("Plots/A673_1/featureplot_time.pdf", height = 12.5, width = 8.25)
print(feature)
dev.off()

vln <- subset(A673_hsv_even,
              idents = c("2 - oHSV", "4 - oHSV+Trabectedin")) %>%
  VlnPlot(features = c("Immediate_Early_Genes",
                       "Early_Genes",
                       "Late_Genes"),
          group.by = "treatment",
          same.y.lims = T,
          log = T) +
  NoLegend()
print(vln)

pdf("Plots/A673_1/vlnplot_time.pdf", height = 12.5, width = 8.25)
print(feature)
dev.off()
```

```{r dot_ab, height = 2, width = 10}
DotPlot(A673_hsv_even,
        features = "oHSV_abundance",
        scale = F,
        group.by = "treatment",
        cols = c("lightgoldenrod", "darkred")
)
```


# Antiviral Analysis

## Gene List

```{r}
DotPlot(A673_hsv_even,
        features = c(
          "ISG15",
          "IFIT2",
          "ZC3HAV1",
          "MAFA",
          "DDX58",
          "IFIT3",
          "HELZ2",
          "IFITM3"),
        cols = c("lightgoldenrod", "darkred"),
        scale = FALSE)
```


## Differential Gene Expression

Differential gene expression in combination compared to oHSV alone (gene list shown in the table below). GO and KEGG analyses are performed.

```{r deg-trt}
trtdiff <- FindMarkers(A673_hsv_even,
                       ident.1 = "4 - oHSV+Trabectedin",
                       ident.2 = "2 - oHSV",
                       logfc.threshold = 0.15,
                       only.pos = F) %>%
  arrange(-abs(avg_log2FC))

trtdiff <- as.data.frame(trtdiff)
trtdiff <- arrange(trtdiff, desc(avg_log2FC))
degs <- as.vector(trtdiff$avg_log2FC)
names(degs) <- rownames(trtdiff) %>%
  str_replace("HSV1-","")
degs

trt <- list()
trt[["GO-BP"]] <- clusterProfiler::gseGO(
  geneList = degs,
  OrgDb = org.Hs.eg.db::org.Hs.eg.db,
  ont = "BP",
  keyType = "SYMBOL",
  nPermSimple = 10000,
  eps = 0)

trt[["GO-BP"]] <- mutate(trt[["GO-BP"]],
      p.adjust = -log10(p.adjust))

plots <- list()
if (nrow(trt[["GO-BP"]]) > 0) {
  plots[["GO-BP"]] <- enrichplot::dotplot(trt[["GO-BP"]],
                                          x = "NES",
                                          showCategory = 30,
                                          font.size = 8) +
    labs(title = "Upregulated Gene Sets in Combination Therapy",
         subtitle = "GO-BP Gene Sets")
} else {
  plots[["GO-BP"]] <-
    ggplot(tibble(x = "A", y = "A",
                  text = "No significant hits"),
           aes(x = x, y = y, label = text)) +
    geom_text() +
    theme(axis.text = element_text(size = 5)) +
    labs(title = "Upregulated Gene Sets in Combination Therapy",
         subtitle = "GO-BP Gene Sets")
}

plots[["GO-BP"]]

# enrichplot::dotplot(trt[[1]][["GO-BP"]],
#                     x = "NES",
#                     showCategory = 15,
#                     font.size = 9) +
#   labs(title = paste(t, c),
#        subtitle = "GO-BP Gene Sets")
# 
# enrichplot::dotplot(trt[[2]][["GO-BP"]],
#                     x = "NES",
#                     showCategory = 15,
#                     font.size = 9) +
#   labs(title = paste(t, c),
#        subtitle = "GO-BP Gene Sets")

```

```{r kegg, dependson='deg-trt', eval = F}
trt[["KEGG"]] <- clusterProfiler::gseKEGG(
  geneList = degs,
  organism = c("hsa"),
  keyType = "kegg",
  eps = 0)

trt[["KEGG"]] <- mutate(trt[["KEGG"]],
      p.adjust = -log10(p.adjust))

plots <- list()
if (nrow(trt[["KEGG"]]) > 0) {
  plots[["KEGG"]] <- enrichplot::dotplot(trt[["KEGG"]],
                                         x = "NES",
                                         showCategory = 15,
                                         font.size = 8) +
    labs(title = "Upregulated Pathways in Combination Therapy",
         subtitle = "KEGG Gene Sets")
} else {
  plots[["KEGG"]] <-
    ggplot(tibble(x = "A", y = "A",
                  text = "No significant hits"),
           aes(x = x, y = y, label = text)) +
    geom_text() +
    theme(axis.text = element_text(size = 5)) +
    labs(title = "Upregulated Pathways in Combination Therapy",
         subtitle = "KEGG Gene Sets")
}

plots[["KEGG"]]
```


```{r trtdiff, dependson="deg-trt"}
# Create interactive datatable showing DEG list
#trtdiff <- as.data.table(trtdiff)
DT::datatable(trtdiff, 
              rownames = TRUE,
              extensions = c('FixedColumns',
                             'Buttons'),
              options = list(
                pageLength = 10,
                scrollX = TRUE,
                scrollCollapse = TRUE,
                dom = 'Bfrtip',
                buttons = c('copy',
                            'csv',
                            'excel')
              ))
```


```{r gobp, dependson="deg-trt"}
# Create interactive datatable showing DEG list
gobp <- as.data.table(trt[["GO-BP"]])
DT::datatable(gobp,
              rownames = TRUE,
              extensions = c('FixedColumns',
                             'Buttons'),
              options = list(
                pageLength = 10,
                scrollX = TRUE,
                scrollCollapse = TRUE,
                dom = 'Bfrtip',
                buttons = c('copy',
                            'csv',
                            'excel')
              ))
```


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