analysis_markdown.Rmd

---
title: "nCounter_RNA_Analysis"
author: "Ned Cauley"
date: "4/16/24"
output: html_document
---

# This markdown and analysis based on https://github.com/bhattacharya-a-bt/CBCS_normalization/blob/master/sabry_analysis.R
# The above analysis is from the publication https://academic.oup.com/bib/article/22/3/bbaa163/5891144?login=true

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

## Set up libraries and QC functions libraries

```{r Libraries and RUVseq functions, message=FALSE, warning=FALSE}
library(NanoStringQCPro)
library(ggplot2)
library(EnvStats)
library(MASS)
library(RUVSeq)
library(DESeq2)
library(limma)
library(matrixStats)
library(limma)
library(NanoStringNorm)
library(readxl)
library(plotly)
library(pheatmap)
library(ggrepel)
library(gridExtra)
library(grid)


# The library in RSW to use is located at
# /home/cauleyes/ccbr-data/renv_cache/ncounter_rna/Snapshot-environment_method/renv/library/R-4.1/x86_64-pc-linux-gnu

# Set the working directory to your project folder if it is not set already
setwd("/home/cauleyes/ccbr-data/users/Ned/nCounter_rna/CCBR1245_Ramaswami_nCounter")

# Functions for RUVseq analysis
#source("R/nanostring_RUV_functions.R")
source("R/ncounter_analysis_functions.R")

# Load the annotation file
#annotation_file_path <- "ncounter_rna_annotation_edit_correction_RR.xlsx"
annotation_file_path <- "annotations/5_20_24_ncounter_rna_annotation_edit_NC.csv"
#annotation_df <- read_excel(annotation_file_path)
annotation_df <- read.csv(annotation_file_path)

```

```{r Define the contrast for the DEG analysis, message=FALSE, warning=FALSE}

# Define the contrast and levels from the annotation columns for QC and DEG analysis
# Levels are: Effect, Reference. 
# This means a positive log2foldchnage = lower reference, higher effect
contrasts.list <- list(
  "CD4_count_200" = c("High", "Low"), 
  "CD4_count_175" = c("High", "Low"), 
  "CD4_count_150" = c("High", "Low"), 
  "KS_alone" = c("N", "Y"), 
  "HIV_VL" = c("Y", "N"), 
  "response" = c("Y", "N"), 
  "kshv_expression_median" = c("above", "below")
)



# For running a single contrast
#contrasts.list <- list(
  "kshv_expression_median" = c("above", "below")
#)

#contrast <- "CD4_group"
#contrast <- "KS_alone"
#contrast <- "KSHV_VL_Bin"


# Levels are: Effect, Reference. 
# This means a positive log2foldchnage = lower reference, higher effect
#contrast_levels <- c("High", "Low")
#contrast_levels <- c("N", "Y")
#contrast_levels <- c("0", "1")





```


## Setup the initial dataframes with read counts and metadata/QC

```{r Initial Dataframes and QC,  message=FALSE, warning=FALSE}

rcc.files.path <- "/rstudio-files/ccbr-data/users/Ned/nCounter_rna/data/CCBR1245/rcc"

# Gather the names of the RCC files
#files.RCC = list.files("/rstudio-files/ccbr-data/users/Ned/nCounter_rna/data/CCBR1245/rcc")

# Process the RCC files as count, metadata, and feature data DFs
rcc.dfs <- process_RCC_files(rcc.files.path = rcc.files.path, 
                  annotation.df = annotation_df)

# Print a sumnmary of the DFs
head(rcc.dfs$raw.counts)
head(rcc.dfs$pData)
head(rcc.dfs$fData)

```

## Verify that housekeeping genes are not DEGs

```{r, message=FALSE, warning=FALSE}
###
#### CHECK IF HK ARE ASSOCIATED WITH PRIMARY PHENO
###

# Gather the read counts for housekeeping genes
hk_raw = rcc.dfs$raw.counts[rcc.dfs$hk.gene.list,]
pval = vector(length = nrow(hk_raw))

print("Housekeeping genes:")
print(rownames(hk_raw))

# Run general linear model to get pvalue for housekeeping
for (contrast in names(contrasts.list)) {
  
  for (i in 1:nrow(hk_raw)){
  
  reg = glm.nb(as.numeric(hk_raw[i,]) ~ as.factor(rcc.dfs$pData[[contrast]]))
  pval[i] = coef(summary(reg))[2,4]
  

  
  print("Number of housekeeping genes that are DE:")
  sum(pval <= .05)
  
}
  
  
}


#print("Housekeeping genes:")
#print(rownames(hk_raw))

# Number of housekeeping that are DE
#print("Number of housekeeping genes that are DE:")
#sum(pval <= .05)

```

## Run RUVseq differential expression analysis

