P6.2-wound_response_GeoMx_spatial.Rmd

---
title: "Wound-response CAF score on stromal spots of GeoMx RNA assay"
author: "Josep Garnica"
date: "`r format(Sys.Date(),'%e de %B, %Y')`"
output: 
  rmdformats::downcute:
    lightbox: true
    thumbnails: false
    self_contained: true
    gallery: true
    code_folding: show
  pkgdown:
    as_is: true
---

The goal of this script is to evaluate the wound-responding signature in CAF obtained in script `P5-wounding_TME_CAFsignature.Rmd` in GeoMx-RNA-derived stromal spots.

GeoMx data data was obtained and produced by [Yerly et al., 2022](https://www.nature.com/articles/s41467-022-32670-w) and can be accessed on [**GSE210648** on GEO](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE210648)


Load needed packages
```{r Load_packages, include=FALSE}
renv::restore()
library(ggplot2)
library(UCell)
library(tidyr)
library(dplyr)
library(tibble)
library(ggprism)
```


# Download data
We will download SOFT formatted family file, which include processed counts for the GeoMx counts. Next we will adapt this files into a matrix, prepared for our analysis. 

Processed file can be found [here](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL32550)

```{r download_data}
# create directory to store data
ddir <- "cache"
dir.create(ddir)

#download file
dfile <- file.path(ddir, "GPL32550_family.soft.gz")
soft_link <- "https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL32nnn/GPL32550/soft/GPL32550_family.soft.gz"
download.file(url = soft_link,
              destfile = dfile)
```

We adapted Suppl. table 1A from [Yerly et al., 2022](https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-022-32670-w/MediaObjects/41467_2022_32670_MOESM4_ESM.xlsx) to get the type of BCC: either Nodular of infiltrative based on histological diagnosis.

```{r get_metadata}
ann <- read.csv("_aux/BCC_GeoMx_spots_labels_type.csv")
```


# Adapt file into a matrix for further processing

```{r adapt_softfile}
require(purrr)

soft <- readLines(dfile) 

# ids indicating that table begings
tb_begin <- grep("\\!sample_table_begin", soft)
tb_end <- grep("\\!sample_table_end", soft)
sample_tit <- grep("\\!Sample_title", soft)


li <- list()
for(r in 1:length(tb_begin)){
  # acomodat title sample
  tit <- soft[sample_tit[r]]
  tit <- strsplit(tit, " ")[[1]][7]
  
  # build table
  tb <- soft[(tb_begin[r]+2):(tb_end[r]-1)]
  tb <- as.data.frame(tb) %>% 
    separate(tb, sep = "\t", into = c("gene", tit))
  # convert to numeric
  tb[[2]] <- as.numeric(tb[[2]])
  
  li[[tit]] <- tb
}

# join by gene
mat <- li %>%
  purrr::reduce(full_join, by = "gene") %>% 
  column_to_rownames("gene") %>%
  as.matrix()

```

## Filter only stromal spots

```{r filter_stromal}
# adapt matrix names, as they are not consistent
colnames(mat) <- gsub("[-]", "_", colnames(mat))

# produce comprehensive metadata based on colnames of the matrix and other metadata
md <- data.frame(sample = colnames(mat)) %>% 
  separate(sample, sep = "_", into = c("ID", "AOI", "type"), remove = F) %>% 
  # join with metadata to have the BCC type
  left_join(., ann, by = "ID") %>% 
  filter(!is.na(AOI)) %>% 
  # further adapt spots names for consistency
  mutate(type = ifelse(grepl("troma", type, ignore.case = T),
                "stroma", 
                ifelse(is.na(type),
                       "CK", type)
                ) 
         ) %>% 
  # get only stroma spots
  filter(type == "stroma") 

mat2 <- mat[, md$sample]
```


# Wound-response score on stromal spots
## Load wound-response CAF signatures

Load signatures previously produced in `P5-wounding_TME_CAFsignature.Rmd`.

```{r Load_signatures}
# Load file with signatures
lys <- readRDS("cache/CAF_signatures.rds")
```


## UCell scoring

### Run UCell
```{r UCell_scoring}
# run scoring
u.scores <- ScoreSignatures_UCell(matrix = mat2,
                                  features = lys,
                                  ncores = 4)
```

### Compute wound-response scoring
The “wound-responding CAF” score was calculated by subtracting the unwound/baseline signature score from the wound signature score in each spot.

```{r process_scoring}
# adapt resulting scores for ploting and statistical analysis
u.scores <- u.scores %>%
  as.data.frame() %>% 
  rownames_to_column("sample") %>% 
  mutate(CAF_score = Wound_UCell - Unwound_UCell) %>% 
  left_join(., md, by = "sample")

# acommodate Diagnosis or BCC type
u.scores <-  u.scores %>% 
  mutate(Diagnosis = factor(Diagnosis,
                            levels = c("Nodular","Infiltrative"),
                            labels = c("Non-invasive BCC","Invasive BCC")))
```

## Plot results

```{r plot_CAF_score}
# set colors
lycols <- c("Non-invasive BCC" = "#FFFF6C",
            "Invasive BCC" = "#F44400")

# fix seed
set.seed(22)

# do plot
pl <- u.scores %>% 
  ggplot(aes(Diagnosis, CAF_score, fill = Diagnosis)) +
  geom_boxplot(outlier.colour = NA,
               show.legend = F,
               width = 0.5) +
  geom_jitter(width = 0.1,
              size = 3,
              shape = 21,
              fill = "white",
              show.legend = F) +
  scale_fill_manual(values = lycols) +
  ggprism::theme_prism() +
  labs(title = "Intratumoral stromal\nGeoMx spots",
       y = "Wound-response CAF score",
       x = "")
pl
```

## Stats analysis
```{r t-test_comparison}
# Separate the values for each category
inf <- u.scores$CAF_score[u.scores$Diagnosis == "Invasive BCC"]
nod <- u.scores$CAF_score[u.scores$Diagnosis == "Non-invasive BCC"]

# Perform t-test
t.test(inf, nod)
```