05_cmap_data.Rmd

---
title: "CMap data"
author: "Clemens Hug"
date: "3/19/2019"
output: html_document
---

```{r setup, include=FALSE}
library(tidyverse)
library(cmapR)
#library(webchem)
library(synExtra)
library(here)
library(lspcheminf)
library(batchtools)
library(furrr)

synapser::synLogin()
syn <- synExtra::synDownloader(here("data"))
```


```{r compound_info}
compound_dict <- syn("syn20835543") %>%
  read_rds()

fingerprints <- syn("syn21042105") %>%
  read_rds()

# Using previously done pertubation meta and signature meta
# pm <- syn("syn21547097") %>%
#   read_csv()
# 
# im <- syn("syn21547100") %>%
#   read_csv()
# 
# sm <- syn("syn21547101") %>%
#   read_csv()

```


## CMap data

There are two big CMap datasets, Phase I and II with two different GEO accession
numbers. Loading metadata and Z-score matrix from both.


```{r perturbation_meta}
# Warning OK
pertubation_meta_raw <- list(
  GSE92742 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_pert_info.txt.gz",
  GSE70138 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_pert_info_2017-03-06.txt.gz",
  LINCS_2020 = syn("syn25049483")
) %>%
  map(read_tsv, na = "-666") %>%
  bind_rows(.id = "dataset")

run_canonicalize_job <- function(compound_df, timeout = 5*60) {
  input_tmp <- tempfile()
  output_tmp <- tempfile()
  readr::write_csv(compound_df, input_tmp)
  script_path <- tempfile("shell_script_", fileext = ".sh")
  script <- c(
    "unset PYTHONPATH",
    "unset PYTHONHOME",
    "source ~/miniconda3/etc/profile.d/conda.sh",
    "conda activate lspcheminf_env",
    "which -a python",
    "env",
    paste(
      "tautomers canonicalize",
      "--standardizer molvs",
      "--timeout", timeout,
      "--compound-encoding", "smiles",
      input_tmp,
      output_tmp
    )
  )
  message("SCRIPT:\n", paste(paste0("# ", script), collapse = "\n"), "\n")
  writeLines(
    script,
    con = script_path
  )
  out <- sys::exec_wait("bash", script_path)
  if (out != 0)
    stop("Canonicalization failed with error code ", out)
  readr::read_csv(output_tmp)
}

unique_compounds <- pertubation_meta_raw %>%
  drop_na(canonical_smiles) %>%
  distinct(pert_id, compound = canonical_smiles)

unique_compounds_chunked <- unique_compounds %>%
  split(rep_len(1:50, nrow(.)))

reg <- makeRegistry(
  file.dir = here(paste0("registry_", gsub(" ", "_", Sys.time()))),
  seed = 1
)

batchMap(
  fun = run_canonicalize_job,
  compound_df = unique_compounds_chunked,
  # Limiting canonicalization to 5 min per compound
  more.args = list(timeout = 5*60)
)

job_table <- findJobs() %>%
  # Chunk jobs into a single array job
  mutate(chunk = 1)

submitJobs(
  job_table,
  resources = list(
    memory = "6gb",
    ncpus = 1L,
    partition = "short",
    walltime = 15*60,
    chunks.as.arrayjobs = TRUE,
    # For some reason these nodes fail to execute R because of an "illegal instruction"
    exclude = "compute-f-17-[09-25]"
  )
)

unique_compounds_canonical <- unique_compounds_chunked %>%
  enframe("chunk", "raw") %>%
  mutate(
    canonical = map2(
      chunk, raw,
      ~.y %>%
        mutate(row = 1:n() - 1) %>%
        left_join(
          loadResult(as.integer(.x)),
          by = "row"
        )
    )
  ) %>%
  pull(canonical) %>%
  bind_rows() %>%
  select(
    pert_id, inchi
  )

write_csv(
  unique_compounds_canonical,
  here("wrangled", "cmap_compounds_canonical.csv.gz")
)

# unique_compounds_canonical <- read_rds(
#   here("wrangled", "cmap_compounds_canonical.rds")
# )

unique_compounds_fingerprints <- unique_compounds_canonical %>%
  drop_na(inchi) %>%
  crossing(
    fingerprint_type = c("morgan", "topological")
  ) %>%
  group_nest(fingerprint_type) %>%
  mutate(
    data = map2(
      data, fingerprint_type,
      function(df, fpt) {
        q <- set_names(df[["inchi"]], df[["pert_id"]])
        lspcheminf::calculate_fingerprints(q, fpt)
      }
    )
  ) %>%
  unnest(data)

write_csv(
  unique_compounds_fingerprints,
  here("wrangled", "cmap_compound_fingerprints.csv.gz")
)

# unique_compounds_fingerprints <- read_rds(
#   here("wrangled", "cmap_compound_fingerprints.rds")
# )

lsp_fp <- fingerprints %>%
  filter(fp_name == "morgan_normal") %>%
  chuck("data", 1) %>%
  filter(fp_name != "morgan_chiral") %>%
  inner_join(
    compound_dict %>%
      filter(fp_name == "morgan_normal") %>%
      chuck("data", 1) %>%
      select(lspci_id, inchi),
    by = "lspci_id"
  ) %>%
  drop_na() %>%
  group_nest(fingerprint_type = fp_type, .key = "reference")

compound_lspci_id_mapping <- unique_compounds_fingerprints %>%
  inner_join(
    unique_compounds_canonical,
    by = c("names" = "pert_id")
  ) %>%
  drop_na() %>%
  group_nest(fingerprint_type) %>%
  left_join(
    lsp_fp,
    by = "fingerprint_type"
  ) %>%
  mutate(
    data = map2(
      data, reference,
      function(d, r) {
        dq <- set_names(d[["inchi"]], d[["names"]]) %>%
          compounds(fingerprints = d[["fingerprints"]])
        rq <- set_names(r[["inchi"]], as.character(r[["lspci_id"]])) %>%
          compounds(fingerprints = r[["fingerprint"]])
        # browser()
        lspcheminf::chemical_similarity_threshold(
          dq, rq, threshold = 1,
          n_threads = 1
        )
      }
    )
  )

compound_dict_chunks <- compound_dict %>%
  filter(fp_name == "morgan_normal") %>%
  chuck("data", 1) %>%
  drop_na(inchi) %>%
  select(lspci_id, inchi) %>%
  chunk_df(15, seed = 1)

plan(sequential)
compound_dict_mw <- compound_dict_chunks %>%
  future_map(
    ~lspcheminf::molecular_mass(
      set_names(.x$inchi, .x$lspci_id)
    ),
    .progress = TRUE
  )

compound_dict_mw_df <- compound_dict_mw %>%
  bind_rows() %>%
  transmute(
    lspci_id = as.integer(compound),
    mass
  )

write_csv(
  compound_dict_mw_df,
  here("wrangled", "compound_dict_mw.csv.gz")
)

cmap_mw <- unique_compounds_canonical %>%
  drop_na(inchi) %>%
  {set_names(.[["inchi"]], .[["pert_id"]])} %>%
  lspcheminf::molecular_mass()

compound_lspci_id_mapping_df <- compound_lspci_id_mapping %>%
  select(data, fingerprint_type) %>%
  unnest(data) %>%
  mutate_at(vars(target), as.integer) %>%
  left_join(
    select(compound_dict_mw_df, lspci_id, mass_lspci_id = mass),
    by = c("target" = "lspci_id")
  ) %>%
  left_join(
    select(cmap_mw, pert_id = compound, mass_pert_id = mass),
    by = c("query" = "pert_id")
  ) %>%
  # Remove matches with unqual mass
  filter(near(mass_lspci_id, mass_pert_id)) %>%
  group_by(query) %>%
  filter(
    n() == 2,
    all(c("morgan", "topological") %in% fingerprint_type),
    length(unique(target)) == 1
  ) %>%
  ungroup() %>%
  distinct(pert_id = query, lspci_id = target)

write_csv(
  compound_lspci_id_mapping_df,
  here("wrangled", "cmap_pertubation_compound_mapping.csv")
)

pertubation_meta <- pertubation_meta_raw %>%
  mutate(
    pert_iname_norm = str_replace_all(pert_iname, "[^a-zA-Z\\d]", "") %>% str_to_lower()
  ) %>%
  left_join(
    compound_lspci_id_mapping_df,
    by = "pert_id"
  )
  # left_join(pubchem_ids, by = "inchi_key") %>%
  # # Some compounds in the original dataset already have pubchem ids, merging
  # # them with the fetched ones
  # mutate(
  #   pubchem_cid = map2(
  #     pubchem_cid.x, pubchem_cid.y,
  #     ~unique(c(.x, .y)) %>%
  #       # Remove NA values if there is a non-NA value
  #       {if (sum(!is.na(.)) > 0) na.omit(.) else .}
  #   )
  # ) %>%
  # select(-pubchem_cid.x, -pubchem_cid.y) %>%
  # unnest(pubchem_cid)

# pertubation_meta %>% filter(!is.na(lspci_id))

write_csv(pertubation_meta, here("wrangled", "cmap_perturbation_meta.csv.gz"))
```