```{r, message=FALSE, warning=FALSE}

# Create a new folder for the run date and results
run.folder <- "7_17_2024/"

results.folder <- "/rstudio-files/ccbr-data/users/Ned/nCounter_rna/CCBR1245_data_results/Results/"

for(contrast in names(contrasts.list)){
  
  contrast.levels <- contrasts.list[[contrast]]
  
  # Run DESEQ with RUVseq normalization for k sources of variability
  deseq.output <- run_RUV_DESEQ2(contrast = contrast, 
                 contrast.levels = contrast.levels, 
                 annotation.data = rcc.dfs$pData, 
                 raw.counts = rcc.dfs$raw.counts, 
                 feature.data = rcc.dfs$fData, 
                 k = 1, 
                 additional.covar = NULL, 
                 hk.genes = rcc.dfs$hk.gene.list, 
                 exclude.hk.genes = NULL)
  
  
  ## Make a volcano plot for the DEGs
  
  # Identify the NEG and POS probes
  NEG.indices <- grep("NEG_", rownames(deseq.output$all.de.results))
  POS.indices <- grep("POS_", rownames(deseq.output$all.de.results))
  control.indices <- c(NEG.indices, POS.indices)
  
  # Remove the control probes
  deg.results.df <- deseq.output$all.de.results[-control.indices,]
  
  # Create a column for direction of DEGs
  deg.results.df$de_direction <- "NONE"
  deg.results.df$de_direction[deg.results.df$padj < 0.05 & 
                               deg.results.df$log2FoldChange > 1] <- "UP"
  deg.results.df$de_direction[deg.results.df$padj < 0.05 & 
                               deg.results.df$log2FoldChange < -1] <- "DOWN"
  
  deg.results.df$deglabel <- ifelse(deg.results.df$de_direction == "NONE", NA, deg.results.df$gene)
  
  # Compute the scale for the volcano x-axis
  log2.scale <- max(abs(deg.results.df$log2FoldChange))
  
  # 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 = deg.results.df, aes(x = log2FoldChange, 
                                                    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(contrast, contrast.levels[[1]], "versus", contrast.levels[[2]]),
         y = "-log10 adjusted p-value", 
         title = paste0(contrast, " DE Analysis")) + 
    geom_point(size = 2) +
    scale_color_manual("DE Direction", 
                       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(results.folder, run.folder, "volcano_", contrast, "_adjpval.png"),
           plot = volcano.plot, 
           width = 8, 
           height = 8, 
           units = "in")
  }
  
  # Export Normalized read counts
  export_norm_counts <- FALSE
  if(export_norm_counts == TRUE){
    # Grab the normalized_counts
    normalized_counts <- as.data.frame(counts(deseq.output$dds, normalized = TRUE))
    
    # Add a gene column
    gene <- rownames(normalized_counts)
    normalized_counts <- cbind(gene, normalized_counts) 
    
    # Export
    write.csv(normalized_counts, file = paste0(results.folder, run.folder, "RUV_normalized_counts.csv"), row.names = FALSE)
    
    
  }
  
  # Export the DEG list
  export_sig_DEG_list <- TRUE
  if(export_sig_DEG_list == TRUE){
    write.csv(deseq.output$deg.list, file = paste0(results.folder, run.folder, "DEG_list_",contrast,".csv"), row.names =   FALSE)
  }
  
  # Export all DE results
  export_DE_results <- FALSE
  if(export_DE_results == TRUE){
    write.csv(deg.results.df, file = paste0(results.folder, run.folder, "DE_results_",contrast,".csv"), row.names =   FALSE)
  }
  
  
  # Create a heatmap for the DEGs from the DEG contrast

  # Get a list of the DEGS
  rownames(deseq.output$deg.list) <- deseq.output$deg.list$gene
  contrast_deg_list <-   rownames(deseq.output$deg.list[deseq.output$deg.list$padj < 0.05,   ])
  
  subtitle <- "padj < 0.05"
  
  # Loosen the cutoff if only 1 DEG found (causes an error in pheatmap)
  if(length(contrast_deg_list) < 1){ 
    
    print("No DEGs found for padj < 0.05, trying pvalue < 0.05")
    
    contrast_deg_list <-   rownames(deseq.output$deg.list[deseq.output$deg.list$pvalue < 0.05,
                  ])
    subtitle <- "p-value < 0.05"
    
  }
  
  # Remove the positive and negative probes
  # Identify the NEG and POS probes
  NEG.indices <- grep("NEG_", contrast_deg_list)
  POS.indices <- grep("POS_", contrast_deg_list)
  control.indices <- c(NEG.indices, POS.indices)
  
  # Remove the control probes
  if(length(control.indices) > 0){
    goi <- contrast_deg_list[-control.indices]
  } else{
    goi <- contrast_deg_list
  }
  
  
  # Gather the annotation for the contrast
  anno.columns <- c(contrast)
  annotation_df <- as.data.frame(annotation_df)
  
  # Ensure the row names are the sample names
  rownames(annotation_df) <- annotation_df$SampleID
  contrast.annotation <- annotation_df[ ,anno.columns, drop = FALSE]
  
  # Subset heatmap annotation for any samples taken out in counts
  # Gather the samples in both
  anno.samples.rownames <- rownames(contrast.annotation)
  counts.samples.colnames <- colnames(deseq.output$dds)

  # Find the samples that are present in both
  matching_names <- intersect(anno.samples.rownames, counts.samples.colnames)

  # Subset counts to include only overlapping samples
  heatmap.annotation <- contrast.annotation[matching_names, , drop = FALSE]
  
  # Order by contrast level
  heatmap.annotation <- heatmap.annotation[order(heatmap.annotation[1]), , drop = FALSE]
  
  # Reorder the data columns to they match the annotation row order
  row.order <- rownames(heatmap.annotation)
  
  log.norm.counts <- log2(counts(deseq.output$dds, normalized = TRUE) + 1)
  log.norm.counts <- log.norm.counts[, row.order]
  
  # Define the annotation colors for the contrast
  anno.colors <- setNames(list(c("cyan", "magenta")), contrast)
  
  # Assign the levels to the colors
  anno.colors[[contrast]] <- setNames(anno.colors[[contrast]], contrast.levels)

  
  heatmap.contrast <- pheatmap(log.norm.counts[goi, ], 
           main = paste0(contrast, "DEGs (", subtitle, ")"), 
           show_rownames = TRUE, 
           scale = "row", 
           show_colnames = TRUE,
           border_color = NA, 
           cluster_rows = TRUE, 
           cluster_cols = FALSE, 
           clustering_method = "average", 
           clustering_distance_rows = "correlation", 
           clustering_distance_cols = "correlation", 
           color = colorRampPalette(c("blue", "white", "red"))(120),   
           annotation_row = NA, 
           annotation_col = heatmap.annotation, 
           annotation_colors = anno.colors
           )
  
  
  export.heatmap.contrast <- TRUE
  if(export.heatmap.contrast == TRUE){
    
    ggsave(paste0(results.folder, run.folder, "heatmap_", contrast,   ".png"),
           plot = heatmap.contrast, 
           width = 12, 
           height = 12, 
           units = "in")
    
  }
 
}



```

