{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Applying `MCube` to the rotated dataset" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "set.seed(20240709)\n", "\n", "library(MCube)\n", "library(ggplot2)\n", "\n", "max_cores <- 36" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "DATA_PATH <- \"/import/home/share/zw/pql/data/mouse_brain\"\n", "RESULT_PATH <- \"/import/home/share/zw/pql/results/mouse_brain\"\n", "\n", "if (!dir.exists(file.path(RESULT_PATH, \"visium_1\"))) {\n", " dir.create(file.path(RESULT_PATH, \"visium_1\"), recursive = TRUE)\n", "}\n", "if (!dir.exists(file.path(RESULT_PATH, \"visium_2\"))) {\n", " dir.create(file.path(RESULT_PATH, \"visium_2\"), recursive = TRUE)\n", "}\n", "if (!dir.exists(file.path(RESULT_PATH, \"ST_3D\"))) {\n", " dir.create(file.path(RESULT_PATH, \"ST_3D\"), recursive = TRUE)\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visium Slice 1" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# reference_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_1\", \"reference_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))\n", "# dim(reference_ST)\n", "# num_celltypes <- nrow(reference_ST)\n", "\n", "# proportions_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_1\", paste0(\"prop_slice\", 0, \".csv\")),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))\n", "# dim(proportions_ST)\n", "\n", "# counts <- as.data.frame(readr::read_csv(\n", "# file.path(DATA_PATH, \"visium_1\", \"counts.csv\")\n", "# ))\n", "# rownames(counts) <- counts[, 1]\n", "# counts[, 1] <- NULL\n", "# counts <- data.matrix(counts)\n", "# dim(counts)\n", "\n", "# coordinates <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_1\", \"3D_coordinates.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, c(\"x\", \"y\")]\n", "# dim(coordinates)\n", "\n", "# spot_effects_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_1\", \"spot_effects_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, 1]\n", "\n", "# library_sizes_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_1\", \"library_sizes_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, 1]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# # Real case\n", "\n", "# mcube_object_1 <- createMCube(\n", "# counts = counts, coordinates = coordinates, \n", "# proportions = proportions_ST, library_sizes = library_sizes_ST, \n", "# covariates = NULL,\n", "# reference = reference_ST,\n", "# spot_effects = spot_effects_ST, platform_effects = NULL,\n", "# project = \"mouse_brain_visium_1\"\n", "# )\n", "# mcube_object_1 <- mcubeFitNull(\n", "# mcube_object_1,\n", "# num_workers = 70, num_threads = 1\n", "# )\n", "# mcube_object_1 <- mcubeTest(\n", "# mcube_object_1,\n", "# num_workers = 70, num_threads = 1, shared_memory = TRUE\n", "# )\n", "# saveRDS(\n", "# mcube_object_1,\n", "# file = file.path(\n", "# RESULT_PATH, \"visium_1\",\n", "# paste0(\"mcube\", \".rds\")\n", "# )\n", "# )\n", "mcube_object_1 <- readRDS(\n", " file = file.path(\n", " RESULT_PATH, \"visium_1\",\n", " paste0(\"mcube\", \".rds\")\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "mcubeKernel: length_scale is set as 0.101059327818091 for the Gaussian kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.142919472004653 for the Gaussian kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.101059327818091 for the Gaussian_transformed kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.142919472004653 for the Gaussian_transformed kernel.\n", "\n", "Number of physical cores: 72.\n", "\n", "Number of workers: 70.\n", "\n", "Number of thread(s) on BLAS per worker: 1.\n", "\n" ] } ], "source": [ "# Rotate 45 degree\n", "\n", "degrees <- 45\n", "coordinates_rotation <- mcube_object_1@coordinates %*% rotation_matrix_2d(degrees)\n", "rownames(coordinates_rotation) <- rownames(mcube_object_1@coordinates)\n", "mcube_object_1_rotation <- mcube_object_1\n", "mcube_object_1_rotation@coordinates <- coordinates_rotation\n", "\n", "mcube_object_1_rotation <- mcubeTest(\n", " mcube_object_1_rotation,\n", " num_workers = 70, num_threads = 1, shared_memory = TRUE\n", ")\n", "pvalues_1_rotation <- mcube_object_1_rotation@pvalues\n", "rm(mcube_object_1_rotation)\n", "saveRDS(\n", " pvalues_1_rotation,\n", " file = file.path(\n", " RESULT_PATH, \"visium_1\",\n", " paste0(\"pvalues_rotation\", \".rds\")\n", " )\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visium Slice 2" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# reference_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_2\", \"reference_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))\n", "# dim(reference_ST)\n", "# num_celltypes <- nrow(reference_ST)\n", "\n", "# proportions_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_2\", paste0(\"prop_slice\", 0, \".csv\")),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))\n", "# dim(proportions_ST)\n", "\n", "# counts <- as.data.frame(readr::read_csv(\n", "# file.path(DATA_PATH, \"visium_2\", \"counts.csv\")\n", "# ))\n", "# rownames(counts) <- counts[, 1]\n", "# counts[, 1] <- NULL\n", "# counts <- data.matrix(counts)\n", "# dim(counts)\n", "\n", "# coordinates <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_2\", \"3D_coordinates.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, c(\"x\", \"y\")]\n", "# dim(coordinates)\n", "\n", "# spot_effects_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_2\", \"spot_effects_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, 1]\n", "\n", "# library_sizes_ST <- data.matrix(read.csv(\n", "# file.path(DATA_PATH, \"visium_2\", \"library_sizes_ST.csv\"),\n", "# header = TRUE, row.names = 1, check.names = FALSE\n", "# ))[, 1]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# Real case\n", "\n", "# mcube_object_2 <- createMCube(\n", "# counts = counts, coordinates = coordinates, \n", "# proportions = proportions_ST, library_sizes = library_sizes_ST, \n", "# covariates = NULL,\n", "# reference = reference_ST,\n", "# spot_effects = spot_effects_ST, platform_effects = NULL,\n", "# project = \"mouse_brain_visium_2\"\n", "# )\n", "# mcube_object_2 <- mcubeFitNull(\n", "# mcube_object_2,\n", "# num_workers = 70, num_threads = 1\n", "# )\n", "# mcube_object_2 <- mcubeTest(\n", "# mcube_object_2,\n", "# num_workers = 70, num_threads = 1, shared_memory = TRUE\n", "# )\n", "# saveRDS(\n", "# mcube_object_2,\n", "# file = file.path(\n", "# RESULT_PATH, \"visium_2\",\n", "# paste0(\"mcube\", \".rds\")\n", "# )\n", "# )\n", "mcube_object_2 <- readRDS(\n", " file = file.path(\n", " RESULT_PATH, \"visium_2\",\n", " paste0(\"mcube\", \".rds\")\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "mcubeKernel: length_scale is set as 0.104577495815249 for the Gaussian kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.147894912900941 for the Gaussian kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.104577495815249 for the Gaussian_transformed kernel.\n", "\n", "mcubeKernel: length_scale is set as 0.147894912900941 for the Gaussian_transformed kernel.\n", "\n", "Number of physical cores: 72.\n", "\n", "Number of workers: 70.\n", "\n", "Number of thread(s) on BLAS per worker: 1.\n", "\n" ] } ], "source": [ "# Rotate 45 degree\n", "\n", "degrees <- 45\n", "coordinates_rotation <- mcube_object_2@coordinates %*% rotation_matrix_2d(degrees)\n", "rownames(coordinates_rotation) <- rownames(mcube_object_2@coordinates)\n", "mcube_object_2_rotation <- mcube_object_2\n", "mcube_object_2_rotation@coordinates <- coordinates_rotation\n", "\n", "mcube_object_2_rotation <- mcubeTest(\n", " mcube_object_2_rotation,\n", " num_workers = 70, num_threads = 1, shared_memory = TRUE\n", ")\n", "pvalues_2_rotation <- mcube_object_2_rotation@pvalues\n", "rm(mcube_object_2_rotation)\n", "saveRDS(\n", " pvalues_2_rotation,\n", " file = file.path(\n", " RESULT_PATH, \"visium_2\",\n", " paste0(\"pvalues_rotation\", \".rds\")\n", " )\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing the results before and after rotation" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# pvalues_1_rotation <- readRDS(\n", "# file = file.path(\n", "# RESULT_PATH, \"visium_1\",\n", "# paste0(\"pvalues_rotation\", \".rds\")\n", "# )\n", "# )\n", "# pvalues_2_rotation <- readRDS(\n", "# file = file.path(\n", "# RESULT_PATH, \"visium_2\",\n", "# paste0(\"pvalues_rotation\", \".rds\")\n", "# )\n", "# )" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "vscode": { "languageId": "r" } }, "outputs": [ { "data": { "text/html": [ "
| celltype | gene | pvalue | pvalue_rotation | slice | |
|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <chr> | |
| 1 | Ext_Hpc_CA1 | 1110051M20Rik | 2.469750e-01 | 3.437970e-01 | Slice 1 |
| 2 | Ext_Hpc_CA1 | 2010300C02Rik | 7.794415e-07 | 2.162494e-08 | Slice 1 |
| 3 | Ext_Hpc_CA1 | 2900026A02Rik | 5.424490e-01 | 5.456476e-01 | Slice 1 |
| 4 | Ext_Hpc_CA1 | 4930402H24Rik | 1.067728e-05 | 8.406705e-06 | Slice 1 |
| 5 | Ext_Hpc_CA1 | 4932438A13Rik | 3.835030e-01 | 5.071124e-01 | Slice 1 |
| 6 | Ext_Hpc_CA1 | Aak1 | 5.320344e-02 | 6.501347e-02 | Slice 1 |