```{r instance_meta}
instance_meta_raw <- list(
  GSE92742 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_inst_info.txt.gz",
  GSE70138 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_inst_info_2017-03-06.txt.gz",
  LINCS_2020 = syn("syn25050165")
) %>%
  map(read_tsv, na = "-666") %>%
  bind_rows(.id = "dataset")

instance_meta <- instance_meta_raw %>%
  left_join(distinct(pertubation_meta, pert_id, pubchem_cid, lspci_id), by = "pert_id")

write_csv(instance_meta, here("wrangled", "cmap_instance_meta.csv.gz"))
```




```{r signature_meta}
signature_meta_raw <- list(
  GSE92742 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_sig_info.txt.gz",
  GSE70138 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_sig_info_2017-03-06.txt.gz",
  LINCS_2020 = syn("syn25050186")
) %>%
  map(read_tsv, na = c("-666", "-666.00")) %>%
  bind_rows(.id = "dataset")
# Parsing failures correspond to malformed fields in input data, are in irrelevant
# signatures, so ignore...

signature_meta <- signature_meta_raw %>%
  mutate(
    cell_id = fcoalesce(cell_id, cell_iname),
    pert_dose = if_else(is.na(pert_dose), as.double(str_match(pert_idose, "([\\d\\.]+) um")[, 2]), pert_dose) %>%
      {if_else(. < 0, NA_real_, .)},
    pert_dose_unit = if_else(str_detect(pert_idose, "([\\d\\.]+) um"), "µM", pert_dose_unit),
    pert_time = if_else(is.na(pert_time), as.double(str_match(pert_itime, "([\\d]+) h")[, 2]), pert_time),
    pert_time_unit = if_else(is.na(pert_itime), NA_character_, "h")
  ) %>%
  left_join(distinct(pertubation_meta, pert_id, lspci_id), by = "pert_id")

write_csv(signature_meta, here("wrangled", "cmap_signature_meta.csv.gz"))
```


