Fig1_additional_unused_panels.Rmd

---
title: "Misc plots-Unused Fig 1"
output: html_document
date: "2024-08-03"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.height = 10,fig.width = 7,include = FALSE)
#### sets tab autocompletion for directories to begin from the directory containing the .Rproj session, instead of the script's directory if different
knitr::opts_knit$set(root.dir = here::here())


library(data.table) # Preferred data manipulation package
library(ggplot2) # Dependency for several plotting functions
library(ggtext) # more character types in ggplots
library(ggrastr) # avoid making raster images one can't even open
library(SpatialExperiment) # Visium data framework
library(SpatialFeatureExperiment) # Xenium data framework
library(spatialLIBD) # one option for plotting cluster assignments, but breaks when the in-tissue portion of a visium area is very un-square.
library(escheR) # alternative spotplotting function, at least for visium
library(viridis) # palettes
library(Polychrome) # better palettes
require(colorout) # Utility for RStudio
library(sf) # define the polygonal boundaries of xenium domains
ColorOut()

# code reformatting in Rstudio
options("styler.addins_style_transformer" = "biocthis::bioc_style()")

# set plotting defaults for ggplot
theme_set(theme_bw() + theme(axis.text.x = element_text(size = 8), axis.title.x = element_text(size = 9), axis.text.y = element_text(size = 8), axis.title.y = element_text(size = 9), plot.title = element_blank(), strip.text = element_text(size = 10), legend.text = element_text(size = 8), legend.title = element_text(size = 8, hjust = 0.5)))
```

## load visium, xenium experiments, cluster labels
```{r}
#### VISIUM ####
hyp2 <- readRDS("spatial_HYP/processed-data/03-QC_filters/hypN10_umi210_gene126_chrm50_spotsweeped_lognorm_rotsNmirrors_072224.RDS")

bscl <- fread("spatial_HYP/processed-data/06-BayesSpace/01-bayesspace60kiter_k15-20-31_out/BSpace_k15_HARMONYlmbna_nnsvg10.txt") # main visium clustering (bs=bayesspace)

setnames(bscl,2,"cl")
bscl[,cl:=paste0("Vis",cl)]
bscl[cl=="Vis7",cl:="VMH.1"]
bscl[cl=="Vis12",cl:="VMH.2"]
bscl[cl=="Vis4",cl:="ARC.1"]
bscl[cl=="Vis6",cl:="ARC.2"]
# WM (optic tract, OT): Clusters 3 and 10
bscl[cl=="Vis3",cl:="OT.1"]
bscl[cl=="Vis10",cl:="OT.2"]

bscl[,dom:=cl]
bscl[dom %in% c("VMH.1","VMH.2"),dom:="VMH"]
bscl[dom %in% c("ARC.1","ARC.2"),dom:="ARC"]
bscl[dom %in% c("OT.1","OT.2"),dom:="OT"]
bscl[!(dom %in% c("VMH","ARC","OT")),dom:="Other"]

bscl<-DataFrame(bscl,row.names=bscl$rn)[colnames(hyp2),]
hyp2$k15dom <- bscl$cl
hyp2$`Visium Domain` <- factor(bscl$dom,levels=c("ARC","VMH","OT","Other"))

#### XENIUM ####
hypx <- readRDS("xenium_HYP/processed-data/05_Banksy_M0lam0_res2_multisamp/01-sfe-in_staggered-coords_genetarg-only_NONlog-norm.RDS")

bksmooth <- fread("xenium_HYP/processed-data/07_Domains-subdomains_by_knnSmoothing_Banksytypes/03b-ARCVMHdomains_2stepsmooth_VMH-k10-0.2-VMH-k200-0.2_ARC-k50-0.1_ARC-k500-0.5.txt") ## xenium domain assignments after smoothing, by xenium cell

bksmooth[dualVMHARC4=="other",dualVMHARC4:="Other"]
bksmooth <- DataFrame(bksmooth,row.names=bksmooth$rn)[colnames(hypx),]
hypx$dom <- bksmooth$dualVMHARC4
```

### previously made a working 15-color palette, so we can use this for both Xen/Vis VMH/ARC at subdomain/cluster levels (max 8 VM+ARC combined, in xenium) and for overarching domains (2 colors + a gray for other)

that prior palette, with slots 1-2 for ARC, 3-4 for VMH, and 5-6 for WM or other, was: "#269422", "#169482", "#4a2881", "#6A0D53", "#3D4047", "#A9A7A9", "#F70022", "#551CFD", "#F59B26", "#00B3FD", "#FF0DF7", "#D193FE", "#D1CA1C", "#FB479F", "#AD260D"
The colors slotted in 1 and 3 are also picked to have the same luminance so that gene expression coloration fill in escheR grayscale doesn't look falsely darker or lighter in one spot type vs. the other 
```{r}
pal <- c("#8c63cf","#5dd959","#F59B26","#A9A7A9")
names(pal) <- c("VMH","ARC","OT","Other")

