# Processinglibrary(Seurat)library(data.table)library(spacexr)library(tidyverse)##### The following create a "minimal" Seurat object that can be used for filtering.metadata<-as.data.frame(data.table::fread(paste0(dir.path,"cells.csv.gz")))# Read cell data###### Read the gene expression data and create Seurat objectdata<-Read10X(paste0(dir.path,"cell_feature_matrix/"))srat<-CreateSeuratObject(data$`Gene Expression`, assay ="Xenium", meta.data =metadata)##### Subset for cells with at least 10 detected transcripts and a positive cell areasrat<-subset(srat, cells =WhichCells(srat, expression =transcript_counts>=10&cell_area>1))##### Construct the reference for RCTD annotation. In this specific case we are using the ATRT Multiome data##### multiome<-readRDS("Multiome.rds")# Seurat object with the Multiome datacounts<-GetAssayData(multiome, assay ="RNA", slot ="counts")cluster<-as.factor(multiome$Tumor_Annotation)names(cluster)<-colnames(multiome)nUMI<-multiome$nCount_RNAnames(nUMI)<-colnames(multiome)nUMI<-colSums(counts)reference<-Reference(counts, cluster, nUMI)##### Construct the RCTD query objectcells<-read.csv("cells_stats.csv")# Vector of cell ID's corresponding to the sample of interest. Excludes cells manually assigned as "Necrotic"srat<-subset(srat, cells =cells$Cell.ID)query.counts<-GetAssayData(srat, assay ="Xenium", slot ="counts")coords<-FetchData(srat, c("x_centroid","y_centroid")); colnames(coords)<-c("x","y")query<-SpatialRNA(coords, query.counts, colSums(query.counts))rm(cells, coords, query.counts, srat);gc()RCTD<-create.RCTD(query, reference, max_cores =8, UMI_min =10, UMI_min_sigma =10)##### Run RCTD and add annotations to the Seurat objectRCTD<-run.RCTD(RCTD, doublet_mode ="doublet")annotations.df<-RCTD@results$results_dfannotations<-annotations.df$first_typenames(annotations)<-rownames(annotations.df)annotations<-as.data.frame(annotations)colnames(annotations)<-"group"srat<-AddMetaData(srat, metadata =annotations)##### Normalize counts by cell areasrat<-GetAssayData(srat, assay ="Xenium", layer ="counts")norm<-sweep(raw, 2, srat$cell_area, FUN ='/')srat<-SetAssayData(object =srat, layer ="data", new.data =norm, assay ="Xenium")
# Code is commented because of compilation errors with Quarto.# # Neighborhood analysis# # The following code was used to perform Neighborhood Enrichment analysis and export the results to CSV files for downstream analysis in R# # Run independently for each tumor sample# # import scimap as sm# import anndata as ad# import pandas as pd# import scanpy as sc# import squidpy as sq# # # Read expression matrix# adata = sc.read_10x_h5(# filename = "cell_feature_matrix.h5"# )# # # Read the cell info file# df = pd.read_csv(# "cells.csv.gz"# )# # df.set_index(adata.obs_names, inplace = True)# adata.obs = df.copy()# # adata.obsm["spatial"] = adata.obs[["x_centroid","y_centroid"]].copy().to_numpy()# # # Read RCTD annotations# annotations = pd.read_csv(# "ANNOTATIONS.csv"# )# # # Keep only malignant cells# annotations = annotations[~annotations.group.isin(["Necrotic","Mural","Endothelial","OPC","Microglia_Immune","Astrocytes","Neurons"])]# # adata_sub = adata[annotations["cell_id"]].copy()# del(adata)# # annotations.set_index(adata_sub.obs_names, inplace = True)# adata_sub.obs["type"]= annotations["group"].copy().astype('category')# # adata_sub.obs['imageid'] = 'imageid'# # adata_sub = sm.tl.spatial_interaction (adata_sub, # method='knn', # knn=10, # label='spatial_interaction_knn',# x_coordinate = 'x_centroid', # y_coordinate = 'y_centroid',# phenotype = 'type')# # adata_sub.uns['spatial_interaction_knn'].to_csv("Output.csv", index=True)# # # Export the interaction frequencies# sq.gr.spatial_neighbors(adata_sub, coord_type = "generic", delaunay = False, n_neighs = 10)# # pd.DataFrame(adata_sub.uns["type_interactions"]).to_csv("interactions.csv")