Vascular Slideseq

Introduction

This R Markdown is meant to facilitate the reproducibility of the findings in https://docs.google.com/document/d/166X4o_6HegeS0uUHRa4ahKedkaVEMeiRagCm7T9nOVU/edit?usp=sharing. Here we specifically focus on the vascular and Slide-seq data, as detailed in Figure 4.

First, we load relevant packages.

library(ggplot2)
library(tidyr)
library(dplyr)
library(rgl)
library(plotrix)
'%notin%' = Negate('%in%')

Examine Vascular Density in 3-Dimensions

The vascular samples are named BL6J-1, BL6J-2, and BL6J-3. For these plots, we read in each sample of vascular density, filter out low density voxels for visual clarity, and assign a color to each voxels.

Sample BL6J-1

voxels_Intensity_1 = readRDS("Data/Vascular/DataFiles/100um_vesSAP_OurMeasure_BL6J-1.rds")
#We filter out voxels with a vascular density less than or equal to 0.4 for visual transparency
voxels_Intensity_1 = filter(voxels_Intensity_1, voxels_Intensity_1$meanIntensity > 0.4)
colfunc <- colorRampPalette(c("yellow", "red"))
rainbow = colfunc(10)
voxels_Intensity_1$Bin = cut(voxels_Intensity_1$meanIntensity, breaks = 10)
names(rainbow) = levels(voxels_Intensity_1$Bin)
voxels_Intensity_1$cols = rainbow[match(voxels_Intensity_1$Bin, names(rainbow))]

plot3d(voxels_Intensity_1$Voxel_X, voxels_Intensity_1$Voxel_Y, voxels_Intensity_1$Voxel_Z, col = voxels_Intensity_1$cols, add = TRUE)

Sample BL6J-2

voxels_Intensity_2 = readRDS("Data/Vascular/DataFiles/100um_vesSAP_OurMeasure_BL6J-2.rds")
#We filter out voxels with a vascular density less than or equal to 0.4 for visual transparency
voxels_Intensity_2 = filter(voxels_Intensity_2, voxels_Intensity_2$meanIntensity > 0.4)
colfunc <- colorRampPalette(c("yellow", "red"))
rainbow = colfunc(10)
voxels_Intensity_2$Bin = cut(voxels_Intensity_2$meanIntensity, breaks = 10)
names(rainbow) = levels(voxels_Intensity_2$Bin)
voxels_Intensity_2$cols = rainbow[match(voxels_Intensity_2$Bin, names(rainbow))]

close3d()
plot3d(voxels_Intensity_2$Voxel_X, voxels_Intensity_2$Voxel_Y, voxels_Intensity_2$Voxel_Z, col = voxels_Intensity_2$cols, add = TRUE)

Sample BL6J-3

voxels_Intensity_3 = readRDS("Data/Vascular/DataFiles/100um_vesSAP_OurMeasure_BL6J-3.rds")
#We filter out voxels with a vascular density less than or equal to 0.4 for visual transparency
voxels_Intensity_3 = filter(voxels_Intensity_3, voxels_Intensity_3$meanIntensity > 0.4)
colfunc <- colorRampPalette(c("yellow", "red"))
rainbow = colfunc(10)
voxels_Intensity_3$Bin = cut(voxels_Intensity_3$meanIntensity, breaks = 10)
names(rainbow) = levels(voxels_Intensity_3$Bin)
voxels_Intensity_3$cols = rainbow[match(voxels_Intensity_3$Bin, names(rainbow))]

close3d()
plot3d(voxels_Intensity_3$Voxel_X, voxels_Intensity_3$Voxel_Y, voxels_Intensity_3$Voxel_Z, col = voxels_Intensity_3$cols, add = TRUE)

Calculating a Vascular Score for each Molecular Cell Type

First, we have to calculate a vascular score for each molecular cell type across the three vascular samples (BL6J-1, BL6J-2, BL6J-3).