# Unsupervised analysis

```{r}

# Plot PCA to check for batch effects

# Prepare the counts
rlog_dds <- rlog(deseq.output$dds)

# Annotation to view PCA by
PCA.group <- "kshv_expression_median"

# PCA Plot
PCA.plot <- plotPCA(rlog_dds, intgroup = PCA.group, ntop = 700)

# Plotly PCA for interactivity
PCA.plot.data <- plotPCA(rlog_dds, intgroup = PCA.group, ntop = 700, returnData = TRUE)

PCA.plotly <- plot_ly(data = PCA.plot.data, x = ~PC1, y = ~PC2, text = ~name, mode = "markers", color = ~factor(group)) %>%
  layout(title = "PCA Plot",
         xaxis = list(title = "PC1"),
         yaxis = list(title = "PC2"))

print(PCA.plotly)

# Export the PCA plot

export.PCA <- FALSE
if(export.PCA == TRUE){
  ggsave(paste0("Results/QC/PCA.pdf"), plot = PCA.plot)
}

```
# Get the read count information for CXCL8

```{r}

# Grab the normalized counts
CXCL8_counts <- normalized_counts["CXCL8", ]

# Make the rows as SampleIDs
CXCL8_df <- as.data.frame(t(CXCL8_counts))
CXCL8_df$SampleID <- rownames(CXCL8_df)

# Set up the annotations
CXCL8_anno_columns <- c("SampleID","KS_alone","CD4_count","response")
subset_anno_df <- annotation_df[CXCL8_anno_columns]

# Correct SampleID names in annotation df
subset_anno_df <- subset_anno_df %>% mutate(SampleID = str_replace_all(SampleID, "-", " "))

# Create a master df
CXCL8_df <- merge(CXCL8_df, subset_anno_df, by = "SampleID")

CXCL8_mean_KS_Y <- mean(CXCL8_df$CXCL8[CXCL8_df$KS_alone == "Y"], na.rm = TRUE)
CXCL8_mean_KS_N <- mean(CXCL8_df$CXCL8[CXCL8_df$KS_alone == "N"], na.rm = TRUE)

CXCL8_df$KS_alone_mean <- ifelse(CXCL8_df$KS_alone == "Y", CXCL8_mean_KS_Y, CXCL8_mean_KS_N)

CXCL8_mean_CD4_High <- mean(CXCL8_df$CXCL8[CXCL8_df$CD4_count == "High"], na.rm = TRUE)
CXCL8_mean_CD4_Low <- mean(CXCL8_df$CXCL8[CXCL8_df$CD4_count == "Low"], na.rm = TRUE)

CXCL8_df$CD4_count_mean <- ifelse(CXCL8_df$CD4_count == "High", CXCL8_mean_CD4_High, CXCL8_mean_CD4_Low)

CXCL8_mean_response_Y <- mean(CXCL8_df$CXCL8[CXCL8_df$response == "Y"], na.rm = TRUE)
CXCL8_mean_response_N <- mean(CXCL8_df$CXCL8[CXCL8_df$response == "N"], na.rm = TRUE)

CXCL8_df$response_mean <- ifelse(CXCL8_df$response == "Y", CXCL8_mean_response_Y, CXCL8_mean_response_N)



#CXCL8_deg <- ruv_deg_list_export[ruv_deg_list_export$gene == "CXCL8", ]

export_CXCL8_summary <- FALSE
if(export_CXCL8_summary == TRUE){
  write.csv(CXCL8_df, file = "Results/CXCL8_summary.csv", row.names = FALSE)
}


```