```{r gene_meta}
gene_meta_raw <- list(
  GSE92742 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_gene_info.txt.gz",
  GSE70138 = "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_gene_info_2017-03-06.txt.gz",
  LINCS_2020 = syn("syn25050895")
) %>%
  map(read_tsv, na = "-666") %>%
  bind_rows(.id = "dataset")

# Check if the two datasets are identical
gene_meta_identical <- gene_meta_raw %>%
  group_by_at(vars(-dataset)) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  dplyr::count(n)
# They are!

gene_meta <- gene_meta_raw %>%
  filter(dataset == "GSE92742") %>%
  mutate_at(vars(pr_gene_id), as.character) %>%
  genebabel::join_hgnc(
    "pr_gene_id",
    "entrez_id",
    c("symbol", "entrez_id", "ensembl_gene_id")
  ) %>%
  select(-dataset)

write_csv(gene_meta, here("wrangled", "cmap_gene_meta.csv.gz"))
```

Gene metadata for both datasets are identical, can save just a single one.

```{r synapse}
activity <- synapser::Activity(
  name = "Wrangle CMap data",
  used = c(
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_pert_info.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_pert_info_2017-03-06.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_inst_info.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_inst_info_2017-03-06.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_sig_info.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_sig_info_2017-03-06.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92742/suppl/GSE92742_Broad_LINCS_gene_info.txt.gz",
    "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE70nnn/GSE70138/suppl/GSE70138_Broad_LINCS_gene_info_2017-03-06.txt.gz",
    "syn20835543",
    "syn21042105"
  ),
  executed = "https://github.com/clemenshug/dge_comparison/blob/master/wrangling/05_cmap_data.Rmd"
)

c(
  here("wrangled", "cmap_perturbation_meta.csv.gz"),
  here("wrangled", "cmap_instance_meta.csv.gz"),
  here("wrangled", "cmap_signature_meta.csv.gz"),
  here("wrangled", "cmap_gene_meta.csv.gz")
) %>%
  synStoreMany(parentId = "syn21547022", activity = activity)

#rsync -rv --exclude *.fq.gz --exclude fastq --exclude *.fastq.gz --exclude *.bam --exclude rapmap --exclude transcriptome 2018_02_steven_rodriguez_ad_repurposing transfer:/n/files/ImStor/sorger/data/rnaseq/lincs/dge_reanalysis_clemens/
```