### also save the prior palette for quick loadup
fullpal <- c("#8c63cf", "#169482", ,"#5dd959" "#6A0D53", "#3D4047", "#A9A7A9", "#F70022", "#551CFD", "#F59B26", "#00B3FD", "#FF0DF7", "#D193FE", "#D1CA1C", "#FB479F", "#AD260D")
saveRDS(fullpal,"manuscript_plots/15color_palette.RDS")
rm(fullpal)
####
```

(GABRE marking Visium ARC )
```{r}
panh <- hyp2["GABRE",hyp2$sample_id=="V12D05-350_C1"]
colData(panh)$`log counts` <- as.numeric(logcounts(panh)["GABRE",])

### we need to get the hulls for the VMH and ARC domains, which I thought we were done with but hey


p <- make_escheR(panh)
p <- p |> add_fill("log counts",size=0.25,point_size = 0.5)
p <- p |> add_ground(var = "Visium Domain",stroke=0.15,point_size = 0.5)

pdf("manuscript_plots/Fig1/1H-V12D05_350_C1_GABRE.pdf",height=2,width=2)
p+
    
    scale_color_manual(values=pal,na.value = NA)+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    #scale_fill_grey()+
    ggtitle("V13Y24-346_C1 (Visium)<br>GABRE Expression")+
    guides(color=guide_legend(override.aes=list(size=1.5,stroke=1)))+
    guides(fill=guide_colorbar(order=1),color=guide_legend(order=2))+
    labs(color="Visium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin = margin(-0.05,0,0,0,unit="in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0,0.05,0,-0.05,unit = "in"))
dev.off()

rm(p,panh)
```

GABRE marking Xenium ARC

Panel J: GABRE in Xenium X99_8741C
```{r}
panj <- hypx["GABRE",hypx$sample_id=="X99_8741C"]

colData(panj)$`log counts` <- as.numeric(logcounts(panj)["GABRE",])

p <- make_escheR(panj,y_reverse=FALSE)
p <- p |> add_fill("log counts",size=0.125,point_size = 0.125)

pdf("manuscript_plots/Fig1/1J-X99_8741C_GABRE.pdf",height=2,width=2)
p+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    ggtitle("X99_8741C (Xenium)<br>GABRE Expression")+
    geom_path(data=dompoly,aes(x=xpol,y=ypol,group=dompol,col=dompol),linewidth = 0.25)+
    labs(color="Xenium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin = margin(0,0,-0.15,0,unit="in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0.05,0.05,-0.1,-0.05,unit = "in"))
dev.off()

## keep the polygons for xenium domain borders -- we're showing this sample for the VMH markers too. 
```

visium domains in V13Y24-346_C1 (Br1225)
```{r}
pane <- hyp2[,hyp2$sample_id=="V13Y24-346_C1"]

p <- make_escheR(pane)
p <- p |> add_fill("Visium Domain",size=0.25,point_size = 0.6)