# Heatmap of Variable Genes

```{r}

# Prepare the normalized counts
log.norm.counts <- log2(counts(dds, normalized = TRUE) + 1)

# Prepare the annotation columns
anno.columns <- c("KS_alone", "CD4_group")
annotation_df <- as.data.frame(annotation_df)

# Ensure the row names are the sample names
rownames(annotation_df) <- annotation_df$SampleID
heatmap.annotation <- annotation_df[ ,anno.columns]

# Order the annotation
heatmap.annotation <- heatmap.annotation %>% 
  arrange(CD4_group) %>% arrange(KS_alone)

# Define the colors for the annotations
anno.colors = list(
  KS_alone = c(Y = "seagreen", N = "red4"), 
  CD4_group = c(High = "purple", Low = "coral")
  )

# Determine the most variable genes using the coefficient of variation
.calcCv <- function(x) {
  sd(x) / mean(x)
}
cv.vector <- apply(log.norm.counts, 1, .calcCv)

# Identify the NEG and POS probes
NEG.indices <- grep("NEG_", names(cv.vector))
POS.indices <- grep("POS_", names(cv.vector))
control.indices <- c(NEG.indices, POS.indices)

# Print the names of the control probes to be removed
control.probes <- names(cv.vector)[control.indices]
print(control.probes)

# Remove the control probes
cv.vector <- cv.vector[-control.indices]

# Sort the genes based on variability
cv.vector <- sort(cv.vector, decreasing = TRUE)

# Top variable genes
ngenes <- 50
goi <- names(cv.vector)[1:ngenes]

# Reorder the data columns to they match the annotation row order
row.order <- rownames(heatmap.annotation)
log.norm.counts <- log.norm.counts[, row.order]


# Top 50 genes
heatmap.top50 <- pheatmap(log.norm.counts[goi, ], 
         main = "Top 50 Variable Genes", 
         show_rownames = TRUE, 
         scale = "row", 
         show_colnames = TRUE,
         border_color = NA, 
         cluster_rows = TRUE, 
         cluster_cols = FALSE, 
         clustering_method = "average", 
         clustering_distance_rows = "correlation", 
         clustering_distance_cols = "correlation", 
         color = colorRampPalette(c("blue", "white", "red"))(120), 
         annotation_row = NA, 
         annotation_col = heatmap.annotation, 
         annotation_colors = anno.colors
         )

export.top50 <- FALSE
if(export.top50 == TRUE){
  
  ggsave("Results/heatmap_top50.png",
         plot = heatmap.top50, 
         width = 10, 
         height = 10, 
         units = "in")
  
}

# Create a heatmap with all genes
heatmap.all <- pheatmap(log.norm.counts, 
         main = "All Genes", 
         show_rownames = FALSE, 
         scale = "row", 
         show_colnames = TRUE,
         border_color = NA, 
         cluster_rows = TRUE, 
         cluster_cols = FALSE, 
         clustering_method = "average", 
         clustering_distance_rows = "correlation", 
         clustering_distance_cols = "correlation", 
         color = colorRampPalette(c("blue", "white", "red"))(120), 
         annotation_row = NA, 
         annotation_col = heatmap.annotation, 
         annotation_colors = anno.colors
         )

export.heatmap.all <- FALSE
if(export.heatmap.all == TRUE){
  
  ggsave("Results/heatmap_all.png",
         plot = heatmap.all, 
         width = 10, 
         height = 10, 
         units = "in")
  
}


# Create a heatmap for the DEGs from the DEG contrast

# Get a list of the DEGS
contrast_deg_list <- rownames(contrast_ruv_deg_list[contrast_ruv_deg_list$padj < 0.05, ])

# Loosen the cutoff if no DEGs found
if(length(rownames(contrast_deg_list)) == 0){ 
  
  print("No DEGs found for padj < 0.05, trying padj < 0.1")
  
  contrast_deg_list <- rownames(contrast_ruv_deg_list[contrast_ruv_deg_list$padj < 0.1, ])
}

if(length(rownames(contrast_deg_list)) == 0){ 
  
  print("No DEGs found for padj < 0.05, trying pvalue < 0.05")
  
  contrast_deg_list <- rownames(contrast_ruv_deg_list[contrast_ruv_deg_list$pvalue < 0.05, ])
  
}

# Remove the positive and negative probes
# Identify the NEG and POS probes
NEG.indices <- grep("NEG_", contrast_deg_list)
POS.indices <- grep("POS_", contrast_deg_list)
control.indices <- c(NEG.indices, POS.indices)

# Remove the control probes
if(length(control.indices) > 0){
  goi <- contrast_deg_list[-control.indices]
} else{
  goi <- contrast_deg_list
}


# Gather the annotation for the contrast
anno.columns <- c(contrast)
annotation_df <- as.data.frame(annotation_df)

# Ensure the row names are the sample names
rownames(annotation_df) <- annotation_df$SampleID
heatmap.annotation <- annotation_df[ ,anno.columns, drop = FALSE]

# Order the annotation
heatmap.annotation <- heatmap.annotation %>% arrange(.data[[contrast]])

# Reorder the data columns to they match the annotation row order
row.order <- rownames(heatmap.annotation)
log.norm.counts <- log.norm.counts[, row.order]

# Define the annotation colors for the contrast
anno.colors1 = list(
  KS_alone = c(Y = "seagreen", N = "red4")
  )

anno.colors2 = list(
  CD4_group = c(High = "purple", Low = "coral")
  )

heatmap.contrast <- pheatmap(log.norm.counts[goi, ], 
         main = paste0("DEGs for ", contrast), 
         show_rownames = TRUE, 
         scale = "row", 
         show_colnames = TRUE,
         border_color = NA, 
         cluster_rows = TRUE, 
         cluster_cols = FALSE, 
         clustering_method = "average", 
         clustering_distance_rows = "correlation", 
         clustering_distance_cols = "correlation", 
         color = colorRampPalette(c("blue", "white", "red"))(120), 
         annotation_row = NA, 
         annotation_col = heatmap.annotation, 
         annotation_colors = anno.colors2
         )


export.heatmap.contrast <- TRUE
if(export.heatmap.contrast == TRUE){
  
  ggsave(paste0(results.folder, run.folder, "heatmap_", contrast, ".png"),
         plot = heatmap.contrast, 
         width = 12, 
         height = 12, 
         units = "in")
  
}

print(heatmap.top50)
print(heatmap.all)
print(heatmap.contrast)



```


