In dmcable/RCTD: SpatialeXpressionR: Cell type identification and cell type-specific differential expression in spatial transcriptomics

knitr::opts_chunk$set(echo = T, warning = F, message = F, cache = T, results = 'hide')

CSIDE on the MERFISH Hypothalamus

Load in CSIDE results and calculate significant genes

# Load in spatialRNA data and Reference data
library(spacexr)
library(Matrix)
library(devtools)
library(ggplot2)
library(ggpubr)
library(reshape2)
library(dplyr)
library(ggrepel)
library(fields)
library(stringr)
source('../helper_functions/merge_de_helper.R')
dir.exists('../../../slideseq/Cell Demixing/ContentStructure')
load_all()
pwd = getwd()
datadir <- paste0('../../data/moffitt','/')
resultsdir <- paste0('../../../slideseq/Cell Demixing/ContentStructure/DEGLAM/results/ResultsMerfish','/')
myRCTD = readRDS(paste0(resultsdir,'myRCTDde.rds'))
cell_types_present <- myRCTD@internal_vars_de$cell_types_present
cell_types <- myRCTD@internal_vars_de$cell_types
gene_fits <- myRCTD@de_results$gene_fits

Cell type comparison

## Compare excitatory to inhibitory
cell_type_1 <- 'Excitatory'
cell_type_2 <- 'Inhibitory'
ct_ind_1 <- 2*which(cell_types == cell_type_1)
ct_ind_2 <- 2*which(cell_types == cell_type_2)
same_genes <- intersect(get_gene_list_type_wrapper(myRCTD, cell_type_1, cell_types_present),
                        get_gene_list_type_wrapper(myRCTD, cell_type_2, cell_types_present))
gene_fits$con_mat[same_genes,] ### all converged :)
rm_syt2 <- setdiff(same_genes, 'Syt2')
m1 <- gene_fits$mean_val[same_genes,cell_type_1]
m2 <- gene_fits$mean_val[same_genes,cell_type_2]
plot(m1,m2)
cor(m1, m2)^2
cor(gene_fits$mean_val[rm_syt2,cell_type_1], gene_fits$mean_val[rm_syt2,cell_type_2])^2

CSIDE without segmentation

resultsdir <- '../../data/SpatialRNA/MERFISH_24'
resultsdir_seg <- file.path(resultsdir, 'seg')
myRCTD <- readRDS(file.path(resultsdir,'myRCTDde_seg.rds'))
my_barc <- rownames(myRCTD@results$results_df)[which((myRCTD@results$results_df$first_type == 'monocyte/DC' & myRCTD@results$results_df$spot_class != 'reject') | (myRCTD@results$results_df$second_type == 'monocyte/DC' & myRCTD@results$results_df$spot_class == 'doublet_certain' ))]
p1 <- plot_puck_continuous(myRCTD@spatialRNA, colnames(myRCTD@spatialRNA@counts) , myRCTD@internal_vars_de$X2[,2], 
                     title ='') + geom_point(data = myRCTD@spatialRNA@coords[my_barc,], size = 0.6)+ ggplot2::scale_shape_identity() + ggplot2::theme_classic() + ggplot2::scale_size_identity() + coord_fixed() + theme(legend.position="top")+  geom_segment(aes(x = 1934.6, y = -3990, xend = 2184.6, yend = -3990), color = "black") + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())+ scale_color_gradientn(colors = pals::brewer.blues(20)[2:20],name = "Distance from midline", labels = c(0,1),breaks = c(0,1), limits = c(0,1))
results_df <- myRCTD@results$results_df
puck <- myRCTD@spatialRNA
barcodes = rownames(results_df[results_df$spot_class != "reject" & puck@nUMI >= 100,])
my_table = puck@coords[barcodes,]
my_table$class = results_df[barcodes,]$first_type
n_levels = myRCTD@cell_type_info$info[[3]]
my_pal = pals::kelly(n_levels+1)[2:(n_levels+1)]
names(my_pal) = myRCTD@cell_type_info$info[[2]]
my_pal_curr <- my_pal
cell_types <- myRCTD@internal_vars_de$cell_types
my_pal_curr[cell_types[1]] <- "#CC79A7"
my_pal_curr[cell_types[2]] <- "#E69F00"
my_pal_curr[cell_types[3]] <- "#D55E00"
my_pal_curr[cell_types[5]] <- "#009E73"
my_pal_curr[cell_types[4]] <- "#0072B2"
my_pal_curr[cell_types[6]] <- "#56B4E9"
pres = unique(as.integer(my_table$class))
pres = pres[order(pres)]
p2 <- ggplot2::ggplot(my_table, ggplot2::aes(x=x, y=y)) + ggplot2::geom_point(ggplot2::aes(size = .4, shape=19,color=class)) + ggplot2::scale_color_manual("",values = my_pal_curr[pres], breaks = cell_types, labels = cell_types)+ ggplot2::scale_shape_identity() + ggplot2::theme_classic() + ggplot2::scale_size_identity() + coord_fixed() + theme(legend.position="top")+ guides(colour = guide_legend(override.aes = list(size=2)))+ geom_segment(aes(x = -2184.6, y = -4000, xend = -1934.6, yend = -4000), color = "black") + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
ggarrange(p1,p2,nrow = 2)

dmcable/RCTD documentation built on Feb. 24, 2024, 11:03 p.m.