NerveHotspotDetector.Rmd

---
title: "Development of a computational framework to detect and quantify nerve hotspots"
author: "Dimitrios Kleftogiannis"
date: "2024-07-05"
output: html_document
---

### Utility

The utility of this code is to process data from  mass cytometry imaging technologies and identify areas in the images that represent possible nerve elements. 

Such elements will be called "hotspots" and their presence is supported by the expression of nerve specific antibodies. In the analysis presented below we used peripherin as a nerve marker, but any other relevant nerve marker can be used.

To generalise the utility of this code, it essentially merging nuclei-based information of single-cells with any pixel-based data from the same image that have been processed individually. This is useful when we want to combine different elements in the images that go beyond the standard single-cell analysis.  

To be able to follow the pipeline we will provide some toy data and we will generate some basic visualisations and outputs that summarise the nerve hotspots.

The cohort data provided are courtesy of CCBIO, they are breast cancer tissue biopsies and contact person for the full dataset is Dr Kenneth Finne (kenneth.finne@uib.no) 

### Contact

Comments and bug reports are welcome. 
Please email: Dr Dimitrios Kleftogiannis (dimitrios.kleftogiannis@uib.no)

We are also interested to know about how you have used our source code, including any improvements that you have implemented.
 
You are free to modify, extend or distribute our source code, as long as our copyright notice remains unchanged and included in its entirety. 

### License

This code is licensed under the MIT License.

Copyright 2024, NeuroSysMed centre of clinical treatment research, University of Bergen (UiB), Norway

### Load packages. 
We load some packages required for the analysis. Please make sure that these packages are installed either via CRAN or from their relevant resources (github or Bioconductor).

```{r load packages, echo=FALSE, eval=TRUE, error=FALSE, warning=FALSE,cache=TRUE}
library("EBImage")
library("cytomapper")
library("ggplot2")
library(SingleCellExperiment)
library(RColorBrewer)
library(ggforce)
library(SpatialExperiment)
library(imcRtools)
library(sp)
```

### Load the data required and process the images

To run the pipeline is it necessary to execute externally nerve-specific segmentation using the Ilastik or any other compatible pipeline and segmentation masks in tiff format must be provided in a separate folder.

It is also recommended to perform and complete nuclei segmentation and single-cell data must processed and phenotype before running this pipeline. We recommend the Steinbock pipeline for single-cell data processing, and the single-cell data must be provided as an input here from an external file in a spatial experiment  format (saved as rds object).

Also note that in order to match nerve segmentation with nuclei segmentation the images must be indexed with the same name. Otherwise the pipeline will fail to run. 

Please check and download the folder named "input_data". 

For bulk processing of the images we also recommend to use a simple data file in txt format that  lists all file names of the nerve segmentation masks to be processed (see file_list.txt)  