pdf("manuscript_plots/Fig1/1E-V13Y24_346_C1_Visdomains.pdf",height=1.57,width=2)
p+
    scale_fill_manual(values=pal,na.value = NA)+
    ggtitle("V13Y24-346_C1 (Visium)")+
    guides(color="none",fill=guide_legend(override.aes=list(size=1.5)))+
    labs(fill="Visium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.125,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(-0.05,0,-0.05,0,unit = "in"))
dev.off()

rm(pane,p)
```


xenium domains in (X99_1225A)
```{r}
panf <- hypx[,hypx$sample_id=="X99_1225A"]

p <- make_escheR(panf)
p <- p |> add_fill("dom",size=0.05,point_size = 0.2)

pdf("manuscript_plots/Fig1/1F-X99_1225A_Xendomains.pdf",height=1.57,width=2)
p+
    # the reverse_y argument didn't work here, but manually flipping the sign on the y variable designated in p (`.y`) restores the orientation we intended.
    aes(y=-.y)+
    scale_fill_manual(values=pal,na.value = NA)+
    ggtitle("X99_1225A (Xenium)")+
    guides(color="none",fill=guide_legend(override.aes=list(size=1.5)))+
    labs(fill="Xenium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.125,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(-0.05,0,-0.05,0,unit = "in"))
dev.off()

rm(panf,p)

```

FEZF1 marking VMH in xenium

```{r}
pann <- hypx["FEZF1",hypx$sample_id=="X99_8741C"]

colData(pann)$`log counts` <- as.numeric(logcounts(pann)["FEZF1",])

p <- make_escheR(pann,y_reverse=FALSE)
p <- p |> add_fill("log counts",size=0.125,point_size = 0.125)

pdf("manuscript_plots/Fig1/1N-X99_8741C_FEZF1.pdf",height=2,width=2)
p+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    ggtitle("X99_8741C (Xenium)<br>FEZF1 Expression")+
    geom_path(data=dompoly,aes(x=xpol,y=ypol,group=dompol,col=dompol),linewidth = 0.25)+
    labs(color="Xenium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin = margin(0,0,-0.15,0,unit="in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0.05,0.05,-0.1,-0.05,unit = "in"))
dev.off()

rm(pann,dompoly,p)
```


Now for K, L: VMH marker genes in V12D05_350_C1
```{r}
pank <- hyp2["NR5A1",hyp2$sample_id=="V12D05-350_C1"]
colData(pank)$`log counts` <- as.numeric(logcounts(pank)["NR5A1",])

p <- make_escheR(pank)
p <- p |> add_fill("log counts",size=0.25,point_size = 0.5)
p <- p |> add_ground(var = "Visium Domain",stroke=0.15,point_size = 0.5)

pdf("manuscript_plots/Fig1/1K-V12D05_350_C1_NR5A1.pdf",height=2,width=2)
p+
    scale_color_manual(values=pal,na.value = NA)+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    ggtitle("V12D05-350_C1 (Visium)<br>NR5A1 Expression")+
    guides(color=guide_legend(override.aes=list(size=1.5,stroke=1)))+
    labs(color="Visium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin = margin(-0.05,0,0,0,unit="in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0,0.05,0,-0.05,unit = "in"))
dev.off()

rm(p,pank)
```


Panel E: 
Now the VMH markers in Xenium:
```{r}
panm <- hypx["NR5A1",hypx$sample_id=="X99_8741C"]
colData(panm)$`log counts` <- as.numeric(logcounts(panm)["NR5A1",])

p <- make_escheR(panm,y_reverse=FALSE)
p <- p |> add_fill("log counts",size=0.125,point_size = 0.125)

pdf("manuscript_plots/Fig1/1M-X99_8741C_NR5A1.pdf",height=2,width=2)
p+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    ggtitle("X99_8741C (Xenium)<br>NR5A1 Expression")+
    geom_path(data=dompoly,aes(x=xpol,y=ypol,group=dompol,col=dompol),linewidth = 0.25)+
    labs(color="Xenium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin = margin(0,0,-0.15,0,unit="in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0.05,0.05,-0.1,-0.05,unit = "in"))
dev.off()
```


### Set up the xenium palette for cell groups for the next panel.
```{r}
rastumap <- fread("xenium_HYP/processed-data/04_general_QC_and_normalization/02a-genetarg-lctUMAP_rawctCAandUMAP_w_metadata.txt.gz")



## create supercluster labels to group related cell types together
bkcl <- as.data.table(bkcl,keep.rownames=T)
bkcl[,supercl:="ifthisisplottedigoofed"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="VMH",value=T),supercl:="VMH (4)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="ARC",value=T),supercl:="ARC (5)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Tanyc|Alar",value=T),supercl:="Tanycytes, Portal Vasc. (4)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Oligo",value=T),supercl:="Oligo (4)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Astro",value=T),supercl:="Astro (2)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Microg",value=T),supercl:="Microglia (3)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Supraopt",value=T),supercl:="SON (2)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Unsure_AVP",value=T),supercl:="Non-SON AVP+OXT+ (1)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Macrop|Periph|Vascul|Endothel",value=T),supercl:="Vascular and Peripheral Immune (5)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="GABA",value=T),supercl:="Other GABAergic Neurons (2)"]
bkcl[bjm_annot %in% grep(bjm_annot,pattern="Periven",value=T),supercl:="PeriVN (Excitatory) (1)"]
### the parenthetical numbers add up to 33, so we're good.

# key is a field from the colData slot of the SPE that was exported, and is identical to the rownames in case they are lost during switches between d.t. and d.f elsewhere
rastumap <- rastumap[key %in% bkcl$rn]
rastumap <- merge.data.table(rastumap,bkcl,by.x="key",by.y="rn")




