CPTR5_initial_qc_markdown.Rmd

---
title: "DSP Workflow Test Package"
output: html_document
date: "2023-10-24"
---

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

# Increase the time out for downloading the DSP package
options(timeout = max(300, getOption("timeout")))
```

## Install DSP Workflow from github

```{r}

# Install DSPWorkflow package

install.DSP <- FALSE
if(install.DSP == TRUE){
  library(devtools)
  install_github("NIDAP-Community/DSPWorkflow", ref = "dev")
}

library(DSPWorkflow)

```

## Load input data

``` {r Load Data, echo=TRUE}
# Load all inputs

dcc.files <- list.files(file.path("dcc"),
  pattern = ".dcc$",
  full.names = TRUE,
  recursive = TRUE
)

pkc.files <- "Hs_R_NGS_WTA_v1.0.pkc"
pheno.data.file <- "annotation_Krug_CPTR_5_april2024_NC_edit.xlsx"

```

# Run the first step of the DSP Workflow package

# Study Design

```{r Study Design, echo=TRUE}

# Save the output from the study design function into a list
sdesign.list <- studyDesign(dcc.files = dcc.files, 
                                pkc.files = pkc.files,
                                pheno.data.file = pheno.data.file,
                                pheno.data.sheet = "annotation",
                                pheno.data.dcc.col.name = "Sample_ID",
                                protocol.data.col.names = c("ROI"),
                                experiment.data.col.names = c("panel"),
                                slide.name.col = "slide name", 
                                class.col = "class", 
                                region.col = "Region", 
                                segment.col = "segment",
                                area.col = "area",
                                nuclei.col = "nuclei", 
                                sankey.exclude.slide = FALSE, 
                                segment.id.length = 10)

# The output of the study design function is a Geomxset Object and a Plot
# Print out a summary of the object

print(sdesign.list$object)

# Print out the Sankey Plot

print(sdesign.list$sankey.plot)
```

## 2. QC Preprocessing:

```{r QC Preprocessing, echo=TRUE}

qc.output <-  qcProc(object = sdesign.list$object,
                        min.segment.reads = 1000, 
                        percent.trimmed = 80,    
                        percent.stitched = 80,   
                        percent.aligned = 80,    
                        percent.saturation = 50, 
                        min.negative.count = 3,   
                        max.ntc.count = 1000,     
                        min.nuclei = 200,         
                        min.area = 1000,
                        print.plots = TRUE)
    print(qc.output$segments.qc)
    print(qc.output$segment.flags)
    print(qc.output$probe.flags)
    
    
  # Export the flags table
    
  export.flags <- FALSE
  
  if(export.flags == TRUE){
    
    write.csv(qc.output$segment.flags, file =  "qc/segment_qc_flags.csv")
    
    write.csv(qc.output$probe.flags, file =  "qc/probe_qc_flags.csv")
    
  }
    
```

## 3. Filtering

### Segment Filtering by Gene Detection

```{r}

library(GeomxTools)
library(dplyr)
library(knitr)


object <- qc.output$object

# Set up lists of segment IDs
segment.list.total <- pData(object)$segmentID

# Define Modules
modules <- gsub(".pkc", "", pkc.files)

# Calculate limit of quantification (LOQ) in each segment
# LOQ = geomean(NegProbes) * geoSD(NegProbes)^(LOQ cutoff)
# LOQ is calculated for each module (pkc file)
loq <- data.frame(row.names = colnames(object))

loq.min <- 2
loq.cutoff <- 2

for(module in modules) {
  vars <- paste0(c("NegGeoMean_", "NegGeoSD_"),
                 module)
  if(all(vars[1:2] %in% colnames(pData(object)))) {
    
    neg.geo.mean <- vars[1]
    neg.geo.sd <- vars[2]
    
    loq[, module] <-
      pmax(loq.min,
           pData(object)[, neg.geo.mean] * 
             pData(object)[, neg.geo.sd] ^ loq.cutoff)
  }
}

# Store the loq df in the annotation df
pData(object)$loq <- loq

# Setup a master loq matrix
loq.mat <- c()


for(module in modules) {
  # Gather rows with the given module
  ind <- fData(object)$Module == module
  
  # Check if each feature has counts above the LOQ
  mat.i <- t(esApply(object[ind, ], MARGIN = 1,
                     FUN = function(x) {
                       x > loq[, module]
                     }))
  
  # Store results in the master loq matrix
  loq.mat <- rbind(loq.mat, mat.i)
}

# ensure ordering since this is stored outside of the geomxSet
loq.mat <- loq.mat[fData(object)$TargetName, ]