rctd_reduced = readRDS("Data/Vascular/SupportFiles/100um_Celltypes_in_each_Voxel_and_totalsums.rds")
#Filter for only voxels that contain atleast one Slide-seq bead with a confident assignment
rctd_reduced = rctd_reduced[rctd_reduced[,"total_sums"] != 0,] 
vasdata = list.files("Data/Vascular/DataFiles/")
for(x in vasdata){
  #Read in each vascular sample
  voxelsIntensity = readRDS(paste0("Data/Vascular/DataFiles/", x))
  #Use the voxel coordinates to make a data frame of the cell types present in each voxel, as well as the vascular intensity within that voxel
  rctdVascular = merge(voxelsIntensity, rctd_reduced)
  running_list = c()
  cell_types = subset(colnames(rctdVascular), colnames(rctdVascular) %notin% c("Voxel_X", "Voxel_Y","Voxel_Z","Voxel_ID", "meanIntensity","total_sums"))
  for(typer in cell_types){
    ratio = rctdVascular[,typer]/rctdVascular[, "total_sums"]
    intensity = rctdVascular$meanIntensity
    products = ratio*intensity
    vascular_score = sum(products)/ sum(rctdVascular[,typer] != 0) 
    names(vascular_score) = typer
    running_list = c(running_list, vascular_score)
  }
  x = sub(".rds", "", x)
  saveRDS(running_list, paste0("Data/Vascular/Results/VascularScore_", x, "_100um_Slideseq.RDS"))
}

Manhattan Plots

To create the manhattan plots, we first need to get a list of all cell types that occupy more than 15 voxels

  cellCounts = readRDS("Data/Vascular/SupportFiles/NumberVoxels_SlideseqCellTypes_subsetforvascularspace.rds")
  #Filter for cell types where the number of number of voxels confidentally assigned this cell types is greater than 15
  cellCounts = filter(cellCounts, cellCounts$NumVoxels > 15)

  sample = readRDS("Data/Vascular/Results/VascularScore_100um_vesSAP_OurMeasure_BL6J-1_100um_Slideseq.RDS")
  #Reformat
  sample = data.frame(sample)
  colnames(sample) = "VascularScore"
  sample$CellTypes = rownames(sample)
  sample = left_join(cellCounts, sample)
  sample = filter(sample, !is.na(sample$VascularScore))

 manhplot <- ggplot(sample, aes(x = CellTypes, y = VascularScore, color = Broad)) +
    geom_point(alpha = 0.75, size = 0.8) + labs(x = NULL,  y = "Vascular Scores", title = "BL6J-1 vesSAP") + theme_test() + scale_x_discrete( expand = c(-1.03, 1.03) ) + geom_text(aes(label=ifelse(VascularScore >0.12 | VascularScore < 0.01,as.character(CellTypes),'')),hjust=0,vjust=0, show_guide = F) + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) 
    print(manhplot)

sample = readRDS("/nfs/turbo/umms-welchjd/BRAIN_initiative/Final_integration_workflow/Markdowns/Data/Vascular/Results/VascularScore_100um_vesSAP_OurMeasure_BL6J-2_100um_Slideseq.RDS")
  #Reformat
  sample = data.frame(sample)
  colnames(sample) = "VascularScore"
  sample$CellTypes = rownames(sample)
  sample = left_join(cellCounts, sample)
  sample = filter(sample, !is.na(sample$VascularScore))

 manhplot <- ggplot(sample, aes(x = CellTypes, y = VascularScore, color = Broad)) +
    geom_point(alpha = 0.75, size = 0.8) + labs(x = NULL,  y = "Vascular Scores", title = "BL6J-2 vesSAP") + theme_test() + scale_x_discrete( expand = c(-1.03, 1.03) ) + geom_text(aes(label=ifelse(VascularScore >0.12 | VascularScore < 0.01,as.character(CellTypes),'')),hjust=0,vjust=0, show_guide = F) + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) 
    print(manhplot)

 sample = readRDS("/nfs/turbo/umms-welchjd/BRAIN_initiative/Final_integration_workflow/Markdowns/Data/Vascular/Results/VascularScore_100um_vesSAP_OurMeasure_BL6J-3_100um_Slideseq.RDS")
  #Reformat
  sample = data.frame(sample)
  colnames(sample) = "VascularScore"
  sample$CellTypes = rownames(sample)
  sample = left_join(cellCounts, sample)
  sample = filter(sample, !is.na(sample$VascularScore))

 manhplot <- ggplot(sample, aes(x = CellTypes, y = VascularScore, color = Broad)) +
    geom_point(alpha = 0.75, size = 0.8) + labs(x = NULL,  y = "Vascular Scores", title = "BL6J-3 vesSAP") + theme_test() + scale_x_discrete( expand = c(-1.03, 1.03) ) + geom_text(aes(label=ifelse(VascularScore >0.12 | VascularScore < 0.01,as.character(CellTypes),'')),hjust=0,vjust=0, show_guide = F) + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) 
    print(manhplot)

Creating Cell Type Specific Plots

We next examine the spatial distribution of two molecular cell types with a high vascular score, and one molecular cell type with a low vascular score. Specifically, for each 100 um voxel, if our cell type of interest is present, we color that voxel blue.