### lastly, in the case of the xenium spotplot of all cell groups, we need to put them in a specific order for them to fit within the plot space alongside 1225B: VMH, ARC, SON, Astro, Oligo, Microglia, PeriVN, Non-SON AVP+OXT+, Tanycytes n portal, Other GABA, Vasc n periph
bkcl$supercl <- factor(bkcl$supercl,levels=c("VMH (4)","ARC (5)","SON (2)","Astro (2)","Oligo (4)","Microglia (3)","PeriVN (Excitatory) (1)","Non-SON AVP+OXT+ (1)","Tanycytes, Portal Vasc. (4)","Other GABAergic Neurons (2)","Vascular and Peripheral Immune (5)"))
```


### Panel D is possibly the xenium cluster UMAP, which we should re-generate here with matched palette components where possible.
```{r}
xenpal <- pals[["xen_cellgroups"]]
stopifnot(all(unique(bkcl$supercl) %in% names(xenpal)))

p <- ggplot(rastumap,aes(x=UMAP1_CA,y=UMAP2_CA,color=supercl))+ 
   xlab("UMAP dim 1")+
   ylab("UMAP dim 2")+
   labs(color="Cell group (N clusters)")+
   scale_color_manual(values=xenpal)+
   guides(color=guide_legend(ncol=2,override.aes=list(size=1.5)))+
   theme(legend.position = "bottom",legend.direction="vertical",title = element_blank(),axis.title = element_text(size=32),axis.text=element_blank(),legend.title = element_text(size=24),legend.text = element_text(size=20),legend.background = element_blank(),legend.box.margin = margin(0,0.25,0,-0.5,unit = "in"))+
   rasterize(geom_point(size=0.1),dpi=450,dev = "cairo_png")


pdf("manuscript_plots/Fig1/1D-Xenium_cellgroup_UMAP.pdf",width=9,height=9)
p
dev.off()
```


Panel were not including for now is Visium expression of a WM marker, KLK6. (Our options for those are somewhat limited on the xenium side).
```{r}
## for escheR, we put the assay value into the colData for it to use
panh <- hyp2["KLK6",hyp2$sample_id=="V13Y24-346_C1"]
colData(panh)$`log counts` <- as.numeric(logcounts(panh)["KLK6",])

p <- make_escheR(panh)
p <- p |> add_fill("log counts",size=0.25,point_size = 0.5)
p <- p |> add_ground(var = "Visium Domain",stroke=0.15,point_size = 0.5)

#pdf("manuscript_plots/Fig1/1H-V12D05_350_C1_KLK6.pdf",height=2,width=2)
p+
    scale_color_manual(values=pal,na.value = NA)+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    #scale_fill_grey()+
    ggtitle("*KLK6* ",mscriptids[sample_id=="V13Y24-346_C1",manuscript_id])+
    guides(color=guide_legend(override.aes=list(size=1.5,stroke=1)))+
    labs(color="Visium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.075,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(0,0.05,0,-0.05,unit = "in"),legend.margin = margin(0,0,0,-0.14,"in"))
dev.off()

rm(p,panh)
```


1L: KLK6 in xenium X99_1225B
```{r}
panl <- hypx["KLK6",hypx$sample_id=="X99_1225B"]
colData(panl)$`log counts` <- as.numeric(logcounts(panl)["KLK6",])

p <- make_escheR(panl,y_reverse=FALSE)
p <- p |> add_fill("log counts",size=0.2,point_size = 0.225)

pdf("manuscript_plots/Fig1/1L-X99_1225B_KLK6.pdf",height=2,width=2)
p+
    scale_fill_continuous("log\ncounts",low= "#FFFFFF",high="#000000")+
    geom_path(data=dompoly,aes(x=xpol,y=ypol,group=dompol,col=dompol),linewidth = 0.25)+
    scale_color_manual(values=pal,na.value = NA)+
    ggtitle(paste0("*KLK6*",mscriptids[sample_id=="X99_1225B",manuscript_id]))+
    guides(color="none")+
    # labs(col="Xenium\nDomain")+
    theme(plot.title.position = "plot",plot.title = element_markdown(size=9,hjust=0.5,margin=margin(0.12,0,-0.12,0,"in")),legend.text = element_text(size=7),legend.key.size = ggplot2::unit(0.125,"in"),legend.title=element_text(size=7.5,hjust=0.5),plot.margin = margin(-0.11,-0.4,-0.12,-0.4,unit = "in"),legend.margin=margin(0,-0.2,0,0,"in"),legend.background = element_blank(),legend.box = element_blank(),legend.key = element_blank())
dev.off()

rm(panl,p)
```