```{r set wd, read and process data, echo=FALSE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}

#set manually the working dir --> this has to be changed once you download the git repo locally in your computer
myWorkDir <- paste('/Users/kleftogi/Desktop/IMC_new_cohort/NerveHostspotDetector/',sep='')
setwd(myWorkDir)

#define the file containing a list of files to process
filename <- "input_data/file_list.txt"
file_list <- read.table(filename)
colnames(file_list)[1] <- 'Filename'

#show the file names
cat('The file names we will use are: ')
file_list

#load the single cell experiment with the single cell data based on nuclei segmentation
nuclei_single_cell_filename <-paste(myWorkDir,'/input_data/single_cell_data.rds',sep='')
spe <- readRDS(nuclei_single_cell_filename)

#show how many single cells are available per file name.
#please check that this must check the file names of the masks shown just above.
cat('The number of single cells per image are: ')
table(spe$sample_id)

#initialise colors
cols <- brewer.pal(3, "BuGn")
pal <- colorRampPalette(cols)

#########################################################################  
# Process all samples in bulk based on the filenames in the file list 
#########################################################################  


start_time = Sys.time()
#parse the tiff files one by one and count them 
myC <- 1
#store some plots
nerve_plot_list <- list()
combined_plot_list <- list()
plorIdx <- 1
combIdx <- 1
#save some summary stats about the nerve hotspots found
summaryNerveStats <- data.frame()
for(idx in file_list$Filename){
  
  str <- paste('Processing sample: ',idx,' ',myC,'/',nrow(file_list),sep='')  
  print(str)
  #read the image --> check that the complete input file name is correct
  currentFilename <- paste(myWorkDir,'input_data/',idx,sep='')
  img <- readImage(currentFilename,as.is = TRUE)
  #split the image into the channels
  img1 <- getFrame(img,1)
  img2 <- getFrame(img,2)
  #make watershed segmentation
  #we use a small tolerance value to combine scattered nerve pixels together
  #but this parameter is required to be tuned
  myTolerance <- 0.0002
  nmask <- watershed( distmap(img1),myTolerance,ext=1)
  #retrieve spatial features for the mask
  coords <- computeFeatures.moment(nmask,img2)
  coords <- as.data.frame(coords)
  #add more features
  coords1 <- computeFeatures.shape(nmask,ref = img2)
  coords1 <- as.data.frame(coords1)
  #and more features....
  coords2 <- computeFeatures.basic(nmask,ref = img2)
  coords2 <- as.data.frame(coords2)
  #combine everything into one data frame
  coords$Area <- coords1$s.area
  coords$Radius <- coords1$s.radius.mean
  coords$Intensity <- coords2$b.mean
  
  #filter the nerve pixel coords using the Area value
  #use cutoff lower than 5% quantile --> this can be adjusted as well depending on the application
  myCutOff <- quantile(coords$Area,0.05)
  
  #data frame coords contains the X,Y coordinates of the nerve hotspots together with extra information about the Area and the radius of these. The intensity value also shows how bright were the pixels that correspond to the antibody used for detecting these pixels.
  
  coords <- coords[coords$Area>myCutOff,]
  
  #we require at least 2 to plot
  if(nrow(coords)>2){
    
   o <- ggplot() +
      geom_circle(aes(x0 = m.cx, y0 = -m.cy, r = Radius, fill = Area), data = coords)+
      scale_fill_gradientn("Area",colours = pal(20))+
      ggtitle(idx)+
     ylim(-850,0)+
     xlim(0,850)+
      theme_bw()+
      theme(aspect.ratio = 1)
   nerve_plot_list[[plorIdx]] <- o
   plorIdx <- plorIdx + 1
  }
  
  #store the data to be tabulated
  if(nrow(coords)>0){
    df <- data.frame(Sample=idx,
                     NumberOfPixels=nrow(coords),
                     AvgArea=mean(coords$Area),
                     AvgRadius=mean(coords$Radius),
                     AvgPeripherin=mean(coords$Intensity))
  }else{
    df <- data.frame(Sample=idx,
                     NumberOfPixels=0,
                     AvgArea=0,
                     AvgRadius=0,
                     AvgPeripherin=0)
  }
  summaryNerveStats <- rbind(summaryNerveStats,df)
  
  ##############################################################
  #find the common file names between the spe and the list of tiffs from the nerve segmentation
  #intersect(gsub('.tiff','',file_list$Filename),unique(spe$sample_id))
  
  #add the info about cells and their phenotypes
  #only for the ones that have coordinates
  if(nrow(coords)>0){
      
    currentSample <- gsub('.tiff','',idx)
    #subset the complete set of single cells called spe from before
    tmp <- spe[, spe$sample_id == currentSample]
    matCoord <- spatialCoords(tmp)
    axis_major_lenght <- tmp$axis_major_length
    eccentricity <- tmp$eccentricity
    area <- tmp$area
    #retrieve from the spe the single cell coordinates together with some other info required
    #note that the tmp$CombinedAnnotation column contains the single cell phenotyping results
    #that must be performed in advance before running this pipeline. 
    #Any of the state of the art methods for phenotyping can be applied including unsupervised clustering or gating 
    dt <- data.frame(X=matCoord[,1],
                     Y=matCoord[,2],
                     axis_major_lenght=axis_major_lenght,
                     eccentricity=eccentricity,
                     theta=-1,
                     Area=0.8,
                     Radius=0.8,
                     Intensity=-1,
                     CombinedAnnotation=tmp$CombinedAnnotation)
    #make the column names in coords consistent with the column names of dt from before 
    colnames(coords)[1] <- "X"
    colnames(coords)[2] <- "Y"
    colnames(coords)[3] <- 'axis_major_lenght'
    colnames(coords)[4] <- 'eccentricity'
    colnames(coords)[5] <- 'theta'
    coords$CombinedAnnotation <- 'NerveHotspot'
    #bind them together
    dt <- rbind(dt,coords)
    
    #here for over simpliication we merge some of the specific cell types found in cancer patients into one cell type that we call 'Cancer'
    #please note that this step is application-specific and different applications might have different cell types that cannot be merged into one.
    #if that's the case please skip this step, or modify your cell types accordingly
    dt[dt$CombinedAnnotation=='Epithelial_CK14+'| dt$CombinedAnnotation=='Epithelial_CK14+_CK818+' |
         dt$CombinedAnnotation=='Epithelial_CK56+' | dt$CombinedAnnotation=='Epithelial_CK56+_CK14+' |
         dt$CombinedAnnotation=='Epithelial_CK56+_CK818+' | dt$CombinedAnnotation=='Epithelial_CK818+' |
         dt$CombinedAnnotation=='Epithelial_Other' |dt$CombinedAnnotation=='Undefined'  , 'CombinedAnnotation'] <- 'Cancer'
    
    
    #perform some data manipulation to be able to generate a combined spatial experiment that contains pixel-based nerve hotspots and single-cells
    dt$sample_id <- currentSample
    dt$ROI <- currentSample
    all_coords <- dt[,c("X", "Y")]
    sce <- SpatialExperiment(assays=list(counts=matrix(-1,ncol = nrow(dt),nrow=35)),
                             sample_id = currentSample,
                             #spatialCoordsNames=c('X','Y'),
                             spatialCoords = as.matrix(all_coords))
    
    sce$sample_id <- currentSample
    #this command can be also modified with different parameters, depending on specific application used
    sce <- buildSpatialGraph(sce, img_id = "sample_id", type = "expansion", threshold = 20,
                             coords=c('X','Y'))
    sce$Area <- dt$Area
    sce$Radius <- dt$Radius
    sce$celltype <- dt$CombinedAnnotation
    #again this part is application specific, and must be changed depending on the application and the availabel cell types in the dataset
    sce$celltype <- factor(sce$celltype, levels=c("Cancer","Endothelial","Stromal",
                                                  "B_cells","Macrophages_CD163-","Macrophages_CD163+",
                                                  "T_cytotoxic" ,"T_helper" ,"T_regulatory","Immune_other",
                                                  'NerveHotspot'))
    
    #define some colors, more colors must be defined if more cell types are available in future applications
    mycols_basel_meta <- c('tomato3','orange','thistle2',
                           'mediumpurple1','skyblue2','steelblue3',
                           'springgreen2','springgreen3','springgreen4','darkolivegreen4',
                           'black')
    
    o <- plotSpatial(sce,
                     node_color_by = "celltype",
                     img_id = "sample_id", 
                     draw_edges = TRUE, 
                     colPairName = "expansion_interaction_graph", 
                     nodes_first = FALSE,
                     #node_size_by="Radius",
                     node_size_fix = 0.8,
                     edge_width_fix = 0.05,
                     edge_color_fix = "grey",
                     coords=c('X','Y'))+ 
      scale_color_manual(values=mycols_basel_meta)+
      theme(legend.position = 'bottom',
            axis.text.y = element_blank(),
            axis.text.x = element_blank())+
      guides(colour = guide_legend(override.aes = list(size=0.8),nrow=2,title=""))
    
    combined_plot_list[[combIdx]] <- o
    combIdx <- combIdx + 1
  }
  
  ##############################################################
  #combine spatial experiments after adding the nerve hotspots
  #this is useful for future analysis
  ##############################################################
  if(myC==1){
    #this is the first sample so we just initialise
    allHotspotsSpe <- sce
  }else{
    allHotspotsSpe <- cbind(allHotspotsSpe,sce)
  }
  myC <- myC+1
}
end_time = Sys.time()
end_time-start_time

#show the updated spatial experiment object with the nerve hotspots
cat('The total size of the objects including single-cells and nerve hotspots: ')
table(allHotspotsSpe$sample_id)

#show the summary statistics of nerve hotspots found 
cat('Show a summary of all nerve pixels found: ')
summaryNerveStats

#and visualise in bulk all images generated from the previous code
for(idx in 1:length(combined_plot_list)){
  plot(combined_plot_list[[idx]])
}

for(idx in 1:length(nerve_plot_list)){
  plot(nerve_plot_list[[idx]])
}

```