# Evaluate and Filter Segment Gene Detection Rate
# Save detection rate information to pheno data
pData(object)$GenesDetected <- colSums(loq.mat, na.rm = TRUE)
pData(object)$GeneDetectionRate <- 100*(pData(object)$GenesDetected / nrow(object))

# Establish detection bins
detection.bins <- c("less_than_1", "1_5", "5_10", "10_15", "greater_than_15")

# Determine detection thresholds: 1%, 5%, 10%, 15%, >15%
pData(object)$DetectionThreshold <- 
  cut(pData(object)$GeneDetectionRate,
      breaks = c(0, 1, 5, 10, 15, 100),
      labels = detection.bins)

# stacked bar plot of different cut points (1%, 5%, 10%, 15%)
segment.stacked.bar.plot<- ggplot(pData(object),
                          aes(x = DetectionThreshold)) +
  geom_bar(aes(fill = region)) +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  theme_bw() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "Gene Detection Rate",
         y = "Segments, #",
         fill = "Segment Type")
  
# cut percent genes detected at 1, 5, 10, 15
segment.table <- kable(table(pData(object)$DetectionThreshold, 
                             pData(object)$class))

# Make a list of segments with low detection
low.detection.segments <- pData(object) %>% 
  filter(GeneDetectionRate < 5) %>% 
  select(any_of(c("segmentID", "GeneDetectionRate")))

print(low.detection.segments)

# Export a summary of the segment gene detection
segment.detection.summary <- pData(object) %>% 
  select(any_of(c("segmentID", "GeneDetectionRate", "DetectionThreshold")))

export.segment.detection.summary <- FALSE

if(export.segment.detection.summary == TRUE){
  
  write.csv(segment.detection.summary, "qc/segment_detection_summary.csv")
  
}

```

### Remove Segments

```{r}

# Filter the data using the cutoff for gene detection rate
segment.gene.rate.cutoff <- 1

object.segment.filtered <-
    object[, pData(object)$GeneDetectionRate >= segment.gene.rate.cutoff]


```

### Gene Filtering by Detection per Segment

```{r}
library(scales)

# Evaluate and Filter Study-wide Gene Detection Rate 
# Calculate detection rate:
loq.mat <- loq.mat[, colnames(object.segment.filtered)]

fData(object.segment.filtered)$DetectedSegments <- rowSums(loq.mat, na.rm = TRUE)
fData(object.segment.filtered)$DetectionRate <-
  100*(fData(object.segment.filtered)$DetectedSegments / nrow(pData(object)))

# Establish detection bins
detection.bins <- c("0", "less_than_1", "1_5", "5_10", "10_20", "20_30", "30_40", "40_50", "greater_than_50")

# Determine detection thresholds: 1%, 5%, 10%, 15%, >15%
fData(object.segment.filtered)$DetectionThreshold <- 
  cut(fData(object.segment.filtered)$DetectionRate,
      breaks = c(-1, 0, 1, 5, 10, 20, 30, 40, 50, 100),
      labels = detection.bins)


gene.stacked.bar.plot <- ggplot(fData(object.segment.filtered),
                          aes(x = DetectionThreshold)) +
  geom_bar(aes(fill = Module)) +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5) +
  theme_bw() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "Gene Detection Rate",
         y = "Genes, #",
         fill = "Probe Set")


# Gene of interest detection table
goi <- c("A2M", "CD44")

goi.table <- data.frame(Gene = goi,
                        Number = fData(object.segment.filtered)[goi, "DetectedSegments"],
                        DetectionRate = percent(fData(object.segment.filtered)[goi, "DetectionRate"]))
print(goi.table)

# Plot detection rate:
plot.detect <- data.frame(Freq = c(1, 5, 10, 20, 30, 50))
plot.detect$Number <-
  unlist(lapply(c(1, 5, 10, 20, 30, 50),
                function(x) {sum(fData(object.segment.filtered)$DetectionRate >= x)}))

plot.detect$Rate <- plot.detect$Number / nrow(fData(object.segment.filtered))
rownames(plot.detect) <- plot.detect$Freq

genes.detected.plot <- ggplot(plot.detect, aes(x = as.factor(Freq), y = Rate, fill = Rate)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = formatC(Number, format = "d", big.mark = ",")),
            vjust = 1.6, color = "black", size = 4) +
  scale_fill_gradient2(low = "orange2", mid = "lightblue",
                       high = "dodgerblue3", midpoint = 0.65,
                       limits = c(0,1),
                       labels = scales::percent) +
  theme_bw() +
  scale_y_continuous(labels = scales::percent, limits = c(0,1),
                     expand = expansion(mult = c(0, 0))) +
  labs(x = "% of Segments",
       y = "Genes Detected, % of Panel > loq")