# Clustering of Viral Gene Probes

```{r}

# See reference vignette: https://bioconductor.org/packages/release/bioc/vignettes/PCAtools/inst/doc/PCAtools.html#introduction

library(stringr)
library(dplyr)
library(ggfortify)
library(PCAtools)

# Transform DESEQ normalized counts to variance stabilized
vst.dds <- varianceStabilizingTransformation(deseq.output$dds)

# Gather normalized counts as a df
normalized.counts <- as.data.frame(assay(vst.dds))

# Identify the KSHV genes
kshv.genes <- str_detect(rownames(normalized.counts), "^KSHV")

# Subset for KSHV genes
normalized.counts.kshv <- normalized.counts[kshv.genes, ]
normalized.counts.notkshv <- normalized.counts[-kshv.genes, ]

# Transpose gene counts
normalized.counts.kshv.transpose <- as.data.frame(t(normalized.counts.kshv))

export.kshv.gene.counts <- FALSE
if(export.kshv.gene.counts == TRUE){
  
  write.csv(normalized.counts.kshv, file = "Results/kshv_counts_samples_as_columns.csv")
  write.csv(normalized.counts.kshv.transpose, file = "Results/kshv_counts_samples_as_rows.csv")
  
  
# Gather the average expression per sample
kshv.means <- colMeans(normalized.counts.kshv)
notkshv.means <- colMeans(normalized.counts.notkshv)

total.median.kshv <- median(kshv.means)
total.mean.kshv <- mean(kshv.means)

total.median.notkshv <- median(notkshv.means)
total.mean.notkshv <- mean(notkshv.means)

label.data <- data.frame(
  x = c(total.median.kshv, 
        total.mean.kshv, 
        total.median.notkshv, 
        total.mean.notkshv),
  y = c(0.05, 0.07, 0.59, 0.61),
  label = c(paste("KSHV Median:", round(total.median.kshv, 2)), 
            paste("KSHV Mean:", round(total.mean.kshv, 2)), 
            paste("Other Genes Median:", round(total.median.notkshv, 2)), 
            paste("Other Genes Mean:", round(total.mean.notkshv, 2)))
)

# Make a histogram of KSHV gene counts
kshv.distribution.plot <- ggplot(data.frame(kshv.means), aes(x = kshv.means)) +
  geom_histogram(aes(y = ..density..), bins = 14, fill = "lightblue", color = "black", alpha = 0.5) + 
  geom_density(alpha=0.6) +
  geom_vline(xintercept = total.median.kshv, color = "coral", linetype = "dashed", size = 1) + 
  geom_vline(xintercept = total.mean.kshv, color = "maroon", linetype = "dashed", size = 1) +
  geom_vline(xintercept = total.median.notkshv, color = "lightgreen", linetype = "dashed", size = 1) +
  geom_vline(xintercept = total.mean.notkshv, color = "darkgreen", linetype = "dashed", size = 1) + 
  labs(title = "Distribution of Average KSHV Gene Counts per Sample",
       x = "Average Gene Count",
       y = "Frequency") + 
  geom_text_repel(data = label.data, aes(x = x, y = y, label = label), color = "black", vjust = -0.5) +
  theme_minimal()
  
  
}
```

