## ----setup, include=FALSE, echo = FALSE--------------------------------------- knitr::opts_chunk$set(echo = TRUE) library(spammR) ##these are all on cran library(tidyverse) library(tidygraph) library(readxl) library(cowplot) library(ggraph) library(BiocFileCache) library(leapR) ##this is on bioconductor ## ----create cache------------------------------------------------------------- path <- tempfile() bfc <- BiocFileCache(path, ask = FALSE) ## ----get proteomics data, warning=FALSE, message=FALSE------------------------ #brain data finame <- paste0('https://zenodo.org/records/13345212/files/DESI+spatial', '%20proteomics%20rat%20brain%20proteomics%20', 'result.xlsx?download=1') pdat <- bfcadd(bfc, "proteomics", fpath=finame) #read in data, convert to matrix dat <- readxl::read_xlsx(pdat)#'dat.xlsx') feature_data <- select(dat, c(Protein,`Protein ID`,`Entry Name`,Gene)) |> dplyr::distinct() |> tibble::column_to_rownames('Protein') dat <- select(dat, c(Protein, starts_with('ROI'))) |> tibble::column_to_rownames('Protein') ## ----load image, warning=FALSE, message=FALSE--------------------------------- ##read in image img <- system.file("extdata","brain_img_0.png",package = 'spammR') ## ----plot image, eval=FALSE--------------------------------------------------- # # cowplot::ggdraw() + cowplot::draw_image(system.file("extdata", # "brain_img_0.png", # package = "spammR")) # # #ggsave('brain_image.pdf') ## ----format coordinates, warning=FALSE, message=FALSE, eval = FALSE----------- # #remove this once we can include it in spammR # # library(geojsonR) # jsondat <- FROM_GeoJson(system.file('extdata','brain_roi.geojson', # package = 'spammR')) # # ##herr er hry yhr # coords <- do.call(rbind, lapply(jsondat$features,function(x){ # roi <- x$properties$name # xv <- x$geometry$coordinates[,1] # yv <- x$geometry$coordinates[,2] # if (is.null(roi)) # roi = "" # return(list(ID = roi, x_pixels = min(xv), y_pixels = min(yv), # spot_height = max(yv) - min(yv), # spot_width = max(xv) - min(xv))) # })) |> # as.data.frame() |> # subset(ID != "") |> # tidyr::separate(ID,into = c('ROI','Replicate'), # sep = '_', # remove = FALSE) |> # tibble::remove_rownames() |> # tibble::column_to_rownames('ID') # # ##now for each ROI we want an x, y, cell height and cell_width # #y-coordinates are from top, so need to udpate # coords$y_pixels = 263 - unlist(coords$y_pixels) - unlist(coords$spot_height) # coords$x_pixels = unlist(coords$x_pixels) # coords$spot_height = unlist(coords$spot_height) # coords$spot_width = unlist(coords$spot_width) # write.csv(as.data.frame(coords),file='../inst/extdata/bcoords.csv') # ## ----load coords-------------------------------------------------------------- coords <- read.csv(system.file('extdata','bcoords.csv', package = 'spammR')) rownames(coords) <- paste(coords$ROI, coords$Replicate,sep = '_') ## ----create proteomics spe, warning=FALSE, message=FALSE---------------------- #create an SFE spe <- spammR::convert_to_spe(dat = dat, feature_meta = feature_data, sample_meta = coords, spatial_coords_colnames = c('x_pixels','y_pixels'), assay_name = 'proteomics', sample_id = 'rat_brain', image_files = img, image_id = 'rat_brain') # myelin spammR::spatial_heatmap(spe, feature = 'P02688', feature_type = 'Protein ID', sample_id = 'rat_brain', image_id = 'rat_brain') ## ----metaspace import and plotting, warning=FALSE, message=FALSE, eval = FALSE---- # ##get lipid data from metaspace # # mspe <- spammR::retrieve_metaspace_data("2024-02-15_20h37m13s", # fdr = 0.2, # assay_name = 'lipids', # sample_id = 'rat_brain', # rotate = TRUE, # drop_zeroes = TRUE) ## ----udpate ms object, error = FALSE------------------------------------------ mspe.url <- 'https://api.figshare.com/v2/file/download/60993349' mcache <- bfcadd(bfc, "mspe", fpath = mspe.url) mspe <- readRDS(mcache) ## ----add image---------------------------------------------------------------- ## add in image to check mspe <- SpatialExperiment::addImg(mspe, img, scaleFactor = NA_real_, sample_id = 'rat_brain', image_id = 'rat_brain') meta <- rowData(mspe) meta$Name <- vapply(meta$moleculeNames, function(x) { x[1]}, character(1)) rowData(mspe) <- meta spatial_heatmap(mspe, feature = 'C41H78NO7P-H-', assay = 'lipids', plot_log = TRUE, metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain')#plot with an ion #ggsave('lipidExpr.pdf') ## ----reduce, echo=FALSE, warning=FALSE, message=FALSE------------------------- reduced <- spat_reduce(spe_target = spe, spe_origin = mspe, origin_assay = 'lipids') spatial_heatmap(reduced, feature = 'C41H78NO7P-H-', assay = 'lipids', plot_log = TRUE, metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain')#plot with an ion #ggsave('lipidRedExpr.pdf') ## ----gfap cors, warning = FALSE, message = FALSE------------------------------ #let's look at the correlation between the protein and the lipids mbp_cor <- spatial_network(reduced, assay_names = c("proteomics","lipids"), query_feature = 'P47819', #protein target_features = rowData(altExp(reduced))[,'Name'],#lipids feature_names = c('Protein ID','Name')) mbp_cor |> activate(edges) |> arrange(desc(abs(corval))) mbp_cor |> ggraph(layout = 'stress') + geom_edge_link(aes(colour = corval)) + geom_node_point(aes(color = class)) + geom_node_label(aes(label = name, color = class)) + scale_color_manual(values=c(lipids = 'darkorange',proteomics = 'slateblue3')) #ggsave('gfap_cor_graph.pdf') nodes <- mbp_cor |> activate(nodes) |> as_tibble() nodes[c(4,2),] ## ----lipid heatmap------------------------------------------------------------ spatial_heatmap(reduced, feature = 'P47819', assay_name = 'proteomics', feature_type = 'Protein ID', metric_display = 'Protein expression', plot_title = 'Gfap expression', sample_id = 'rat_brain', image_id = 'rat_brain') #ggsave('gfap_brain.pdf') ## ----lipid cor---------------------------------------------------------------- ## anti-correlated a1 <- spatial_heatmap(reduced, feature = '22:6(4Z,7Z,10Z,13Z,16Z,19Z)', assay_name = 'lipids', feature_type = 'Name', metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain') a2 <- spatial_heatmap(mspe, feature = '22:6(4Z,7Z,10Z,13Z,16Z,19Z)', assay_name = 'lipids', feature_type = 'Name', metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain') cowplot::plot_grid(a1, a2, ncol=2) #ggsave('antiCorrelatedLipids.pdf', width = 10) #correlated c1 <- spatial_heatmap(reduced, feature = 'PI(18:1(11Z)/20:3(8Z,11Z,14Z))', assay_name = 'lipids', feature_type = 'Name', metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain') c2 <- spatial_heatmap(mspe, feature = 'PI(18:1(11Z)/20:3(8Z,11Z,14Z))', assay_name = 'lipids', feature_type = 'Name', metric_display = 'Lipid expression', sample_id = 'rat_brain', image_id = 'rat_brain') cowplot::plot_grid(c1, c2, ncol = 2) #ggsave('correlatedLipids.pdf', width = 10) ## ----get protein cors--------------------------------------------------------- lvs <- sort(apply(assay(mspe,'lipids'),1,var,na.rm = TRUE),decreasing = TRUE)[1:10] lvn <- rowData(altExp(reduced))[names(lvs),'Name'] lvn[1:5] ## ----plot lipid of interest--------------------------------------------------- #number 1 seems to have an interesting pattern spatial_heatmap(mspe,feature = lvn[1], assay_name = 'lipids', feature_type = 'Name', image_id = 'rat_brain', sample_id = 'rat_brain', metric_display = 'Lipid expression', plot_log = TRUE) spatial_heatmap(reduced,feature = lvn[1], assay_name = 'lipids', feature_type = 'Name', image_id = 'rat_brain', sample_id = 'rat_brain', metric_display = 'Lipid expression', plot_log = TRUE) ## ----correlation network, warning=FALSE--------------------------------------- #let's look at the correlation between the lipid and proteins fl <- c(lvn[1], rowData(reduced)[,'Gene']) cg <- spatial_network(reduced, assay_names = c("proteomics","lipids"), query_feature = lvn[1], #lipid target_features = rowData(reduced)[,'Gene'], feature_names = c('Gene','Name')) ##we can then use trick above to get the neighborhood of our lipid gt <- cg %>% activate(nodes) |> as_tibble() neigh_weights <- cg |> activate(edges) |> #now sure how to or from is assigned, might want #to check both filter(to == which(gt$name == lvn[1])) #then we can get the edge weights and nodes edges <- neigh_weights |> as_tibble() nodes <- neigh_weights |> activate(nodes) |> as_tibble() ##create a table with the protein correlationm #sep2026: proteins used to be in $to column now in $from prot_cor <- data.frame(nodes[edges$from,],lipid_cor = edges$corval) |> dplyr::rename(Gene = 'name') |> subset(Gene %in% rowData(reduced)[,'Gene']) neigh_weights |> activate(edges) |> filter(abs(corval) > 0.75) |> activate(nodes) |> mutate(degree = centrality_degree(mode = 'all')) |> filter(degree > 0) |> ggraph(layout = 'stress') + geom_edge_link(aes(colour = corval)) + geom_node_point(aes(color = class)) + geom_node_label(aes(label = name, color = class)) ## ----leapR, warning=FALSE----------------------------------------------------- ##add back to data set prot_cor <- prot_cor |> group_by(Gene) |> summarize(lipid_cor = mean(lipid_cor)) rowData(reduced) <- rowData(reduced) |> as.data.frame() |> left_join(as.data.frame(prot_cor)) |> mutate(upperGene = toupper(Gene)) data("krbpaths") enrich <- leapR::leapR(geneset = krbpaths, enrichment_method = 'enrichment_in_order', eset = reduced, primary_column = 'lipid_cor', method = 'ks', id_column = 'upperGene', minsize = 5) paths <- subset(enrich, BH_pvalue < 0.3) #subset(enrich,ingroup_mean > 0.25) |> # subset(pvalue <0.05) DT::datatable(paths,options = list(scrollX = TRUE)) ## ----check protein correlation------------------------------------------------ prots <- unlist(stringr::str_split(paths[1,'ingroupnames'], pattern = ', ')) pv <- assay(subset(reduced,upperGene %in% prots),'proteomics') lv <- assay(subset(altExp(reduced),Name == lvn[1]),'lipids') ##proteins seem anti-correlated mean(cor(t(pv),t(lv), method = 'spearman')) ## ----plot proteins------------------------------------------------------------ p2 <- spatial_heatmap(reduced, assay_name = 'proteomics', feature = prots, feature_type = 'upperGene', sample_id = 'rat_brain', image_id = 'rat_brain', plot_log = TRUE, title_size = 10, plot_title = rownames(paths)[1], metric_display = 'Protein expression') p2 ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()