# Export a summary of the gene detection
gene.detection.summary <- fData(object.segment.filtered) %>% 
  select(any_of(c("segmentID", "DetectionRate", "DetectionThreshold")))

export.gene.detection.summary <- FALSE

if(export.gene.detection.summary == TRUE){
  
  write.csv(gene.detection.summary, "qc/gene_detection_summary.csv")
  
}

```

### Remove Genes

```{r}

# Set the cutoff for gene detection
study.gene.rate.cutoff <- 0.01

# Subset for genes above the study gene detection rate cutoff
# Manually include the negative control probe, for downstream use
negative.probe.fData <- subset(fData(object.segment.filtered), CodeClass == "Negative")
neg.probes <- unique(negative.probe.fData$TargetName)
object.gene.filtered <- object.segment.filtered[fData(object.segment.filtered)$DetectionRate >= study.gene.rate.cutoff |
                   fData(object.segment.filtered)$TargetName %in% neg.probes, ]

```

## 4. Normalization:
  
```{r Normalization, echo=TRUE}
  
    q3.normalization.output <- geomxNorm(
                                  object = object.gene.filtered, 
                                  norm = "q3")
    
    print(q3.normalization.output$multi.plot)
    print(q3.normalization.output$boxplot.raw)
    print(q3.normalization.output$boxplot.norm)
    
    neg.normalization.output <- geomxNorm(
                                  object = object.gene.filtered, 
                                  norm = "neg")
    
    print(neg.normalization.output$multi.plot)
    print(neg.normalization.output$boxplot.raw)
    print(neg.normalization.output$boxplot.norm)
```


## 5. Unsupervised Analysis:

```{r Unsupervised Analysis, echo=TRUE}

    #Test Unsupervised Analysis:
    unsupervised.output <- dimReduct(object = q3.normalization.output$object,
                        point.size = 2,
                        point.alpha = 0.5,
                        color.variable1 = "region",
                        shape.variable = "class"
    )
    
    print(unsupervised.output$plot$PCA)
    print(unsupervised.output$plot$tSNE)
    print(unsupervised.output$plot$UMAP)

```


## 6. Clustering high CV Genes and Heatmap:

 
```{r Clustering high CV Genes, echo=TRUE}
    
  anno.colors <- list(
    region = c(immune = "red", 
               stroma = "seagreen1", 
               'stroma/immune' = "olivedrab", 
               tumor = "darkblue", 
               vessel = "purple"), 
    class = c(Biopsy = "orange", 
              PDX = "springgreen")
  )  

  heatmap.output <- heatMap(object = unsupervised.output$object, 
                        ngenes = 200, 
                        scale.by.row.or.col = "row", 
                        show.rownames = FALSE, 
                        show.colnames = FALSE, 
                        clustering.method = "average", 
                        cluster.rows = TRUE, 
                        cluster.cols = TRUE,
                        clustering.distance.rows = "correlation", 
                        clustering.distance.cols = "correlation", 
                        annotation.row = NA, 
                        annotation.col = c("class", "region"), 
                        breaks.by.values = seq(-3, 3, 0.05), 
                        heatmap.color = colorRampPalette(c("blue", "white", "red"))(120), 
                        norm.method = "quant", 
                        annotation.colors = anno.colors)


     print(heatmap.output$plot)
    
```
# 7. Differential Expression

```{r}

library(GeomxTools)
library(dplyr)

# Source the helper functions
#source("DSP_functions.R")

```

## Biopsy vs. Biopsy

### Infected vs. Uninfected

#### Setup

```{r, warning=FALSE}

# Setup annotation groups for the contrast
tumor.types <- c("Biopsy")
region.types <- c("immune", "tumor", "vessel")
source.types <- c("Skin", "GI")
segment.types <- c("FullROI")

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list(
  "Tumor" = tumor.types, 
  "region" = region.types, 
  "segment" = segment.types, 
  "Source" = source.types)

# Filter data for the chosen annotations
object <- q3.normalization.output$object

subset.indices <- pData(object)$Tumor %in% tumor.types & 
  pData(object)$region %in% region.types &
  pData(object)$segment %in% segment.types &
  pData(object)$Source %in% source.types

subset.object <- object[, subset.indices]

# Convert all annotations to factors
for(column.name in names(contrast.groups.list)){ 
  pData(subset.object)[[column.name]] <- 
    factor(pData(subset.object)[[column.name]], 
           contrast.groups.list[[column.name]])
}

# Factor the slide column
pData(subset.object)[["slide_name"]] <- 
  factor(pData(subset.object)[["slide_name"]])

# Create log2 counts
assayDataElement(object = subset.object, elt = "log_q") <-
    assayDataApply(subset.object, 2, FUN = log, base = 2, elt = "q_norm")

