Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Handling.qmd

---
title: "Data handling for Widefield PDGFR-β and GFAP-stained barain sections (Ipsilateral hemisphere)"
subtitle: "Data analysis notebook"
author: 
  - name: "Daniel Manrique-Castano"
    orcid: 0000-0002-1912-1764
    degrees:
      - PhD
    affiliation: 
      - name: Univerisity Laval 
        department: Psychiatry and Neuroscience
        group: Laboratory of neurovascular interactions 
note: "GitHub: https://daniel-manrique.github.io/"

keywords: 
  - Brain shrinkage
  - Brain injury
  - MCAO
  - Bayesian modeling 
   
license: "CC BY"

format:
   pdf: 
    toc: true
    number-sections: true
    colorlinks: true
   html:
    code-fold: true
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: left
    theme: spacelab

knitr:
  opts_chunk: 
    warning: false
    message: false

csl: science.csl
bibliography: references.bib
---

# Brains with cortico-striatal lesions

## Preview

In this notebook, we handle .tsv files to generate point patterns and summary tables for further analysis. 

**Parent dataset:** PDGFR-β and GFAP-stained ischemic hemispheres imaged at 5x (with stitching). Samples were grouped at 0 (Sham), 3, 7, 14, and 30 days post-ischemia (DPI).

**Working dataset**: .tsv files exported with [QuPath](https://qupath.github.io/) [@bankhead2017]. We performed unbiased detection and quantification of PDGFR-β and GFAP-positive cells in the **whole ipsilateral hemisphere**. QuPath generates `_detections.tsv` files containing the coordinates of individual cells and other measurements. Also, it creates `_annotations.tsv` files summarizing the information by image. These files are located in the `QuPath_Ipsilateral_5x` folder within the working directory.


## I. Install and load required packages

Install and load all required packages. Please uncomment (delete #) the line code if installation is required. Load the installed libraries each time you start a new R session. 
```{r}
#| label: Install_Packages
#| include: true
#| warning: false
#| message: false

#library(devtools)
#install.packages(c("dpylr", "tidyr", "readr","spatstat"))

library(tidyr)
library(dplyr)
library(readr)
library(spatstat)
library(spatstat.geom)

```



## Processing of image **annotations**

The following chunk handles the `_annotations.tsv` files of each image to obtain a single .csv file. The files for PDGFR-β and GFAP cells are processed independently, given that PDGFR-β files contain an additional column corresponding to the employed (QuPath) cell classifier.

### For GFAP

```{r}
#| label: Gfap_Ipsi_5x_Annotations
#| include: true
#| warning: false
#| message: false
#| cache: true

append_annotations <- function(base_path, results_path) {
  
  Gfap_csv_path <- paste0(results_path, "/Raw_Widefield_5x_Ipsilateral_Gfap_QuPathAnnotations.csv")
  Annotations_Path <- paste0(base_path, "/Gfap/")
  process_annotation(results_path = Gfap_csv_path, Annotations_Path)
}
  
process_annotation <- function(results_path, path) {
  
  print (path)
  
  Annotations <- list.files(path = path, pattern = "annotations.tsv", full.names = TRUE) %>% 
    lapply(read_tsv) %>%                              
    bind_rows
 
   print(Annotations)
  
  Annotations <- as.data.frame(Annotations)
  names(Annotations) <- NULL
  
  write.table(Annotations, results_path, append = TRUE, sep=",")
}


basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x"
resultsPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb"

Gfap_csv_path <- paste0(resultsPath, "/Raw_Widefield_5x_Ipsilateral_Gfap_QuPathAnnotations.csv")

Annotations_Header <- c("x", "Image",	"Object", "Name",	"Class",	"Parent",	"ROI",	"Centroid X",	"Centroid Y",	"Tissue_Area","Num Detections",	"Area^2",	"Perimeter")

df_header <- data.frame(matrix(ncol = 13, nrow = 0))
names(df_header) <- Annotations_Header

write.csv(df_header, Gfap_csv_path, row.names = FALSE)

brains <- list.dirs(basePath, full.names = FALSE, recursive = FALSE)

for (brain in brains){
  append_annotations(basePath, resultsPath)
  break
}
```

### For PDGFR-β

```{r}
#| label: Pdgfrb_Ipsi_5x_Annotations
#| include: true
#| warning: false
#| message: false


append_annotations <- function(base_path, results_path) {
  
  Pdgfrb_csv_path <- paste0(results_path, "/Raw_Widefield_5x_Ipsilateral_Pdgfrb_QuPathAnnotations.csv")
  Annotations_Path <- paste0(base_path, "/Pdgfrb/")
  process_annotation(results_path = Pdgfrb_csv_path, Annotations_Path)
}
  
process_annotation <- function(results_path, path) {
  
  print (path)
  
  Annotations <- list.files(path = path, pattern = "annotations.tsv", full.names = TRUE) %>% 
    lapply(read_tsv) %>%                              
    bind_rows
 
   print(Annotations)
  
  Annotations <- as.data.frame(Annotations)
  names(Annotations) <- NULL
  
  write.table(Annotations, results_path, append = TRUE, sep=",")
}


basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x"
resultsPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb"

Pdgfrb_csv_path <- paste0(resultsPath, "/Raw_Widefield_5x_Ipsilateral_Pdgfrb_QuPathAnnotations.csv")

Annotations_Header <- c("x", "Image",	"Object",	"Name", "Class",	"Parent",	"ROI",	"Centroid X",	"Centroid Y",	"Tisue_Area",	"Num Detections",	"Num Negatie", "Num Positive", "Positive %", "Num Positive per mm^2",  "Area^2",	"Perimeter")

df_header <- data.frame(matrix(ncol = 17, nrow = 0))
names(df_header) <- Annotations_Header

write.csv(df_header, Pdgfrb_csv_path, row.names = FALSE)

brains <- list.dirs(basePath, full.names = FALSE, recursive = FALSE)

for (brain in brains){
  append_annotations(basePath, resultsPath)
  break
}
```

Now, the results are stored as `Data_Raw/Raw_Widefield_5x_Ipsilateral_Gfap_QuPathAnnotations.csv` and `Data_Raw/Raw_Widefield_5x_Ipsilateral_Gfap_QuPathAnnotations.csv`.

## Processing of image **detections**

Next, we handle the `_detections.tsv` files containing the coordinates of individual cells to generate several `.csv` files with usable data.As previously, PDGFR-β and GFAP are processed independently. 

Please note that we include DAPI to obtain the limits of the observation window for the point patterns.

### For GFAP and DAPI

```{r}
#| label: Gfap_Ipsi_5x_Detections
#| include: true
#| warning: false
#| message: false
#| cache: true

process_initial_data <- function(basePath, Cells_Path, filename, resultsPath) {
  
  Cells_Raw <- read_tsv(paste0(basePath, "/", Cells_Path))
  
  # Convert to data frame
  Cells <- as.data.frame(Cells_Raw) 
  
  # Subset the date set to keep only relevant columns
  Cells <- subset(Cells, select = c(Image, Name, Parent, `Centroid X µm`, `Centroid Y µm`))
    
  Cells <- subset(Cells, Parent=="Positive", select = 1:5)
  
  # Extract metadata information from image name
  Cells <- cbind(Cells, do.call(rbind , strsplit(Cells$Image , "[_\\.]"))[,1:3])
  colnames(Cells) <- c("Image", "ID", "Class", "X", "Y", "MouseID", "DPI", "Region")
  Cells <- subset(Cells, select = c(MouseID, DPI, ID, X, Y))
  
  # Write a .csv file 
  write.csv(Cells, paste0(resultsPath, "/", Cells_Path, filename))
}
basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x"
resultsPath <- "Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_CellCoordinates"

process_folder <- function(folderPath, filename_suffix) {
  files <- list.files(folderPath, pattern = "_detections.tsv", full.names = FALSE)
  for (file in files) {
    process_initial_data(folderPath, file, filename_suffix, resultsPath)
  }
}
process_folder(paste0(basePath, "/Gfap"), "_Coordinates.csv")
process_folder(paste0(basePath, "/Dapi"), "_Coordinates.csv")
```

### For PDGFR-β

```{r}
#| label: Pdgfrb_Ipsi_5x_Detections
#| include: true
#| warning: false
#| message: false
#| cache: true

process_initial_data <- function(basePath, Cells_Path, filename, resultsPath) {
  
  Cells_Raw <- read_tsv(paste0(basePath, "/", Cells_Path))
  
  # Convert to data frame
  Cells <- as.data.frame(Cells_Raw) 
  
  # Subset the date set to keep only relevant columns
  Cells <- subset(Cells, select = c(Image, Name, Class, Parent, `Centroid X µm`, `Centroid Y µm`))
  Cells <- subset(Cells, Parent=="Positive")
  
  # Extract metadata information from image name
  Cells <- cbind(Cells, do.call(rbind , strsplit(Cells$Image , "[_\\.]"))[,1:3])
  colnames(Cells) <- c("Image", "ID", "Class","Parent", "X", "Y", "MouseID", "DPI", "Condition")
  Cells <- subset(Cells, select = c(MouseID, DPI, Condition, Class, X, Y))
  
  # Write a .csv file 
  write.csv(Cells, paste0(resultsPath, "/", Cells_Path, filename))
}
basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x"
resultsPath <- "Data_Raw/Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_CellCoordinates"

process_folder <- function(folderPath, filename_suffix) {
  files <- list.files(folderPath, pattern = "_detections.tsv", full.names = FALSE)
  for (file in files) {
    process_initial_data(folderPath, file, filename_suffix, resultsPath)
  }
}

process_folder(paste0(basePath, "/Pdgfrb"), "_Coordinates.csv")
```

Now, individual .csv files are stored in `Data_Raw/CellCoordinates/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_CellCoordinates`.
## Create Point Patterns

Using the the newly generated `_Coordinates.csv` files for PDGFR-β, GFAP, and DAPI, we create point patterns using the functions from the `spatstat` package [@spatstat; @baddeley2005]. Additionally, we generate a `.csv` containing the mean spatial intensity of GFAP and PDGFR-β+ cells. 

Here, we set the observation window using the DAPI images and set the scale to 3.407/3753 to re-scale the point patterns to mm. We also extract metadata from the file name to create identifiers for the point patterns. 

```{r}
#| label: Pdgfrb_Gfap_PPP
#| include: true
#| warning: false
#| message: false
#| cache: true

coordinatesPath <- "Data_Processed/Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_CellCoordinates/"

IntensityPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/"

Cells_Intensity_CSV_Path <- paste0(IntensityPath, "/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv")
Cells_Intensity_Header <- c("Brain", "PdgfrbHigh_Intensity", "PdgfrbLow_Intensity", "Gfap_Intensity")

# Results to generate
Result_Hyperframe <- NULL

# Functions

add_to_hyperframe <- function (...) {
    if (is.null(Result_Hyperframe)){
      Result_Hyperframe <<- hyperframe(...)
    } else {
      Result_Hyperframe <<- rbind(Result_Hyperframe, hyperframe(...))
    }
}

create_empty_table <- function (path, header) {
  df_header <- data.frame(matrix(ncol = length(header), nrow = 0))
  names(df_header) <- header

  write.csv(df_header, path)
}

create_empty_table(Cells_Intensity_CSV_Path, Cells_Intensity_Header)

coordinates_manipulation <- function (Raw_Table) {
  Cell_Coor_X <- Raw_Table$Y
  Cell_Coor_Y <- Raw_Table$X

  ## Bind the vectors, rotate and bind to original table
  Coords <- cbind(Cell_Coor_X, Cell_Coor_Y)
  Coords <- secr::rotate(Coords, 180)
  Coords <- as.data.frame(Coords)
  return(cbind(Raw_Table, Coords))
}

# Create a point pattern (PPP) object

create_point_pattern <- function(Subset, ReferenceSubset) {
  # We define the limits of the window according to Dapi coordinates
  xlim <- range(ReferenceSubset$X)
  ylim <- range(ReferenceSubset$Y)

  # Create point pattern for neurons
  Cells_PPP <- with(Subset, spatstat.geom::ppp(x = Subset$X, y = Subset$Y, xrange = xlim, yrange = ylim))
  unitname(Cells_PPP)  <- list("mm", "mm", 3.407/3753)
  Cells_PPP <- spatstat.geom::rescale (Cells_PPP)
  
  ## We rescale the unit to obtain measurements in mm2
  return(Cells_PPP)
  
}

process_file <- function (basePath, path) {

  Dapi_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Dapi_detections.tsv_Coordinates.csv'), header = TRUE)
  Dapi_Raw  <- Dapi_Raw  %>% sample_frac(.5)
  Gfap_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Gfap_detections.tsv_Coordinates.csv'), header = TRUE)
  Gfap_Raw  <- Gfap_Raw  %>% sample_frac(.5)
  Pdgfrb_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Pdgfrb_detections.tsv_Coordinates.csv'), header = TRUE)
  Pdgfrb_Raw  <- Pdgfrb_Raw  %>% sample_frac(.5)
  
  Dapi_Raw <- coordinates_manipulation(Dapi_Raw)
  Gfap_Raw <- coordinates_manipulation(Gfap_Raw)
  Pdgfrb_Raw <- coordinates_manipulation(Pdgfrb_Raw)
  
  # Subset cells
  Pdgfrb_High <- Pdgfrb_Raw[(Pdgfrb_Raw$Class=="Positive"),]
  Pdgfrb_Low <- Pdgfrb_Raw[(Pdgfrb_Raw$Class=="Negative"),]
  
  Dapi_PPP <- create_point_pattern(Dapi_Raw, Dapi_Raw)
  Window(Dapi_PPP) <- convexhull(Dapi_PPP)
  Pdgfrb_High_PPP <- create_point_pattern(Pdgfrb_High, Dapi_Raw)
  Window(Pdgfrb_High_PPP) <- convexhull(Dapi_PPP)
  Pdgfrb_Low_PPP <- create_point_pattern(Pdgfrb_Low, Dapi_Raw)
  Window(Pdgfrb_Low_PPP) <- convexhull(Dapi_PPP)
  Gfap_PPP <- create_point_pattern(Gfap_Raw, Dapi_Raw)
  Window(Gfap_PPP) <- convexhull(Dapi_PPP)

  Pdgfrb_High_Intensity <- summary(Pdgfrb_High_PPP)$intensity
  Pdgfrb_Low_Intensity <- summary(Pdgfrb_Low_PPP)$intensity
  Gfap_Intensity <- summary(Gfap_PPP)$intensity
  
  Intensity_Row <- t(c(path, Pdgfrb_High_Intensity, Pdgfrb_Low_Intensity, Gfap_Intensity))
    write.table(Intensity_Row, Cells_Intensity_CSV_Path, append = TRUE, sep=",", col.names = FALSE)
  
  fragments <- strsplit(path, "_")[[1]]
  len <- length(fragments)
  mouse <- fragments[1]
  dpi <- fragments[2]
  condition <- fragments[3]

add_to_hyperframe(Pdgfrb_High = Pdgfrb_High_PPP, Pdgfrb_Low = Pdgfrb_Low_PPP, Gfap = Gfap_PPP, ID = mouse, DPI=dpi, Condition = condition, stringsAsFactors=TRUE)
}

csv_files <- list.files(coordinatesPath, full.names = FALSE, recursive = FALSE)

brains <- c()

for (csv in csv_files) {
  fragments <- strsplit(csv, "_")[[1]]
  brain_name <- paste(fragments[1:3], collapse="_")
  brains <- append(brains, brain_name)
}

brains <- unique(brains)

for (brain in brains) {
  process_file(coordinatesPath, brain)
}

saveRDS(Result_Hyperframe, "PointPatterns/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_PPP.rds")
```

The point patterns are stored as a hyperframe object in `PointPatterns/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_PPP.rds`, and the file with the summary information in `Data_Raw/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv`. The .rds file hyperframe is structured by columns the following way:

Pdgfr_High: Point patterns for PDGFR-β_High-expressing cells.
Pdgfr_Low: Point patterns for PDGFR-β_low-expressing cells.
Gfap = GFAP-expressing cells.
ID = Mouse unique identifier.
DPI = Days post-ischemia 0 (for sham), 3, 7, 14, or 30 DPI.
Condition = Experimental procedure (Sham or MCAO).

## Handle dataset for analysis

Our reference file in this section is `Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv` generated during the creation of the point patterns. From this file, we obtain the grouping variables and other metadata consigned in the file names. 

```{r}
#| label: Pdgfrb_Gfap_Ipsi_5x_table
#| include: true
#| warning: false
#| message: false

# Read the data tables
Ipsilateral_Intensity_5x <- read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv', header = TRUE)

# We subset the relevant columns
Ipsilateral_Intensity_5x <-  Ipsilateral_Intensity_5x[, -1]

# Extract metadata from the image name and rename columns
Ipsilateral_Intensity_5x <- cbind(Ipsilateral_Intensity_5x, do.call(rbind , strsplit(Ipsilateral_Intensity_5x$Brain, "[_\\.]"))[,1:3])
colnames(Ipsilateral_Intensity_5x) <- c("Brain", "PdgfrbHigh_Intensity", "PdgfrbLow_Intensity", "Gfap_Intensity", "MouseID", "DPI", "Condition")

# We subset the relevant columns
Ipsilateral_Intensity_5x<- subset(Ipsilateral_Intensity_5x, select = c(MouseID, DPI, Condition, PdgfrbHigh_Intensity, PdgfrbLow_Intensity, Gfap_Intensity))

# Convert DPI to numeric
DPI_mapping <- c("0D" = "0", "3D" = "3", "7D" = "7", "14D" = "14", "30D" = "30")
Ipsilateral_Intensity_5x$DPI <- as.numeric(DPI_mapping[as.character(Ipsilateral_Intensity_5x$DPI)])

# We will set the Sham condition as 0 DPI
Ipsilateral_Intensity_5x$DPI[Ipsilateral_Intensity_5x$Condition == "SHAM"] <- 0

# Set explicitly the factor levels for Condition
Ipsilateral_Intensity_5x$Condition <- factor(Ipsilateral_Intensity_5x$Condition, levels = c("SHAM", "MCAO"))

```

Now, we merge the generated `Gfap_Annotations.csv` and `Pdgfr_Annotations.csv` and subset the relevant data.

```{r}
#| label: Pdgfrb_Gfap_Ipsi_5x_Subset
#| include: true
#| warning: false
#| message: false

Gfap_Annotations_Raw <- 
  read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Gfap_QuPathAnnotations.csv', 
           header = TRUE
           )
Pdgfrb_Annotations_Raw <- 
  read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Pdgfrb_QuPathAnnotations.csv', 
           header = TRUE
           )
# We select the rows containing "Positive" tissue
Gfap_Annotations_Raw  <- 
  Gfap_Annotations_Raw [(Gfap_Annotations_Raw$Class=="Positive"),]

Pdgfrb_Annotations_Raw  <- 
  Pdgfrb_Annotations_Raw [(Pdgfrb_Annotations_Raw$Class=="Positive"),]
```

Next, we merged the intensity measurements with the annotations files.

```{r}
#| label: Pdgfr_Gfap_Ipsi_5x_Merge
#| include: true
#| warning: false
#| message: false

# Bind the Gfap dataset
Ipsilateral_Summary_5x <- cbind(Ipsilateral_Intensity_5x, Gfap_Annotations_Raw[, c(9, 10)])
names(Ipsilateral_Summary_5x)[names(Ipsilateral_Summary_5x) == 'Num.Detections'] <- 'Gfap_Num'

# Bind the Pdgfrb dataset
column_names <- c('Pdgfrb_Total', 'Pdgfrb_Neg', 'Pdgfrb_Pos', 'Pdgfrb_PercPositive')
Ipsilateral_Summary_5x <- cbind(Ipsilateral_Summary_5x, Pdgfrb_Annotations_Raw[, 10:13])
names(Ipsilateral_Summary_5x)[names(Ipsilateral_Summary_5x) %in% c('Num.Detections', 'Num.Negatie', 'Num.Positive', 'Positive..')] <- column_names

# Scale the Tissue_Area
Ipsilateral_Summary_5x$Tissue_Area <- Ipsilateral_Summary_5x$Tissue_Area / 1000000

```

Finally, we include the Area and intensity measurements performed in FIJI

```{r}
#| label: Pdgfr_Gfap_Ipsi_5x_Fiji
#| include: true
#| warning: false
#| message: false

# Read the data table
Ipsilateral_AreaIntensity_Summary_Raw <- read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_AreaIntensity.csv', header = TRUE)

# Extract Gfap data
Gfap_Area_Raw <- 
  Ipsilateral_AreaIntensity_Summary_Raw[Ipsilateral_AreaIntensity_Summary_Raw$Label == "Gfap", ]
Gfap_Columns <- c("Gfap_Area", "Gfap_MeanInt", "Gfap_sdInt", "Gfap_IntDen")

Ipsilateral_Summary_5x <- 
  cbind(Ipsilateral_Summary_5x, Gfap_Area_Raw[, 3:6])

names(Ipsilateral_Summary_5x)[names(Ipsilateral_Summary_5x) %in% c("Area", "Mean", "StdDev", "IntDen")] <- Gfap_Columns

# Extract Pdgfrb data
Pdgfrb_Area_Raw <- 
  Ipsilateral_AreaIntensity_Summary_Raw[Ipsilateral_AreaIntensity_Summary_Raw$Label == "Pdgfr", ]

Pdgfrb_Columns <- c("Pdgfrb_Area", "Pdgfbr_MeanInt", "Pdgfrb_sdInt", "Pdgfrb_IntDen")

Ipsilateral_Summary_5x <- 
  cbind(Ipsilateral_Summary_5x, Pdgfrb_Area_Raw[, 3:6])
names(Ipsilateral_Summary_5x)[names(Ipsilateral_Summary_5x) %in% c("Area", "Mean", "StdDev", "IntDen")] <- Pdgfrb_Columns

# Scale the required columns
Ipsilateral_Summary_5x$Gfap_Area <- Ipsilateral_Summary_5x$Gfap_Area /1000000
Ipsilateral_Summary_5x$Gfap_IntDen <- Ipsilateral_Summary_5x$Gfap_IntDen /100000
Ipsilateral_Summary_5x$Pdgfrb_Area <- Ipsilateral_Summary_5x$Pdgfrb_Area /1000000
Ipsilateral_Summary_5x$Pdgfrb_IntDen <- Ipsilateral_Summary_5x$Pdgfrb_IntDen /100000

head(Ipsilateral_Summary_5x)

write.csv(Ipsilateral_Summary_5x, "Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Inten.csv", row.names = FALSE)
```

These whole data handling results in a dataset organized by columns, which will be used for analysis. 

# Brains with (only) striatal lesions

In this section, we perform the same procedure for brains comprising striatal lesions. This data set is processed separately because brain sections were stained in a different day than those comprising cortico-striatal lesions. 

**Parent dataset:** PDGFR-β and GFAP-stained ischemic hemispheres imaged at 5x (with stitching). Samples were grouped at 14, and 30 days post-ischemia (DPI).

**Working dataset**: .tsv files exported with [QuPath](https://qupath.github.io/) [@bankhead2017]. We performed unbiased detection and quantification of PDGFR-β and GFAP-positive cells in the **whole ipsilateral hemisphere**. QuPath generates `_detections.tsv` files containing the coordinates of individual cells and other measurements. Also, it creates `_annotations.tsv` files summarizing the information by image. These files are located in the `QuPath_Ipsilateral_5x_Striatal` folder within the working directory.

## Processing of image **annotations**

### For GFAP

```{r}
#| label: Gfap_Annotations_Striatum_5x
#| include: true
#| warning: false
#| message: false
#| cache: true

append_annotations <- function(base_path, results_path) {
  
  Gfap_csv_path <- paste0(resultsPath, "/Raw_Widefield_5x_Ipsilateral_Gfap_Str_QuPathAnnotations.csv")
  Annotations_Path <- paste0(base_path, "/Gfap/")
  process_annotation(results_path = Gfap_csv_path, Annotations_Path)
  
}

process_annotation <- function(results_path, path) {
  
  print (path)
  
  Annotations <- list.files(path = path, pattern = "_annotations.tsv", full.names = TRUE) %>% 
    lapply(read_tsv) %>%                              
    bind_rows
 
   print(Annotations)
  
  Annotations <- as.data.frame(Annotations)
  names(Annotations) <- NULL
  
  write.table(Annotations, results_path, append = TRUE, sep=",")
}

basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x_Striatal"
resultsPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb"

Gfap_csv_path <- paste0(resultsPath, "/Raw_Widefield_5x_Ipsilateral_Gfap_Str_QuPathAnnotations.csv")

Annotations_Header <- c("x", "Image",	"Object", "Name",	"Class",	"Parent",	"ROI",	"Centroid X ?m",	"Centroid Y ?m",	"Tissue_Area","Num Detections",	"Area ?m^2",	"Perimeter ?m")

df_header <- data.frame(matrix(ncol = 13, nrow = 0))
names(df_header) <- Annotations_Header

write.csv(df_header, Gfap_csv_path, row.names = FALSE)

brains <- list.dirs(basePath, full.names = FALSE, recursive = FALSE)

for (brain in brains){
  append_annotations(basePath, resultsPath)
  break
 
}
```

### For PDGFR-β

```{r}
#| label: PDGFR-β_Annotations_Striatum_5x
#| include: true
#| warning: false
#| message: false


append_annotations <- function(base_path, results_path) {
  
  Pdgfrb_csv_path <- paste0(results_path, "/Raw_Widefield_5x_Ipsilateral_Pdgfrb_Str_QuPathAnnotations.csv")
  Annotations_Path <- paste0(base_path, "/Pdgfrb/")
  process_annotation(results_path = Pdgfrb_csv_path, Annotations_Path)
}
  
process_annotation <- function(results_path, path) {
  
  print (path)
  
  Annotations <- list.files(path = path, pattern = "annotations.tsv", full.names = TRUE) %>% 
    lapply(read_tsv) %>%                              
    bind_rows
 
   print(Annotations)
  
  Annotations <- as.data.frame(Annotations)
  names(Annotations) <- NULL
  
  write.table(Annotations, results_path, append = TRUE, sep=",")
}


basePath <- "QuPath_Ipsilateral_5x_Striatal"
resultsPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/"

Pdgfr_csv_path <- paste0(resultsPath, "/Raw_Widefield_5x_Ipsilateral_Pdgfrb_Str_QuPathAnnotations.csv")

Annotations_Header <- c("x", "Image",	"Object",	"Name", "Class",	"Parent",	"ROI",	"Centroid X ?m",	"Centroid Y ?m",	"Tisue_Area",	"Num Detections",	"Num Negatie", "Num Positive", "Positive %", "Num Positive per mm^2",  "Area ?m^2",	"Perimeter ?m")

df_header <- data.frame(matrix(ncol = 17, nrow = 0))
names(df_header) <- Annotations_Header

write.csv(df_header, Pdgfr_csv_path, row.names = FALSE)

brains <- list.dirs(basePath, full.names = FALSE, recursive = FALSE)

for (brain in brains){
  append_annotations(basePath, resultsPath)
  break
}
```

## Batch processing of image **detections**

### For GFAP and DAPI

```{r}
#| label: GFAP_DAPI_Detections_Striatum_5x
#| include: true
#| warning: false
#| message: false
#| cache: true

process_initial_data <- function(basePath, Cells_Path, filename, resultsPath) {
  
  Cells_Raw <- read_tsv(paste0(basePath, "/", Cells_Path))
  
  # Convert to data frame
  Cells <- as.data.frame(Cells_Raw) 
  
  # Subset the date set to keep only relevant columns
  Cells <- subset(Cells, select = c(Image, Name, Parent, `Centroid X µm`, `Centroid Y µm`))
  Cells <- subset(Cells, Parent=="Positive", select = 1:5)
  
  # Extract metadata information from image name
  Cells <- cbind(Cells, do.call(rbind , strsplit(Cells$Image , "[_\\.]"))[,1:3])
  colnames(Cells) <- c("Image", "ID", "Class", "X", "Y", "MouseID", "DPI", "Region")
  Cells <- subset(Cells, select = c(MouseID, DPI, ID, X, Y))
  
  # Write a .csv file 
  write.csv(Cells, paste0(resultsPath, "/", Cells_Path, filename))
}
basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x_Striatal"
resultsPath <- "Data_Raw/CellCoordinates/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_Striatum_CellCoordinates"

process_folder <- function(folderPath, filename_suffix) {
  files <- list.files(folderPath, pattern = "_detections.tsv", full.names = FALSE)
  for (file in files) {
    process_initial_data(folderPath, file, filename_suffix, resultsPath)
  }
}
process_folder(paste0(basePath, "/Gfap"), "_Coordinates.csv")
process_folder(paste0(basePath, "/Dapi"), "_Coordinates.csv")
```

### For PDGFR-β

```{r}
#| label: PDGFR-β_Detections_5x
#| include: true
#| warning: false
#| message: false
#| cache: true

process_initial_data <- function(basePath, Cells_Path, filename, resultsPath) {
  
  Cells_Raw <- read_tsv(paste0(basePath, "/", Cells_Path))
  
  # Convert to data frame
  Cells <- as.data.frame(Cells_Raw) 
  
  # Subset the date set to keep only relevant columns
  Cells <- subset(Cells, select = c(Image, Name, Class, Parent, `Centroid X µm`, `Centroid Y µm`))
  Cells <- subset(Cells, Parent=="Positive")
  
  # Extract metadata information from image name
  Cells <- cbind(Cells, do.call(rbind , strsplit(Cells$Image , "[_\\.]"))[,1:3])
  colnames(Cells) <- c("Image", "ID", "Class","Parent", "X", "Y", "MouseID", "DPI", "Condition")
  Cells <- subset(Cells, select = c(MouseID, DPI, Condition, Class, X, Y))
  
  # Write a .csv file 
  write.csv(Cells, paste0(resultsPath, "/", Cells_Path, filename))
}
basePath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/QuPath_Ipsilateral_5x_Striatal"
resultsPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_Striatum_CellCoordinates"

process_folder <- function(folderPath, filename_suffix) {
  files <- list.files(folderPath, pattern = "_detections.tsv", full.names = FALSE)
  for (file in files) {
    process_initial_data(folderPath, file, filename_suffix, resultsPath)
  }
}

process_folder(paste0(basePath, "/Pdgfrb"), "_Coordinates.csv")
```

## Create Point Patterns for striatal lesions

Like done previously, we generate point patterns using the functions from the `spatstat` package [@spatstat; @baddeley2005]. 

```{r}
#| label: PDGFR-β_GFAP_Striatum_PPP
#| include: true
#| warning: false
#| message: false
#| cache: true

coordinatesPath <- "Data_Raw/CellCoordinates/Widefield_5x_Ipsilateral_Gfap-Pdgfrb-Dapi_Striatum_CellCoordinates/"

IntensityPath <- "Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb"

Cells_Intensity_CSV_Path <- paste0(IntensityPath, "/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Striatum_Inten.csv")
Cells_Intensity_Header <- c("Brain", "PdgfrbHigh_Intensity", "PdgfrbLow_Intensity", "Gfap_Intensity")

# Results to generate
Result_Hyperframe <- NULL

# Functions

add_to_hyperframe <- function (...) {
    if (is.null(Result_Hyperframe)){
      Result_Hyperframe <<- hyperframe(...)
    } else {
      Result_Hyperframe <<- rbind(Result_Hyperframe, hyperframe(...))
    }
}

create_empty_table <- function (path, header) {
  df_header <- data.frame(matrix(ncol = length(header), nrow = 0))
  names(df_header) <- header

  write.csv(df_header, path)
}

create_empty_table(Cells_Intensity_CSV_Path, Cells_Intensity_Header)

coordinates_manipulation <- function (Raw_Table) {
  Cell_Coor_X <- Raw_Table$Y
  Cell_Coor_Y <- Raw_Table$X

  ## Bind the vectors, rotate and bind to original table
  Coords <- cbind(Cell_Coor_X, Cell_Coor_Y)
  Coords <- secr::rotate(Coords, 180)
  Coords <- as.data.frame(Coords)
  return(cbind(Raw_Table, Coords))
}

# Create a point pattern (PPP) object

create_point_pattern <- function(Subset, ReferenceSubset) {
  # We define the limits of the window according to Dapi coordinates
  xlim <- range(ReferenceSubset$X)
  ylim <- range(ReferenceSubset$Y)

  # Create point pattern for neurons
  Cells_PPP <- with(Subset, spatstat.geom::ppp(x = Subset$X, y = Subset$Y, xrange = xlim, yrange = ylim))
  unitname(Cells_PPP)  <- list("mm", "mm", 3.407/3753)
  Cells_PPP <- spatstat.geom::rescale (Cells_PPP)
  
  ## We rescale the unit to obtain measurements in mm2
  return(Cells_PPP)
  
}

process_file <- function (basePath, path) {

  Dapi_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Dapi_detections.tsv_Coordinates.csv'), header = TRUE)
  Dapi_Raw  <- Dapi_Raw  %>% sample_frac(.5)
  Gfap_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Gfap_detections.tsv_Coordinates.csv'), header = TRUE)
  Gfap_Raw  <- Gfap_Raw  %>% sample_frac(.5)
  Pdgfrb_Raw <- read.csv(file = paste0(basePath, '/', path, '_Str_Pdgfrb_detections.tsv_Coordinates.csv'), header = TRUE)
  Pdgfrb_Raw  <- Pdgfrb_Raw  %>% sample_frac(.5)
  
  Dapi_Raw <- coordinates_manipulation(Dapi_Raw)
  Gfap_Raw <- coordinates_manipulation(Gfap_Raw)
  Pdgfrb_Raw <- coordinates_manipulation(Pdgfrb_Raw)
  
  # Subset neurons
  Pdgfrb_High <- Pdgfrb_Raw[(Pdgfrb_Raw$Class=="Positive"),]
  Pdgfrb_Low <- Pdgfrb_Raw[(Pdgfrb_Raw$Class=="Negative"),]
  
  Dapi_PPP <- create_point_pattern(Dapi_Raw, Dapi_Raw)
  Window(Dapi_PPP) <- convexhull(Dapi_PPP)
  Pdgfrb_High_PPP <- create_point_pattern(Pdgfrb_High, Dapi_Raw)
  Window(Pdgfrb_High_PPP) <- convexhull(Dapi_PPP)
  Pdgfrb_Low_PPP <- create_point_pattern(Pdgfrb_Low, Dapi_Raw)
  Window(Pdgfrb_Low_PPP) <- convexhull(Dapi_PPP)
  Gfap_PPP <- create_point_pattern(Gfap_Raw, Dapi_Raw)
  Window(Gfap_PPP) <- convexhull(Dapi_PPP)

  Pdgfrb_High_Intensity <- summary(Pdgfrb_High_PPP)$intensity
  Pdgfrb_Low_Intensity <- summary(Pdgfrb_Low_PPP)$intensity
  Gfap_Intensity <- summary(Gfap_PPP)$intensity
  
  Intensity_Row <- t(c(path, Pdgfrb_High_Intensity, Pdgfrb_Low_Intensity, Gfap_Intensity))
    write.table(Intensity_Row, Cells_Intensity_CSV_Path, append = TRUE, sep=",", col.names = FALSE)
  
  fragments <- strsplit(path, "_")[[1]]
  len <- length(fragments)
  mouse <- fragments[1]
  dpi <- fragments[2]
  condition <- fragments[3]

add_to_hyperframe(Pdgfrb_High = Pdgfrb_High_PPP, Pdgfrb_Low = Pdgfrb_Low_PPP, Gfap = Gfap_PPP, ID = mouse, DPI=dpi, Condition = condition, stringsAsFactors=TRUE)
}

csv_files <- list.files(coordinatesPath, full.names = FALSE, recursive = FALSE)

brains <- c()

for (csv in csv_files) {
  fragments <- strsplit(csv, "_")[[1]]
  brain_name <- paste(fragments[1:3], collapse="_")
  brains <- append(brains, brain_name)
}

brains <- unique(brains)

for (brain in brains) {
  process_file(coordinatesPath, brain)
}

saveRDS(Result_Hyperframe, "PointPatterns/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Str_PPP.rds")
```

## Buid dataset for analysis

Here we take the `Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Striatum_Inten.csv` generated during the creation of point patterns as a reference file.

```{r}
#| label: DataSet_Ipsilateral_Striatum_5x_1
#| include: true
#| warning: false
#| message: false

# Read the data tables
Ipsilateral_Striatum_Intensity_5x <- read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Striatum_Inten.csv', header = TRUE)

# We subset the relevant columns
Ipsilateral_Striatum_Intensity_5x <-  Ipsilateral_Striatum_Intensity_5x[, -1]

# Extract metadata from the image name and rename columns
Ipsilateral_Striatum_Intensity_5x <- cbind(Ipsilateral_Striatum_Intensity_5x, do.call(rbind , strsplit(Ipsilateral_Striatum_Intensity_5x$Brain, "[_\\.]"))[,1:3])
colnames(Ipsilateral_Striatum_Intensity_5x) <- c("Brain", "PdgfrbHigh_Intensity", "PdgfrbLow_Intensity", "Gfap_Intensity", "MouseID", "DPI", "Condition")

# We subset the relevant columns
Ipsilateral_Striatum_Intensity_5x<- subset(Ipsilateral_Striatum_Intensity_5x, select = c(MouseID, DPI, Condition, PdgfrbHigh_Intensity, PdgfrbLow_Intensity, Gfap_Intensity))

# Convert DPI to numeric
DPI_mapping <- c("14D" = "14", "30D" = "30")
Ipsilateral_Striatum_Intensity_5x$DPI <- as.numeric(DPI_mapping[as.character(Ipsilateral_Striatum_Intensity_5x$DPI)])
```

Now, we merge the generated `GFAP_Annotations.csv` and `PDGFR_Annotations.csv` and subset the relevant data.

```{r}
#| label: DataSet_Ipsilateral_Striatum_5x_2 
#| include: true
#| warning: false
#| message: false

Gfap_Annotations_Raw <- 
  read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Gfap_Str_QuPathAnnotations.csv', 
           header = TRUE
           )
Gfap_Annotations_Raw <- Gfap_Annotations_Raw[, -1]

Pdgfrb_Annotations_Raw <- 
  read.csv(file = 'Data_Raw/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Raw_Widefield_5x_Ipsilateral_Pdgfrb_Str_QuPathAnnotations.csv', 
           header = TRUE
           )
Pdgfrb_Annotations_Raw <- Pdgfrb_Annotations_Raw[, -1]

# We select the rows containing "Positive" tissue
Gfap_Annotations_Raw  <- 
  Gfap_Annotations_Raw [(Gfap_Annotations_Raw$Class=="Positive"),]

Pdgfr_Annotations_Raw  <- 
  Pdgfrb_Annotations_Raw [(Pdgfr_Annotations_Raw$Class=="Positive"),]
```

Next, we merged the intensity measurements with the annotations files.

```{r}
#| label: DataSet_Ipsilateral_Striatum_5x (3)
#| include: true
#| warning: false
#| message: false

# Bind the Gfap dataset
Ipsilateral_Striatum_Summary_5x <- cbind(Ipsilateral_Striatum_Intensity_5x, Gfap_Annotations_Raw[, c(9, 10)])
names(Ipsilateral_Striatum_Summary_5x)[names(Ipsilateral_Striatum_Summary_5x) == 'Num.Detections'] <- 'Gfap_Num'

# Bind the Pdgfrb dataset
column_names <- c('Pdgfrb_Total', 'Pdgfrb_NonReact', 'Pdgfrb_React', 'Pdgfrb_PercPositive')
Ipsilateral_Striatum_Summary_5x <- cbind(Ipsilateral_Striatum_Summary_5x, Pdgfrb_Annotations_Raw[, 10:13])
names(Ipsilateral_Striatum_Summary_5x)[names(Ipsilateral_Striatum_Summary_5x) %in% c('Num.Detections', 'Num.Negatie', 'Num.Positive', 'Positive..')] <- column_names

# Scale the Tissue_Area
Ipsilateral_Striatum_Summary_5x$Tissue_Area <- Ipsilateral_Striatum_Summary_5x$Tissue_Area / 1000000

write.csv(Ipsilateral_Striatum_Summary_5x, "Data_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Str_Inten.csv", row.names = FALSE)
```

Now, the complete dataset from this staining is complete and saved in `Data_Tables_Processed/Widefield_5x_Ipsilateral_Gfap-Pdgfrb_Striatum_Inten.csv`.

# References

::: {#refs}
:::

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