```{r}
# Run PCA
kshv.pca <- pca(normalized.counts.kshv, 
                metadata = merged_pData, 
                removeVar = 0.1)

# Plot PCA

kshv.pca.plot.by.sample <- biplot(kshv.pca, 
       colby = contrast, 
       legendPosition = 'top', 
       legendLabSize = 16, 
       legendIconSize = 8.0)

screeplot <- screeplot(kshv.pca)

# Visualize clustering of genes

kshv_annotation <- as.data.frame(read_excel("annotations/kshv_genes_annotation.xlsx"))

rownames(kshv_annotation) <- kshv_annotation$gene

# Ensure the rows and columns of the metadata and counts are aligned
kshv_annotation <- kshv_annotation[order(rownames(kshv_annotation)), ]
normalized.counts.kshv.transpose <- normalized.counts.kshv.transpose[ order(colnames(normalized.counts.kshv.transpose))]

# Create PCA for genes 

kshv.pca.by.gene <- pca(normalized.counts.kshv.transpose, 
                        metadata = kshv_annotation, 
                        removeVar = 0.1)


kshv.pca.plot.by.gene <- biplot(kshv.pca.by.gene, 
       colby = "type", 
       legendPosition = 'top', 
       legendLabSize = 16, 
       legendIconSize = 8.0)

export.PCA <- FALSE
if(export.PCA == TRUE){
  
  ggsave(paste0("Results/QC/kshv_pca_plot_by_sample_", contrast, ".png"),
         plot = kshv.pca.plot.by.sample, 
         width = 8, 
         height = 8, 
         units = "in")
  
  ggsave(paste0("Results/QC/kshv_pca_plot_by_gene.png"),
         plot = kshv.pca.plot.by.gene, 
         width = 8, 
         height = 8, 
         units = "in")
  
}

```

```{r Divide samples by KSHV average expression}

# Determine the averge per sample
average.kshv.counts <- colMeans(normalized.counts.kshv)

# Set up a data frame for determining groups
avg.kshv.df <- data.frame(SampleID = names(average.kshv.counts), 
                          meanKSHV = unname(average.kshv.counts))

# Find the median of all samples' average KSHV gene expression
median.kshv <- median(avg.kshv.df$meanKSHV)

# Assign samples to high or low groups based on the median KSHV average
avg.kshv.df <- avg.kshv.df %>% 
  mutate(kshv_expression_median = ifelse(meanKSHV > median.kshv, "above", "below"))




# Merge with the annotation
annotation_df_cleaned <- annotation_df %>% 
  mutate(SampleID = str_replace_all(SampleID, "-", " "))

annotation.update <- merge(annotation_df_cleaned, avg.kshv.df, by = "SampleID")

# Export a new annotation with kshv group assignments
export.annotation <- FALSE
if(export.annotation == TRUE){
  
  write.csv(annotation.update,
             paste0("annotations/5_20_24_ncounter_rna_annotation_edit_NC.csv"), 
            row.names = FALSE)
  
}


```

# Boxplots of Viral Gene Probes