# Edit columns with bad characters

pData(subset.object)$kshv_infection <- pData(subset.object)$'infected with KSHV'

```

#### DE contrast: Skin and GI, Tumor (uninfected) & Vessel (infected)

```{r}

# Select the column with levels for the contrast
contrast.column <- "region"
contrast.levels <- c("tumor", "vessel")

# Further subset if needed
subset.column <- "region"
subset.levels <- c("tumor", "vessel")

contrast.indices <- pData(subset.object)[[subset.column]] %in% subset.levels

contrast.object <- subset.object[, contrast.indices]

# Create summary table of group counts
summary.table.df <- pData(contrast.object) %>%  
  select(c(names(contrast.groups.list), contrast.column))

summary.table <- table(summary.table.df)

print("Sample Numbers per Annotation Group")
print("-----------------------------------")
print(summary.table)

# Choose groups to make the contrast within

column.name <- "Source"
groups.list <- list()
groups.list[[column.name]] <- c("Skin", "GI")

# Setup a vector for the results
results <- c()

# Run DE analysis for each group
for(group in groups.list[[column.name]]){ 
  
  
  lmm.results.list <- run_lmm()
  
  # Subset for the current source
  group.indices <- pData(contrast.object)[[column.name]] == group
  subset.group.object <- contrast.object[, group.indices]
  
  # Run the linear model
  lm.results <- mixedModelDE(subset.group.object, 
                             elt = "log_q", 
                             modelFormula = ~ region + (1 | slide_name), 
                             groupVar = contrast.column, 
                             nCores = parallel::detectCores(), 
                             multiCore = TRUE)
  
  # Gather the results into an output table
  r_test <- do.call(rbind, lm.results["lsmeans", ])
  tests <- rownames(r_test)
  r_test <- as.data.frame(r_test)
  r_test$Contrast <- tests
  
  # use lapply in case you have multiple levels of your test factor to
  # correctly associate gene name with it's row in the results table
  r_test$Gene <- 
      unlist(lapply(colnames(lm.results),
                    rep, nrow(lm.results["lsmeans", ][[1]])))
  r_test$Subset <- group
  r_test$FDR <- p.adjust(r_test$`Pr(>|t|)`, method = "fdr")
  r_test <- r_test[, c("Gene", "Subset", "Contrast", "Estimate", 
                       "Pr(>|t|)", "FDR")]
  results <- rbind(results, r_test)
  
}




```

##### Volcano Plot

```{r}

library(ggplot2)
library(ggrepel)
library(plotly)


# Create a column for direction of DEGs
results$de_direction <- "NONE"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate > 1] <- "UP"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate < -1] <- "DOWN"

results$deglabel <- ifelse(results$de_direction == "NONE", NA, results$Gene)

# Compute the scale for the volcano x-axis
log2.scale <- max(abs(results$Estimate))

# Make the volcano plot
volcano.plot <- ggplot(data = results, aes(x = Estimate, 
                                           y = -log10(FDR), 
                                           col = de_direction, 
                                           text = Gene)) +
  geom_vline(xintercept = c(-1, 1), col = "gray", linetype = 'dashed') +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = 'dashed') + 
  xlim(-7.5, 7.5) + 
  labs(x = "Tumor vs. Vessel", 
       y = "-log10 adjusted p-value", 
       title = paste0("Biopsy DE Analysis")) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("steelblue4", "grey", "violetred4"), 
                     labels = c("Downregulated in Tumor", "Not significant", "Upregulated in Tumor")) + 
  xlim(-log2.scale-1, log2.scale+1) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  facet_grid(cols = vars(Subset))


  #geom_text_repel(max.overlaps = Inf) +

# Plotly volcano

plotly.volcano <- ggplotly(volcano.plot)

export.volcano <- FALSE
if(export.volcano == TRUE){
    ggsave(paste0("Results/volcano_", contrast, "_adjpval.png"),
         plot = volcano.plot, 
         width = 8, 
         height = 8, 
         units = "in")
}

```

#### DE contrast: Skin and GI Immune, Infected vs Uninfected

```{r}

# Select the column with levels for the contrast
contrast.column <- "kshv_infection"
contrast.levels <- c("T", "F")

# Further subset if needed
subset.column <- "region"
subset.levels <- c("immune")

contrast.indices <- pData(subset.object)[[subset.column]] %in% subset.levels

contrast.object <- subset.object[, contrast.indices]

# Create summary table of group counts
summary.table.df <- pData(contrast.object) %>%  
  select(c(names(contrast.groups.list), contrast.column))

summary.table <- table(summary.table.df)

print("Sample Numbers per Annotation Group")
print("-----------------------------------")
print(summary.table)