#Read in the voxelized Slide-seq data
rctd_reduced = readRDS("Data/Vascular/SupportFiles/100um_Celltypes_in_each_Voxel_and_totalsums.rds")

We first look at the molecular cell type Ex_Sall3_Cd44_2

celltype = "Ex_Sall3_Cd44_2"
cd44 = rctd_reduced[,c(celltype, "Voxel_X", "Voxel_Y", "Voxel_Z")]
cd44 = filter(cd44, cd44[,1] > 0)

Cell Type: Ex_Sall3_Cd44_2 Sample: BL6J-1

close3d()
plot3d(cd44$Voxel_X, cd44$Voxel_Y, cd44$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_1$Voxel_X, voxels_Intensity_1$Voxel_Y, voxels_Intensity_1$Voxel_Z, col = voxels_Intensity_1$cols, add = TRUE)

Cell Type: Ex_Sall3_Cd44_2 Sample: BL6J-2

close3d()
plot3d(cd44$Voxel_X, cd44$Voxel_Y, cd44$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_2$Voxel_X, voxels_Intensity_2$Voxel_Y, voxels_Intensity_2$Voxel_Z, col = voxels_Intensity_2$cols, add = TRUE)

Cell Type: Ex_Sall3_Cd44_2 Sample: BL6J-3

close3d()
plot3d(cd44$Voxel_X, cd44$Voxel_Y, cd44$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_3$Voxel_X, voxels_Intensity_3$Voxel_Y, voxels_Intensity_3$Voxel_Z, col = voxels_Intensity_3$cols, add = TRUE)

We next look at the molecular cell type Ex_Fezf2_Il22

celltype = "Ex_Fezf2_Il22"
fezf2 = rctd_reduced[,c(celltype, "Voxel_X", "Voxel_Y", "Voxel_Z")]
fezf2 = filter(fezf2, fezf2[,1] > 0)

Cell Type: Ex_Fezf2_Il22 Sample: BL6J-1

close3d()
plot3d(fezf2$Voxel_X, fezf2$Voxel_Y, fezf2$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_1$Voxel_X, voxels_Intensity_1$Voxel_Y, voxels_Intensity_1$Voxel_Z, col = voxels_Intensity_1$cols, add = TRUE)

Cell Type: Ex_Fezf2_Il22 Sample: BL6J-2

close3d()
plot3d(fezf2$Voxel_X, fezf2$Voxel_Y, fezf2$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_2$Voxel_X, voxels_Intensity_2$Voxel_Y, voxels_Intensity_2$Voxel_Z, col = voxels_Intensity_2$cols, add = TRUE)

Cell Type: Ex_Fezf2_Il22 Sample: BL6J-3

close3d()
plot3d(fezf2$Voxel_X, fezf2$Voxel_Y, fezf2$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_3$Voxel_X, voxels_Intensity_3$Voxel_Y, voxels_Intensity_3$Voxel_Z, col = voxels_Intensity_3$cols, add = TRUE)

Finally, we look at the molecular cell type Ex_Pou4f1_Htr5b_3

celltype = "Ex_Pou4f1_Htr5b_3"
htr5b = rctd_reduced[,c(celltype, "Voxel_X", "Voxel_Y", "Voxel_Z")]
htr5b = filter(htr5b, htr5b[,1] > 0)

Cell Type: Ex_Pou4f1_Htr5b_3 Sample: BL6J-1

close3d()
plot3d(htr5b$Voxel_X, htr5b$Voxel_Y, htr5b$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_1$Voxel_X, voxels_Intensity_1$Voxel_Y, voxels_Intensity_1$Voxel_Z, col = voxels_Intensity_1$cols, add = TRUE)

Cell Type: Ex_Pou4f1_Htr5b_3 Sample: BL6J-2

close3d()
plot3d(htr5b$Voxel_X, htr5b$Voxel_Y, htr5b$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_2$Voxel_X, voxels_Intensity_2$Voxel_Y, voxels_Intensity_2$Voxel_Z, col = voxels_Intensity_2$cols, add = TRUE)

Cell Type: Ex_Pou4f1_Htr5b_3 Sample: BL6J-3

close3d()
plot3d(htr5b$Voxel_X, htr5b$Voxel_Y, htr5b$Voxel_Z, col = "blue", add = TRUE)
plot3d(voxels_Intensity_3$Voxel_X, voxels_Intensity_3$Voxel_Y, voxels_Intensity_3$Voxel_Z, col = voxels_Intensity_3$cols, add = TRUE)