```{r}

library(tidyr)

# Boxplot of Lytic versus Latent for individual samples

# Set up a column for the viral gene in the counts df
normalized.counts.kshv$gene <- rownames((normalized.counts.kshv))

# Merge into a master DF
kshv.df <- merge(normalized.counts.kshv, kshv_annotation, by = "gene")


# Melt the df to a long version
kshv.df.long <- kshv.df %>%
    pivot_longer(cols = starts_with("SB"), names_to = "sample", values_to = "count")

# Gather the sample names
sample.names <- colnames(kshv.df)[grep("SB", colnames(kshv.df))]

for(ID in sample.names){
  
  # Filter for a specific sample
  kshv.df.sample <- kshv.df.long %>% filter(sample == ID)
  
  # Make a boxplot of lytic versus latent
  kshv.boxplot.type <- ggplot(kshv.df.sample, aes(x = type, y = count)) +
    geom_boxplot() +
    labs(title = paste0("Lytic vs Latent KSHV genes for ", ID),
         x = NULL,
         y = "Count")
  
  kshv.boxplot.all <- ggplot(kshv.df.sample, aes(x = sample, y = count)) +
    geom_boxplot() +
    labs(title = paste0("All KSHV gene counts for ", ID),
         x = NULL,
         y = "Count")
  
  # Combine into a single image
  ID.combined.boxplot <- grid.arrange(kshv.boxplot.all, kshv.boxplot.type, ncol = 2)
  
  # Save the plot
  ggsave(paste0(results.folder, "kshv_genes/", run.folder, ID, "_kshv_genes_boxplot.png"), 
        plot = ID.combined.boxplot, 
        width = 8, 
        height = 8, 
        units = "in")
}

kshv.boxplot.all.samples <- ggplot(kshv.df.long, aes(x = sample, y = count)) +
  geom_boxplot() +
  labs(title = "Lytic vs Latent Gene Counts by Sample",
         x = "Gene Type",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(paste0(results.folder, "kshv_genes/", run.folder, "kshv_all_boxplot.png"), 
        plot = kshv.boxplot.all.samples, 
        width = 14, 
        height = 8, 
        units = "in")

# Combine sample and type
kshv.df.long$sample_type <- paste0(kshv.df.long$sample, "_", kshv.df.long$type)

kshv.boxplot.sample.types <- ggplot(kshv.df.long, aes(x = sample, y = count, color = type)) +
  geom_boxplot() +
  labs(title = "Lytic vs Latent Gene Counts by Sample",
         x = "Gene Type",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 5))


ggsave(paste0(results.folder, "kshv_genes/", run.folder, "kshv_type_boxplot.png"), 
        plot = kshv.boxplot.sample.types, 
        width = 14, 
        height = 8, 
        units = "in")

kshv.genes.boxplot <- ggplot(kshv.df.long, aes(x = gene, y = count, color = type)) +
  geom_boxplot() +
  labs(title = "KSHV Gene Counts for all Samples",
         x = "KSHV Gene",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 5))


ggsave(paste0(results.folder, "kshv_genes/", run.folder, "kshv_genes_boxplot.png"), 
        plot = kshv.genes.boxplot, 
        width = 14, 
        height = 8, 
        units = "in")

write.csv(kshv.df.long, file = paste0(results.folder, "kshv_genes/", run.folder, "kshv_gene_df.csv"), row.names = FALSE)

```

# DSP Samples KSHV Gene Analysis

```{r}

dsp.samples <- c("SB 18 4726", "SB 22 2794", "SB 22 2963", "SB 22 696")

kshv.df.long.dsp <- kshv.df.long %>% 
  filter(sample %in% dsp.samples)


kshv.boxplot.all.dsp.samples <- ggplot(kshv.df.long.dsp, aes(x = sample, y = count)) +
  geom_boxplot() +
  labs(title = "Lytic vs Latent Gene Counts for DSP Overlap Samples",
         x = "Gene Type",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


kshv.boxplot.dsp.sample.types <- ggplot(kshv.df.long.dsp, aes(x = sample, y = count, color = type)) +
  geom_boxplot() +
  labs(title = "Lytic vs Latent Gene Counts for DSP Overlap Samples",
         x = "Gene Type",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 5))


kshv.genes.boxplot <- ggplot(kshv.df.long.dsp, aes(x = gene, y = count, color = type)) +
  geom_boxplot() +
  labs(title = "KSHV Gene Counts for all Samples",
         x = "KSHV Gene",
         y = "Count") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 5))

```

# DSP Samples CD4 Analysis

```{r}

CD4.counts <- normalized.counts %>% 
  filter(gene == "CD4") %>% 
  select(all_of(dsp.samples))

```


# DSP Samples DE Analysis