# Choose groups to make the contrast within
group.column <- "Source"
group.names <- c("Skin", "GI")

# Setup a vector for the results
results <- c()

# Run DE analysis for each source
for(group in group.names){
  
  
  
  # Subset for the current source
  group.indices <- pData(contrast.object)[[group.column]] == group
  subset.group.object <- contrast.object[, group.indices]
  
  # Run the linear model
  lm.results <- mixedModelDE(subset.group.object, 
                             elt = "log_q", 
                             modelFormula = ~ kshv_infection + (1 | slide_name), 
                             groupVar = contrast.column, 
                             nCores = parallel::detectCores(), 
                             multiCore = TRUE)
  
  # Gather the results into an output table
  r_test <- do.call(rbind, lm.results["lsmeans", ])
  tests <- rownames(r_test)
  r_test <- as.data.frame(r_test)
  r_test$Contrast <- tests
  
  # use lapply in case you have multiple levels of your test factor to
  # correctly associate gene name with it's row in the results table
  r_test$Gene <- 
      unlist(lapply(colnames(lm.results),
                    rep, nrow(lm.results["lsmeans", ][[1]])))
  r_test$Subset <- group
  r_test$FDR <- p.adjust(r_test$`Pr(>|t|)`, method = "fdr")
  r_test <- r_test[, c("Gene", "Subset", "Contrast", "Estimate", 
                       "Pr(>|t|)", "FDR")]
  results <- rbind(results, r_test)
  
  }

```

##### Volcano Plot

```{r}

library(ggplot2)
library(ggrepel)
library(plotly)


# Create a column for direction of DEGs
results$de_direction <- "NONE"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate > 1] <- "UP"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate < -1] <- "DOWN"

results$deglabel <- ifelse(results$de_direction == "NONE", NA, results$Gene)

# Compute the scale for the volcano x-axis
log2.scale <- max(abs(results$Estimate))

# Make the volcano plot
volcano.plot <- ggplot(data = results, aes(x = Estimate, 
                                           y = -log10(FDR), 
                                           col = de_direction, 
                                           text = Gene)) +
  geom_vline(xintercept = c(-1, 1), col = "gray", linetype = 'dashed') +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = 'dashed') + 
  xlim(-7.5, 7.5) + 
  labs(x = "Tumor vs. Vessel", 
       y = "-log10 adjusted p-value", 
       title = paste0("Biopsy DE Analysis")) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("steelblue4", "grey", "violetred4"), 
                     labels = c("Downregulated in Tumor", "Not significant", "Upregulated in Tumor")) + 
  xlim(-log2.scale-1, log2.scale+1) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  facet_grid(cols = vars(Subset))


  #geom_text_repel(max.overlaps = Inf) +

# Plotly volcano

plotly.volcano <- ggplotly(volcano.plot)

export.volcano <- FALSE
if(export.volcano == TRUE){
    ggsave(paste0("Results/volcano_", contrast, "_adjpval.png"),
         plot = volcano.plot, 
         width = 8, 
         height = 8, 
         units = "in")
}

```


#### DE contrast: Skin and GI, Immune vs Tumor

```{r}

# Select the column with levels for the contrast
contrast.column <- "region"
contrast.levels <- c("tumor", "immune")

# Further subset if needed
subset.column <- "region"
subset.levels <- c("tumor", "immune")

contrast.indices <- pData(subset.object)[[subset.column]] %in% subset.levels

contrast.object <- subset.object[, contrast.indices]

# Create summary table of group counts
summary.table.df <- pData(contrast.object) %>%  
  select(c(names(contrast.groups.list), contrast.column))

summary.table <- table(summary.table.df)

print("Sample Numbers per Annotation Group")
print("-----------------------------------")
print(summary.table)

# Choose groups to make the contrast within
group.column <- "Source"
group.names <- c("Skin", "GI")

# Setup a vector for the results
results <- c()

# Run DE analysis for each source
for(group in group.names){
  
  
  
  # Subset for the current source
  group.indices <- pData(contrast.object)[[group.column]] == group
  subset.group.object <- contrast.object[, group.indices]
  
  # Run the linear model
  lm.results <- mixedModelDE(subset.group.object, 
                             elt = "log_q", 
                             modelFormula = ~ kshv_infection + (1 | slide_name), 
                             groupVar = contrast.column, 
                             nCores = parallel::detectCores(), 
                             multiCore = TRUE)
  
  # Gather the results into an output table
  r_test <- do.call(rbind, lm.results["lsmeans", ])
  tests <- rownames(r_test)
  r_test <- as.data.frame(r_test)
  r_test$Contrast <- tests
  
  # use lapply in case you have multiple levels of your test factor to
  # correctly associate gene name with it's row in the results table
  r_test$Gene <- 
      unlist(lapply(colnames(lm.results),
                    rep, nrow(lm.results["lsmeans", ][[1]])))
  r_test$Subset <- group
  r_test$FDR <- p.adjust(r_test$`Pr(>|t|)`, method = "fdr")
  r_test <- r_test[, c("Gene", "Subset", "Contrast", "Estimate", 
                       "Pr(>|t|)", "FDR")]
  results <- rbind(results, r_test)
  
  }

