exp3_cellTypeSignatures.Rmd

---
title: "exp3_pulpDeconvolution"
output: html_document
author: Sara Gosline
date: "2023-01-27"
---

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

library(ggplot2)
library(ggfortify)
library(cowplot)
#library(leapR)
library(dplyr)

source('spleenDataFormatting.R')
library(spammer)
#source('spatialProtUtils.R')
```

The assignment of pulp type from the basic data was not robust. Here we will use the cell type signatures from the sorted data to label the voxels.


## Cell type signatures

Here we get the differential expression from the sorted cells.

```{r cell type signatures}


##get fit values for diffex
pumap<-scater::runPCA(spat.prot)

phumap<-spat.phos%>%
  scater::runPCA()

fullmap<-scater::runPCA(global.sorted)

fullmap<-spatialDiffEx(fullmap)

full<- fullmap%>%
  rowData(.)%>%
  as.data.frame()%>%
  dplyr::select(featureID='X',
                logFC='pulpAnnotation.limma.logFC',
                adj.P.Val='pulpAnnotation.limma.adj.P.Val',
                AveExpr='pulpAnnotation.limma.AveExpr')

upsig<-full%>%
  subset(adj.P.Val<0.01)%>%
  subset(logFC>1)

downsig<-full%>%
  subset(adj.P.Val<0.01)%>%
  subset(logFC<(-1))

print(paste0('we now have a white pulp signature of ',nrow(upsig)))
print(paste0('we now have a red pulp signature of ',nrow(downsig)))

```

## Plot values in heatmap
Let's now plot the signatures in heatmaps of the original and spatial data.


```{r pressure, echo=FALSE}
library(pheatmap)


pheatmap(as.matrix(exprs(fullmap)[rownames(upsig),]),annotation_col = as.data.frame(colData(fullmap)),cellwidth = 10,main='White pulp signature')


pheatmap(as.matrix(exprs(pumap)[intersect(rownames(upsig),rownames(rowData(pumap))),]),annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],main='White pulp signature')


pheatmap(as.matrix(exprs(fullmap)[rownames(downsig),]),annotation_col = as.data.frame(colData(fullmap)),cellwidth=10,main='Red pulp signature')


pheatmap(as.matrix(exprs(pumap)[intersect(rownames(downsig),rownames(rowData(pumap))),]),annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],main='Red pulp signature')


```

These signatures cluster voxels, but don't serve as great markers because tehy are still very much 'down' regulated. We want to pick the top expressed proteins that represent each pulp type.

```{r reduced signature}


whitesig<-full%>%
  subset(adj.P.Val<0.01)%>%
  subset(logFC>1)%>%
  subset(AveExpr>1)

redsig<-full%>%
  subset(adj.P.Val<0.01)%>%
  subset(logFC<(-1))%>%
  subset(AveExpr>(1))

print(paste0('we now have a reduced white pulp signature of ',nrow(whitesig)))
print(paste0('we now have a reduced red pulp signature of ',nrow(redsig)))

pheatmap(as.matrix(exprs(fullmap)[rownames(whitesig),]),annotation_col = as.data.frame(colData(fullmap)),cellwidth = 10,
         main='White pulp signature',filename='whitePulpsorted.png')


pheatmap(as.matrix(exprs(pumap)[intersect(rownames(whitesig),rownames(rowData(pumap))),]),annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],
         main='White pulp signature',filename='whitePulp2d.png')


pheatmap(as.matrix(exprs(fullmap)[rownames(redsig),]),annotation_col = as.data.frame(colData(fullmap)),cellwidth=10,main='Red pulp signature',
         filename='redPulpsorted.png')


pheatmap(as.matrix(exprs(pumap)[intersect(rownames(redsig),rownames(rowData(pumap))),]),annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],main='Red pulp signature',
         filename='redPulp2d.png')

```

So we have a handful of proteins that are marking the signature. Let's see if we can link them in a biological fashion.

```{r pathway enrichment}
#library(clusterProfiler)
library(leapR)
library(org.Hs.eg.db)
data('krbpaths')
map<-as.list(org.Hs.egSYMBOL2EG)

redgenes<-lapply(rownames(redsig),function(x) return(map[[x]][1]))
whitegenes<-lapply(rownames(whitesig),function(x) return(map[[x]][1]))

universe<-unlist(lapply(rownames(full),function(x) return(map[[x]][1])))



order.enrich<-leapR::leapR(geneset=krbpaths,
                         enrichment_method='enrichment_in_order',id_column='featureID',
                    datamatrix=full,primary_columns='logFC')#,greaterthan=T,threshold=1)
sig.enrich<-order.enrich%>%
  subset(BH_pvalue<0.01)

print(sig.enrich)

groups<-stringr::str_split(sig.enrich[1,]$ingroupnames,pattern=', ')%>%unlist()


pheatmap(as.matrix(exprs(fullmap)[intersect(groups,rownames(rowData(fullmap))),]),
         annotation_col = as.data.frame(colData(fullmap)),
         main='CD3 and TCR signature')

pheatmap(as.matrix(exprs(pumap)[intersect(groups,rownames(rowData(pumap))),]),
         annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],
         main='CD3 and TCR signature')

pheatmap(as.matrix(exprs(fullmap)[intersect(groups,rownames(rowData(fullmap))),]),
         annotation_col = as.data.frame(colData(fullmap)),
         main='CD3 and TCR signature',filename='cd3tcrSorted.png')

pheatmap(as.matrix(exprs(pumap)[intersect(groups,rownames(rowData(pumap))),]),
         annotation_col = as.data.frame(colData(pumap))[,c('pulpAnnotation','TMT.set')],
         main='CD3 and TCR signature',filename='cd3tcr2d.png')

```

There is only one kegg/reactome pathway that is enriched, and it is related to TCR signaling. 

```{r calculate signature scores}

pumap<-calcSigScore(pumap,rownames(redsig),'RedPulp')%>%
  calcSigScore(rownames(whitesig),'WhitePulp')
  
p1<-plotPCA(pumap,color_by='WhitePulp')

p2<-plotPCA(pumap,color_by='RedPulp')

p3<-cowplot::plot_grid(p2,p1)
p3

ggsave('pcaLabeledByScore.png',p3,width=12)

```

This isn't quite right - the 'green' items are definitely on the right, but we need some sort of threshold to 'call' the white vs. red pulp items.


### Plot in image

Lastly we want to plot a set of scores, in this case the signature scores, in an image. 


```{r plot image}


 p <- plotGrid(pumap,'RedPulp')
 p1<-plotGrid(pumap,'WhitePulp')
    

 pc<-cowplot::plot_grid(p,p1)
 pc
 ggsave('scoredImage.png',pc,width=12)
```

Now that we have a rank score, we can "call" each voxel. 


```{r voxel scoring}

cols<-list(red='darkred',white='white',None='darkgrey')
newDat<-colData(pumap)%>%
  as.data.frame()%>%
  mutate(isRed=RedPulp>0.5,isWhite=WhitePulp>0.5)%>%
  mutate(pulp=ifelse(isRed,'red',ifelse(isWhite,'white','None')))

p1<-ggplot(newDat)+geom_point(aes(x=RedPulp,y=WhitePulp,col=pulp))+scale_color_manual(values=cols)

p2<-ggplot(newDat,aes(x=Xcoord,y=Ycoord,fill=pulp))+geom_raster()+scale_fill_manual(values=cols)

res<-cowplot::plot_grid(p1,p2)
res
ggsave('assignedVals.png',res,width=12)
```