### Downstream analyses - part I

Now that the pixel-based detection of nerve hotspots has been incorporated to the single-cell experiment, we can perform different types of downstream analysis.

First we will summarise the hotspots per image to engineer features and we will perform some visualisations of the nerve hotspots. 


```{r downstream analysis part I, echo=FALSE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}
library(sp)
#since we have to generate grid with need a step that describes the size of hotspots
myStep <- 25
GridEnd <- 850 - myStep
#the analysis might have to be repreated for different size of grid steps

myC <- 1
nerve_hotspot_list <- list()
allHotSpots <- data.frame()
for(idx in file_list$Filename){
  
  str <- paste('Processing sample: ',idx,' ',myC,'/',nrow(file_list),sep='')  
  print(str)
  currentSample <- gsub('.tiff','',idx)
  tmpSpe <- allHotspotsSpe[,allHotspotsSpe$sample_id==currentSample]
  #focus only on the NerveHotspot
  tmpSpe <- tmpSpe[,tmpSpe$celltype=='NerveHotspot']
  #retrieve the coordinates
  tmpCoord <- spatialCoords(tmpSpe)
  dt <- data.frame(Area=tmpSpe$Area,
                   Radius=tmpSpe$Radius,
                   tmpCoord)

  #generate grid coordinates with function
  #griddf <- expand.grid( Y = seq(from = -850, to=0,by = myStep),
  #                     X = seq(from = 0, to=850,by = myStep))

  #generate grid squares and find overlap with points from nerve hotspots
  HotSpotSummary <- data.frame()
  for(myX in seq(from = 0, to=GridEnd,by = myStep)){
    X_start <- myX
    X_end <- myX+myStep
    for(myY in seq(from = 0, to=GridEnd,by = myStep)){
      Y_start <- myY
      Y_end <- myY+myStep
      #at this point we have the coordinates and we have to screen our nerve pixels coordinates to find if they are inside the square
      tmpHotspot <- data.frame()
      for(K in 1:nrow(dt)){
        X_nerve <- dt[K,'X']
        Y_nerve <- dt[K,'Y'] 
        #very careful here on how you write the coordinates
        res <- point.in.polygon(c(X_nerve),c(Y_nerve),
                                c(X_start,X_end,X_end,X_start),
                                c(Y_start,Y_start,Y_end,Y_end))
        if(res==0){
          radius <- 0
          area <- 0
          inside <- 0
        }else{
          radius <- dt[K,'Radius']
          area <- dt[K,'Area']
          inside <- 1
        }
        tmpRes <- data.frame(X=X_start+(myStep/2),
                                 Y=Y_start+(myStep/2),
                                 Sx1=X_start,
                                 Sx2=X_end,
                                 Sy1=Y_start,
                                 Sy2=Y_end,
                                 X_nerve=X_nerve,
                                 Y_nerve=Y_nerve,
                                 Inside=inside,
                                 Radius=radius,
                                 Area=area)
        tmpHotspot <- rbind(tmpHotspot,tmpRes)
      }
      
      #aggregate and store
      if(sum(tmpHotspot$Inside)==0){
        #this means that there is no pixel inside the square so we add zero to everything
        tmpRes <- data.frame(X=X_start+(myStep/2),
                                 Y=Y_start+(myStep/2),
                                 Sx1=X_start,
                                 Sx2=X_end,
                                 Sy1=Y_start,
                                 Sy2=Y_end,
                                 X_nerve=-1,
                                 Y_nerve=-1,
                                 AvgDist=-1,
                                 PixelsFound=0,
                                 AvgRadius=0,
                                 AvgArea=0,
                                 MaxRadius=0,
                                 MaxArea=0)
      }else{
        #here means that we have at least one success
        a <- which(tmpHotspot$Inside==1)
        tmpHotspot <- tmpHotspot[a,]
        tmp <- tmpHotspot[,c(7,8)]
        if(nrow(tmpHotspot)==1){
          myDist <- 0
        }else{
          myDist <- mean(dist(tmp),na.rm=T)
        }
        tmpRes <- data.frame(X=X_start+(myStep/2),
                                 Y=Y_start+(myStep/2),
                                 Sx1=X_start,
                                 Sx2=X_end,
                                 Sy1=Y_start,
                                 Sy2=Y_end,
                                 X_nerve=mean(tmpHotspot$X_nerve),
                                 Y_nerve=mean(tmpHotspot$Y_nerve),
                                 AvgDist=myDist,
                                 PixelsFound=nrow(tmpHotspot),
                                 AvgRadius=mean(tmpHotspot$Radius,na.rm=T),
                                 AvgArea=mean(tmpHotspot$Area,na.rm=T),
                                 MaxRadius=max(tmpHotspot$Radius,na.rm=T),
                                 MaxArea=max(tmpHotspot$Area,na.rm=T))
        
      }
      HotSpotSummary <- rbind(HotSpotSummary,tmpRes)
    }
  }
  
  # plot heatmap of expression for 100 clusters
breaks <- seq(0, 1, by = 0.05)
white.red <- colorRampPalette(c("white", "red"))(n = 20)
  
  o <- ggplot()+
    geom_raster(data=HotSpotSummary,aes(x=X,y=-Y,fill=AvgArea),size=0.3,shape=4)+
    scale_fill_gradientn("Nerve hotspots (avg.Area) ",colours = white.red)+
    geom_point(data=HotSpotSummary,aes(x=X,y=-Y),size=0.4,color='gray88')+
    theme_minimal()+
    ggtitle(idx)+
    theme(axis.text.y = element_text(size = 8 ),
          axis.text.x = element_text(size = 8,angle = 0, vjust = 0.5, hjust = 0.5),
          axis.title.x = element_text( size = 8),
          axis.title.y = element_text(size = 8),
          strip.text = element_text(size = 8,face='bold',lineheight=1),
          legend.position = "bottom",aspect.ratio = 1)

  nerve_hotspot_list[[myC]] <- o 
  myC <- myC + 1
  
  #gather info from all samples
    HotSpotSummary$Sample <- currentSample
    HotSpotSummary$TotalPixels <- nrow(dt)
    HotSpotSummary$AvgPixelArea <- mean(dt$Area,na.rm=T)
    HotSpotSummary$AvgPixelRadius <- mean(dt$Radius,na.rm=T)
    allHotSpots <- rbind(allHotSpots,HotSpotSummary)

}

for(idx in 1:length(nerve_hotspot_list)){
  plot(nerve_hotspot_list[[idx]])
}

#show the summary statistics of nerve hotspots found 
cat('Show info about nerve hotspots: ')
head(allHotSpots)


```