```

##### Volcano Plot

```{r}

library(ggplot2)
library(ggrepel)
library(plotly)


# Create a column for direction of DEGs
results$de_direction <- "NONE"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate > 1] <- "UP"
results$de_direction[results$FDR < 0.05 & 
                             results$Estimate < -1] <- "DOWN"

results$deglabel <- ifelse(results$de_direction == "NONE", NA, results$Gene)

# Compute the scale for the volcano x-axis
log2.scale <- max(abs(results$Estimate))

# Make the volcano plot
volcano.plot <- ggplot(data = results, aes(x = Estimate, 
                                           y = -log10(FDR), 
                                           col = de_direction, 
                                           text = Gene)) +
  geom_vline(xintercept = c(-1, 1), col = "gray", linetype = 'dashed') +
  geom_hline(yintercept = -log10(0.05), col = "gray", linetype = 'dashed') + 
  xlim(-7.5, 7.5) + 
  labs(x = "Tumor vs. Vessel", 
       y = "-log10 adjusted p-value", 
       title = paste0("Biopsy DE Analysis")) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("steelblue4", "grey", "violetred4"), 
                     labels = c("Downregulated in Tumor", "Not significant", "Upregulated in Tumor")) + 
  xlim(-log2.scale-1, log2.scale+1) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  facet_grid(cols = vars(Subset))


  #geom_text_repel(max.overlaps = Inf) +

# Plotly volcano

plotly.volcano <- ggplotly(volcano.plot)

export.volcano <- FALSE
if(export.volcano == TRUE){
    ggsave(paste0("Results/volcano_", contrast, "_adjpval.png"),
         plot = volcano.plot, 
         width = 8, 
         height = 8, 
         units = "in")
}

```

## Noemi Requested Comparisons

### Setup

```{r}


library(limma)
library(edgeR)
library(ggrepel)

# Setup annotation groups for the contrast
#tumor.types <- c("Biopsy")
#region.types <- c("immune", "tumor", "vessel")
#source.types <- c("Skin", "GI")
#segment.types <- c("FullROI")

# Gather the column names to be used in defining contrasts
#contrast.groups.list <- list(
#  "Tumor" = tumor.types, 
#  "region" = region.types, 
#  "segment" = segment.types, 
#  "Source" = source.types)

# Filter data for the chosen annotations
object <- q3.normalization.output$object

#subset.indices <- pData(object)$Tumor %in% tumor.types & 
#  pData(object)$region %in% region.types &
#  pData(object)$segment %in% segment.types &
#  pData(object)$Source %in% source.types

#subset.object <- object[, subset.indices]

comparison.column.names <- pData(object) %>% 
  select(starts_with("comparison")) %>% 
  colnames()

print(comparison.column.names)

# Convert all annotations to factors
for(column.name in comparison.column.names){ 
  pData(object)[[column.name]] <- factor(pData(object)[[column.name]])
}

# Factor the slide column
#pData(subset.object)[["slide_name"]] <- 
#  factor(pData(subset.object)[["slide_name"]])

# Create log2 counts
assayDataElement(object = object, elt = "log_q") <-
    assayDataApply(object, 2, FUN = log, base = 2, elt = "q_norm")

```

### Run LM for each comparison and create volcano plot for DEGs

```{r}

# Set up the final results df for adding all results

final.results.df <- data.frame()

# Run for each comparison

