GridDBを用いたRでの遺伝子発現データセットの取り込みとクエリー

git clone https://github.com/griddbnet/Blogs.git --branch gene_analysis

attachment::att_from_rmd("01.Rmd")

deps <- attachment::att_from_rmd("01.Rmd")
install.packages(deps)

# replace `PKG1/2/3` with the required package names
sudo apt install r-cran-PKG1 r-cran-PKG2 r-cran-PKG3 ...

suppressPackageStartupMessages({
  library(RJDBC)
  library(DBI)
  library(dplyr)
  library(tidyr)
  library(purrr)
  library(qs)
})
connect_to_griddb <- function() {
  drv <- JDBC(driverClass = "com.toshiba.mwcloud.gs.sql.Driver",
              # Point this to your gridstore jar
              classPath = "~/src/jdbc/bin/gridstore-jdbc.jar")
  
  conn <-
    dbConnect(drv,
              "jdbc:gs://127.0.0.1:20001/myCluster/public",
              "admin",
              "admin")
  
  return(conn)
}

TCGA <- qs::qread("data/expressions.qs")
sapply(TCGA, nrow) |> summary()

TCGA[[1]] |> head() |> knitr::kable()

sapply(TCGA, "[[", 1, 5) |> unlist() |> knitr::kable()

griddb <- connect_to_griddb()
create_table <- function(table_name) {
  dbSendUpdate(
    griddb,
    sprintf(
      "CREATE TABLE IF NOT EXISTS %s (Gene STRING, Sample STRING, Expression INTEGER, Code STRING, Name STRING);",
      table_name
    )
  )
}
insert_table <- function(conn, name, df, append = TRUE) {
  for (i in seq_len(nrow(df))) {
    dbWriteTable(conn, name, df[i, ], append = append)
  }
}
# using `invisible` here to hide non-informative output
invisible(lapply(names(TCGA), create_table))
DBI::dbListTables(griddb)
# 6 minutes for 10K rows per table. So 6 min for 100K rows
# 1h for 100K rows per table. So 1h for 1M rows.
invisible(pbapply::pblapply(seq_along(TCGA), function(i)
  insert_table(
    conn = griddb,
    name = names(TCGA)[[i]],
    df = TCGA[[i]][1:100000, ],
    append = TRUE
  ))
)
DBI::dbDisconnect(griddb)

griddb <- connect_to_griddb()
iter <- setNames(names(TCGA), names(TCGA))
map_dfr(iter, function(x) {
  griddb |> tbl(x) |> collect() |> nrow()
}) |> knitr::kable()

high_expression_genes <- purrr::map_dfr(iter, .id = "Cancer", 
  function(x) {
    griddb |>
      tbl(x) |>
      group_by(Gene) |>
      summarise(Mean_expression = mean(Expression, na.rm = TRUE)) |>
      filter(Mean_expression >= 2000) |> 
      collect()
  }
)
dim(high_expression_genes)
head(high_expression_genes) |> knitr::kable()

low_expression_genes <- purrr::map_dfr(iter, .id = "Cancer",
  function(x) {
    griddb |>
      tbl(x) |>
      group_by(Gene) |>
      summarise(Mean_expression = mean(Expression, na.rm = TRUE)) |>
      filter(Mean_expression <= 300) |> 
      collect() 
  }
)
dim(low_expression_genes)
tail(low_expression_genes) |> knitr::kable()

unique_genes_per_dataset <- purrr::map(iter, 
  function(x) {
    griddb |>
      tbl(x) |>
      collect() |>
      pull(Gene) |> 
      unique()
  }
)
shared_genes <- Reduce(intersect, unique_genes_per_dataset)
shared_genes

first10 <- shared_genes[1:10]
iter2 <- expand.grid(iter, first10, stringsAsFactors = FALSE)
first10_expression <-
  purrr::map2_dfr(iter2$Var1, iter2$Var2, function(.x, .y) {
    griddb |>
      tbl(.x) |>
      filter(Gene == .y) |> 
      collect() 
  })
first10_expression |>
  group_by(Code, Name, Gene) |>
  summarise(Mean_expression = mean(Expression, na.rm = TRUE)) |>
  ungroup() |> 
  select(-Name) |>
  pivot_wider(names_from = Code, values_from = Mean_expression) |> 
  knitr::kable()

DBI::dbDisconnect(griddb)

# This has a lot of dependencies and can take a few minutes to install
# BiocManager::install("GSEABenchmarkeR")
library(GSEABenchmarkeR)
tcga <- loadEData("tcga", nr.datasets = 10)
cancer_abbreviations <- names(tcga)
cancer_names <- c("BLCA" = "Bladder Urothelial Carcinoma",
                  "BRCA" = "Breast Invasive Carcinoma",
                  "COAD" = "Colon Adenocarcinoma",
                  "HNSC" = "Head and Neck Squamous Cell Carcinoma",
                  "KICH" = "Kidney Chromophobe",
                  "KIRC" = "Kidney Renal Clear Cell Carcinoma",
                  "KIRP" = "Kidney Renal Papillary Cell Carcinoma",
                  "LIHC" = "Liver Hepatocellular Carcinoma",
                  "LUAD" = "Lung Adenocarcinoma",
                  "LUSC" = "Lung Squamous Cell Carcinoma")
expressions <- lapply(seq_along(tcga), function(i) {
  tcga[[i]]@assays@data@listData |> 
    as.data.frame() |> 
    tibble::rownames_to_column(var = "Gene") |> 
    tidyr::pivot_longer(cols = -Gene, names_to = "Sample", values_to = "Expression") |> 
    mutate(cancer = cancer_abbreviations[[i]], cancer_name = cancer_names[[i]])
})
qs::qsave(expressions, "expressions.qs")

sessionInfo()