### Downstream analyses - part II

With the objects we have created in the previous subsection the spatial analysis pipelines published here https://bodenmillergroup.github.io/IMCDataAnalysis/performing-spatial-analysis.html are readily applicable and the results of this analysis can be associated with clinical info if available.

We continue the analysis by computing the cellular abundance per hotspot. 

With this we will be able to perform differential abundance analysis and discover potential associations with nerve elements.

For example, the code introduced below can be fed to the diffcyt pipeline to find differences in the abundance between areas (squares) that have nerve elements, versus the ones that do not. 

TODO: more applications and scenarios will be presented in the future

```{r downstream analysis part II, echo=FALSE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}


library(sp)
#since we have to generate grid with need a step that describes the size of hotspots
myStep <- 25
GridEnd <- 850 - myStep
#the analysis might have to be repreated for different size of grid steps
myC <- 1
myFilenames <- unique(allHotspotsSpe$sample_id)
HotSpotAbundance <- data.frame()

start_time = Sys.time()
for(idx in myFilenames){
  
  str <- paste('Processing sample: ',idx,' ',myC,'/',length(myFilenames),sep='')  
  print(str)
  #currentSample <- gsub('.tiff','',idx)
  currentSample <- idx
  tmpSpe <- allHotspotsSpe[,allHotspotsSpe$sample_id==currentSample]
  #retrieve the coordinates of all elements in the image including the cell type
  tmpCoord <- spatialCoords(tmpSpe)
  dt <- data.frame(Area=tmpSpe$Area,
                   Radius=tmpSpe$Radius,
                   CellType=tmpSpe$celltype,
                   tmpCoord)

  
  #generate grid squares and find overlap with points from nerve hotspots
  for(myX in seq(from = 0, to=GridEnd,by = myStep)){
    X_start <- myX
    X_end <- myX+myStep
    for(myY in seq(from = 0, to=GridEnd,by = myStep)){
      Y_start <- myY
      Y_end <- myY+myStep
      
      tmpHotspot <- data.frame()
      res <- point.in.polygon(c(dt$X),c(dt$Y),
                                c(X_start,X_end,X_end,X_start),
                                c(Y_start,Y_start,Y_end,Y_end))
      a <- which(res!=0)
      tmpHotspot <- dt[a,]
      tmpHotspot$CellType <- factor(tmpHotspot$CellType,levels=c("Cancer","Endothelial","Stromal",
                                                         "B_cells","Macrophages_CD163-","Macrophages_CD163+",
                                                            "T_cytotoxic" ,"T_helper" ,"T_regulatory","Immune_other",
                                                            'NerveHotspot'))
      tmpAbund <- as.data.frame(table(tmpHotspot$CellType))
      n_cells <- sum(tmpAbund$Freq)
      n_nerves <- length(which(tmpHotspot$CellType=='NerveHotspot'))
      #tmpAbund$Freq <- tmpAbund$Freq/sum(tmpAbund$Freq)
      tmpAbund <- tmpAbund[order(tmpAbund$Var1, decreasing = F),]
      tmpAbund <- as.data.frame(t(tmpAbund))
      myColNames <- tmpAbund[1,]
      tmpAbund <- tmpAbund[2,]
      colnames(tmpAbund) <- myColNames
      
      tmpAbund$NCells <- n_cells
      tmpAbund$HasNerve <- ifelse(n_nerves>0,1,0)
      tmpAbund$X <- X_start+(myStep/2)
      tmpAbund$Y <- Y_start+(myStep/2)
      tmpAbund$Sample <- currentSample
      rownames(tmpAbund) <- NULL
      HotSpotAbundance <- rbind(HotSpotAbundance,tmpAbund)
    }
  }
}

end_time = Sys.time()
end_time-start_time

```