for(comparison in comparison.column.names){ 
  
  print(comparison)
  
  subset.indices <- pData(object)[[comparison]] %in% c("A", "B")
  
  subset.object <- object[, subset.indices]

  # Create summary table of group counts
  summary.table.df <- pData(subset.object) %>%  
    select(c(names(comparison.column.names), comparison))
  
  summary.table <- table(summary.table.df)
  
  
  # Gather the normalized log counts (columns are sample IDs)
  norm.log.counts <- subset.object@assayData$log_q
  
  # Gather the annotations (rows are sample IDs)
  annotation.subset <- pData(subset.object) %>% select(all_of(comparison))
  
  # Create the DGE object
  DGE.list <- DGEList(counts = norm.log.counts, 
                    samples = annotation.subset)
  
  
  # Create the LM model design
  design <- model.matrix(formula(paste0("~ 0 + ", comparison)), data = DGE.list$samples)
  
  # Create the fit for the model
  fit <- lmFit(DGE.list$counts, design)
  
  # Set up the contrast
  contrast <- makeContrasts(paste0(comparison, "A - ", comparison, "B"), levels = colnames(coef(fit)))
  
  # Generate the estimate of the contrast
  contrast.estimate <- contrasts.fit(fit, contrast)
  
  # Run Empirical Bayes smoothing of standard errors 
  fit.eb <- eBayes(contrast.estimate, robust = TRUE)
  
  # Generate the results table
  results <- topTable(fit.eb, sort.by = "P", n=Inf)
  
  # Columns to keep in the final results data frame
  result.columns <- c("logFC", "P.Value", "adj.P.Val")
  
  # Reformat the results df
  results <- results %>% 
    select(all_of(result.columns)) %>% 
    rename_with(~ paste0("logFC_", comparison), "logFC") %>% 
    rename_with(~ paste0("pval_", comparison), "P.Value") %>% 
    rename_with(~ paste0("adj.pval_", comparison), "adj.P.Val")
  
  # Create a gene column
  results$gene <- rownames(results)
  
  # Add the results to the master data frame
  if(nrow(final.results.df) > 0){
    # Add to the final data frame
    final.results.df <- merge(final.results.df, results, by = "gene", all = TRUE)
    
  } else{
    # Start the data frame if this is the first comparison
    final.results.df <- results
  }
  
  ## Make a volcano plot for the comparison
  
  # Define the columns for the volcano plot data
  logfc.column.name <- paste0("logFC_", comparison)
  padj.column.name <- paste0("adj.pval_", comparison)
  
  results$logfc <- results[[logfc.column.name]]
  results$padj <- results[[padj.column.name]]
  
  # Create a column for direction of DEGs
  results$de_direction <- "NONE"
  results$de_direction[results$padj < 0.1 & 
                               results$logfc > 0.9] <- "UP"
  results$de_direction[results$padj < 0.1 & 
                               results$logfc < -0.9] <- "DOWN"
  
  # Create a label for DEGs
  results$deglabel <- ifelse(results$de_direction == "NONE", NA, results$gene)
  
  # Compute the scale for the volcano x-axis
  log2.scale <- max(abs(results$logfc))
  
  # Establish the color scheme for the volcano plot
  contrast.level.colors <- c("steelblue4", "grey", "violetred4")
  names(contrast.level.colors) <- c("DOWN", "NONE", "UP")
  
  # Make the volcano plot
  volcano.plot <- ggplot(data = results, aes(x = logfc, 
                                                    y = -log10(padj), 
                                                    col = de_direction, 
                                                    label = deglabel)) +
    geom_vline(xintercept = c(-1, 1), col = "gray", linetype = 'dashed') +
    geom_hline(yintercept = -log10(0.05), col = "gray", linetype = 'dashed') + 
    xlim(-7.5, 7.5) + 
    labs(x = paste("LANA+ versus Full ROI - log FC"),
         y = "-log10 adjusted p-value", 
         title = " DE Analysis", 
         subtitle = comparison) + 
    geom_point(size = 2) +
    scale_color_manual("DE Direction in LANA+", 
                       values = contrast.level.colors) + 
    geom_text_repel(max.overlaps = Inf) + 
    xlim(-log2.scale-1, log2.scale+1) + 
    theme(plot.title = element_text(hjust = 0.5))
  
  # Export the volcano plot
  export.volcano <- TRUE
  if(export.volcano == TRUE){
      ggsave(paste0("volcano_", comparison, ".png"),
           plot = volcano.plot, 
           width = 8, 
           height = 8, 
           units = "in")
  }
  
}


# Export the final data frame
export.final.df <- TRUE
if(export.final.df == TRUE){
  
  write.csv(final.results.df, 
            "CPTR5_Krug_DE_results_noemi_comparisons.csv", 
            row.names = FALSE)
  
}

```

```{r}

library(dplyr)
library(stringr)

# Remove the bad characaters from cleaned annotation treatment column

cleaned.annotation.df <- cleaned.annotation.df %>%
  mutate(treatment = str_replace_all(treatment, "-", ""))

# Select the log counts file for the normalization type of interest

norm.type <- "neg"

if(norm.type == "hk"){
  
  norm.log.counts <- hk_log_counts
  
} else if(norm.type == "neg"){ 
  
  norm.log.counts <- neg_log_counts
  norm.count <- neg.norm_counts
  
} else if(norm.type == "q3"){
  
  norm.log.counts <- q3_log_counts
  
}

# Set up the DGEList object for age groups
norm.log.counts <- norm.log.counts %>% select(-feature)