```{r}

# Create a subset DF with only the DSP overlap samples
dsp.samples.df <- rcc.dfs

# Subset the annotation
dsp.samples.df$pData <- dsp.samples.df$pData %>% 
  filter(SampleID %in% dsp.samples)

# Subset the counts
dsp.samples.df$raw.counts <- raw.count.df %>% 
  dplyr::select(all_of(dsp.samples))

  
# Run DESEQ with RUVseq normalization for k sources of variability
deseq.output <- run_RUV_DESEQ2(contrast = "CD4_count_200", 
                               contrast.levels = c("High", "Low"), 
                               annotation.data = dsp.samples.df$pData, 
                               raw.counts = dsp.samples.df$raw.counts, 
                               feature.data = dsp.samples.df$fData, 
                               k = 1, 
                               additional.covar = NULL, 
                               hk.genes = dsp.samples.df$hk.gene.list, 
                               exclude.hk.genes = NULL)
  
  
  ## Make a volcano plot for the DEGs
  
# Identify the NEG and POS probes
NEG.indices <- grep("NEG_", rownames(deseq.output$all.de.results))
POS.indices <- grep("POS_", rownames(deseq.output$all.de.results))
control.indices <- c(NEG.indices, POS.indices)

# Remove the control probes
deg.results.df <- deseq.output$all.de.results[-control.indices,]

# Create a column for direction of DEGs
deg.results.df$de_direction <- "NONE"
deg.results.df$de_direction[deg.results.df$padj < 0.05 & 
                             deg.results.df$log2FoldChange > 1] <- "UP"
deg.results.df$de_direction[deg.results.df$padj < 0.05 & 
                             deg.results.df$log2FoldChange < -1] <- "DOWN"

deg.results.df$deglabel <- ifelse(deg.results.df$de_direction == "NONE", NA, deg.results.df$gene)
  
# Compute the scale for the volcano x-axis
log2.scale <- max(abs(deg.results.df$log2FoldChange))

# 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 = deg.results.df, aes(x = log2FoldChange, 
                                                  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 = "CD4 High versus Low",
       y = "-log10 adjusted p-value", 
       title = "DSP Sample Overlap DE Analysis, CD4 High versus Low") + 
  geom_point(size = 2) +
  scale_color_manual("DE Direction", 
                     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(results.folder, run.folder, "overlap_ncounter_volcano_plot.png"),
           plot = volcano.plot, 
           width = 8, 
           height = 8, 
           units = "in")
  }
  
  # Export Normalized read counts
  export_norm_counts <- FALSE
  if(export_norm_counts == TRUE){
    # Grab the normalized_counts
    normalized_counts <- as.data.frame(counts(deseq.output$dds, normalized = TRUE))
    
    
    # Add a gene column
    gene <- rownames(normalized_counts)
    normalized_counts <- cbind(gene, normalized_counts) 
    
    # Export
    write.csv(normalized_counts, file = paste0(results.folder, run.folder, "RUV_normalized_counts.csv"), row.names = FALSE)
    
  }
  
  # Export the DEG list
  export_sig_DEG_list <- FALSE
  if(export_sig_DEG_list == TRUE){
    write.csv(deseq.output$deg.list, file = paste0(results.folder, run.folder, "DEG_list_",contrast,".csv"), row.names =   FALSE)
  }
  
  # Export all DE results
  export_DE_results <- FALSE
  if(export_DE_results == TRUE){
    write.csv(deg.results.df, file = paste0(results.folder, run.folder, "overlap_ncounter_de_results.csv"), row.names =   FALSE)
  }
  
  


```

# CIBERSORT

```{r}

# Load the reference file
# LM22 reference: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5895181/
lm22.file.path <- "/rstudio-files/ccr-dceg-data/users/ned/CPTR/CPTR-6 Annunziata/cell_decon_reference/LM22.txt"
lm22.reference <- read.csv(lm22.file.path, sep = "\t")

# Compare probe set to CIBERSORT reference
lm22.reference.genes <- lm22.reference$Gene.symbol

save(lm22.reference.genes, file = "/rstudio-files/ccbr-data/users/Ned/nCounter_rna/CCBR1245_data_results/cibersort/reference_genes.rda")


# Load the normalized counts
norm.counts <- read.csv("/rstudio-files/ccbr-data/users/Ned/nCounter_rna/CCBR1245_data_results/Results/5_20_2024/RUV_normalized_counts.csv")

# Determine genes missing from the reference in the ncounter set
ncounter.gene.list <- norm.counts$gene

save(lm22.reference.genes, ncounter.gene.list, file = "/rstudio-files/ccbr-data/users/Ned/nCounter_rna/CCBR1245_data_results/cibersort/cibersort_gene_lists.rda")
  
missing.from.ref <- setdiff(ncounter.gene.list, lm22.reference.genes)

missing.from.ncounter <- setdiff(lm22.reference.genes, ncounter.gene.list)

found.in.ref <- intersect(ncounter.gene.list, lm22.reference.genes)


# Genes to map, from Bahman's analysis using geneSynonym
mapping.list <- c("EMR1" = "ADGRE1", 
                  "FZD3" = "FZD9", 
                  "KIR2DL1" = "KIR2DL3", 
                  "SIRPG" = "SIRPB2", 
                  "TEP1" = "PTEN")

norm.counts.mapped <- norm.counts %>%
  mutate(gene = ifelse(gene %in% mapping.list, names(mapping.list)[match(gene, mapping.list)], gene))

# Write out inout file for cibersort
write.table(norm.counts.mapped,
          paste0("/rstudio-files/ccbr-data/users/Ned/nCounter_rna/CCBR1245_data_results/cibersort/CCBR1245_ciber_sort_input.tsv"), sep = "\t")



```