CPTR5_noemi_contrasts_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

## Noemi Requested Comparisons

### Setup

```{r}


library(limma)
library(edgeR)
library(ggrepel)
library(GeomxTools)
library(dplyr)

# 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)
  
}

```