colnames(norm.log.counts) <- gsub("\\.dcc$", "", colnames(norm.log.counts))

annotation.old <- cleaned.annotation.df[cleaned.annotation.df$age == "old", ]
annotation.young <- cleaned.annotation.df[cleaned.annotation.df$age == "young", ]

# Set up the count tables for old versus young

old.sample.IDs <- annotation.old$Sample_ID
young.sample.IDs <- annotation.young$Sample_ID

old.counts <- norm.log.counts[, old.sample.IDs]
young.counts <- norm.log.counts[, young.sample.IDs]

# Create a DGE list object for each age group

DGE.list.old <- DGEList(counts = old.counts, 
                    samples = annotation.old)

DGE.list.young <- DGEList(counts = young.counts, 
                    samples = annotation.young)

# Set up the design matrix for each age group

design.old <- model.matrix(~ 0 + treatment, data = DGE.list.old$samples)
design.young <- model.matrix(~ 0 + treatment, data = DGE.list.young$samples)

# Run the linear model for each age group

fit.old <- lmFit(DGE.list.old$counts, design.old)
fit.young <- lmFit(DGE.list.young$counts, design.young)

# Set up the contrasts of interest for each treatment type

#old
blz.old.contrast <- makeContrasts(treatmentBLZ945 - treatmentcontrol, 
                                  levels = colnames(coef(fit.old)))

ip549.old.contrast <- makeContrasts(treatmentIP549 - treatmentcontrol, 
                                  levels = colnames(coef(fit.old)))

#young
blz.young.contrast <- makeContrasts(treatmentBLZ945 - treatmentcontrol, 
                                  levels = colnames(coef(fit.young)))

ip549.young.contrast <- makeContrasts(treatmentIP549 - treatmentcontrol, 
                                  levels = colnames(coef(fit.young)))


# Estimate contrast

#old
blz.old.contrast.estimate <- contrasts.fit(fit.old, blz.old.contrast)
ip549.old.contrast.estimate <- contrasts.fit(fit.old, ip549.old.contrast)

#young
blz.young.contrast.estimate <- contrasts.fit(fit.young, blz.young.contrast)
ip549.young.contrast.estimate <- contrasts.fit(fit.young, ip549.young.contrast)

# Run Empirical Bayes smoothing of standard errors

#old
fit.blz.old.eb <- eBayes(blz.old.contrast.estimate, robust = TRUE)
fit.ip549.old.eb <- eBayes(ip549.old.contrast.estimate, robust = TRUE)

#young
fit.blz.young.eb <- eBayes(blz.young.contrast.estimate, robust = TRUE)
fit.ip549.young.eb <- eBayes(ip549.young.contrast.estimate, robust = TRUE)

# Results for each contrast

blz.old.results <- topTable(fit.blz.old.eb, sort.by = "P", n=Inf)
ip549.old.results <- topTable(fit.ip549.old.eb, sort.by = "P", n=Inf)
blz.young.results <- topTable(fit.blz.young.eb, sort.by = "P", n=Inf)
ip549.young.results <- topTable(fit.ip549.young.eb, sort.by = "P", n=Inf)

# Create deg lists for all contrasts
all.deg.lists <- list(blz.old.degs = blz.old.results, 
                  ip.old.degs = ip549.old.results, 
                  blz.young.degs = blz.young.results, 
                  ip.young.degs = ip549.young.results)

# Filter for adj p-value and add column for significance level
for(deg.list in names(all.deg.lists)){ 
  all.deg.lists[[deg.list]] <- all.deg.lists[[deg.list]]  %>% 
    filter(adj.P.Val < 0.05) %>% 
    mutate(sig.level = case_when(
      adj.P.Val > 0.01 & adj.P.Val < 0.05 ~ "*",
      adj.P.Val > 0.001 & adj.P.Val <= 0.01 ~ "**",
      adj.P.Val <= 0.001 ~ "***"))
  }


# Write all four DE results to csv files

export.de <- FALSE

if(export.de == TRUE){
  
  write.csv(blz.old.results, file.path(paste0(results.folder, "DE/", norm.type,    "_DE_Old_Blz.csv")), 
                row.names = TRUE)
  
  write.csv(ip549.old.results, file.path(paste0(results.folder, "DE/", norm.type,    "_DE_Old_IP549.csv")), 
                row.names = TRUE)
  
  write.csv(blz.young.results, file.path(paste0(results.folder, "DE/", norm.type,    "_DE_Young_Blz.csv")), 
                row.names = TRUE)
  
  
  write.csv(ip549.young.results, file.path(paste0(results.folder, "DE/", norm.type,    "_DE_Young_IP549.csv")), 
                row.names = TRUE)
  
}

```