### Downstream analyses - part III

Here we will perform the "classic" neighborhood analysis presented in the pipeline published here https://bodenmillergroup.github.io/IMCDataAnalysis/performing-spatial-analysis.html 

The code it also implements some visualisations about the interactions found between cell types and nerve pixels.

Remember that if a pair of cell types is significantly interacting we have sigval = 1, if a pair of cell types is significantly avoiding we have sigval = -1 and if no significant interaction or avoidance was detected we have sigval = 0.

```{r downstream analysis part III, echo=FALSE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}

library(scales)
library(BiocParallel)

 mycols_basel_meta <- c('tomato3','orange','thistle2',
                           'mediumpurple1','skyblue2','steelblue3',
                           'springgreen2','springgreen3','springgreen4','darkolivegreen4',
                           'black')

interaction_spatial_plots <- list()
interaction_heatmaps <- list()
outAll <- data.frame()
myC <- 1 

start_time = Sys.time()

for(idx in myFilenames){
  
  str <- paste('Processing sample: ',idx,' ',myC,'/',length(myFilenames),sep='')  
  print(str)
  currentSample <- idx
  tmpSpe <- allHotspotsSpe[,allHotspotsSpe$sample_id==currentSample]
  #retrieve the coordinates of all elements in the image including the cell type
  tmpCoord <- spatialCoords(tmpSpe)
  sce <- SpatialExperiment(assays=list(counts=matrix(-1,ncol = nrow(tmpCoord),nrow=35)),
                             sample_id = idx,
                             #spatialCoordsNames=c('X','Y'),
                             spatialCoords = as.matrix(tmpCoord))
    
    sce$sample_id <- idx
    sce$celltype=tmpSpe$celltype
  
  
    sce <- buildSpatialGraph(sce, img_id = "sample_id", type = "knn", k = 20,
                             coords=c('X','Y'))
    sce <- buildSpatialGraph(sce, img_id = "sample_id", type = "expansion", threshold = 20,
                             coords=c('X','Y'))
    sce <- buildSpatialGraph(sce, img_id = "sample_id", type = "delaunay", max_dist = 20,
                             coords=c('X','Y'))
    
    #find the interactions per sample 
    out <- testInteractions(sce, 
                          group_by = "sample_id",
                          label = "celltype", 
                          colPairName = "knn_interaction_graph",
                          BPPARAM = SerialParam(RNGseed = 221029))
    out <- as.data.frame(out)
    #visualise 
    o1 <- out %>% as_tibble() %>%
         group_by(from_label, to_label) %>%
        summarize(sum_sigval = sum(sigval, na.rm = TRUE)) %>%
        ggplot() +
        geom_tile(aes(from_label, to_label, fill = sum_sigval)) +
        ggtitle(idx)+
        scale_fill_gradient2(low = muted("blue"), mid = "white", high = muted("red")) +
        theme(axis.text.x = element_text(angle = 45, hjust = 1))
    
    #merge all interaction outs
    outAll <- rbind(outAll,out)
    
    sce$celltype <- factor(sce$celltype, levels=c("Cancer","Endothelial","Stromal",
                                                         "B_cells","Macrophages_CD163-","Macrophages_CD163+",
                                                            "T_cytotoxic" ,"T_helper" ,"T_regulatory","Immune_other",
                                                            'NerveHotspot'))
      o2 <- plotSpatial(sce,
                  node_color_by = "celltype",
                  img_id = "sample_id", 
                  draw_edges = TRUE, 
                  colPairName = "knn_interaction_graph", 
                  nodes_first = FALSE,
                  node_size_fix = 0.8,
                  edge_width_fix = 0.05,
                  edge_color_fix = "grey",
                  coords=c('X','Y'))+ 
        scale_color_manual(values=mycols_basel_meta)+
        theme(legend.position = 'bottom',
              axis.text.y = element_blank(),
              axis.text.x = element_blank())+
        guides(colour = guide_legend(override.aes = list(size=0.8),nrow=2,title=""))
    
    interaction_spatial_plots[[myC]] <- o2  
    interaction_heatmaps[[myC]] <- o1
    
    myC <- myC + 1
}
  
end_time = Sys.time()
end_time-start_time

#visualise the images with interactions 
for(idx in 1:length(interaction_spatial_plots)){
  plot(interaction_spatial_plots[[idx]])
}

#and the heatmaps per sample
#for(idx in 1:length(interaction_heatmaps)){
#  plot(interaction_heatmaps[[idx]])
#}

########################################################################
#visualise the interactions for the full cohort
########################################################################

    outAll <- as.data.frame(outAll)
    
    #visualise 
    outAll %>% as_tibble() %>%
         group_by(from_label, to_label) %>%
        summarize(sum_sigval = sum(sigval, na.rm = TRUE)) %>%
        ggplot() +
        geom_tile(aes(from_label, to_label, fill = sum_sigval)) +
        ggtitle(idx)+
        scale_fill_gradient2(low = muted("blue"), mid = "white", high = muted("red")) +
        theme(axis.text.x = element_text(angle = 45, hjust = 1))


```