[WIP] Option for isolated and colony profiles #30
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -3,56 +3,78 @@ | |
| 'aggregate | ||
|
|
||
| Usage: | ||
| aggregate.R -s <sqlite_file> -o <file> [-t <sc_type>] | ||
|
|
||
| Options: | ||
| -h --help Show this screen. | ||
| -s <sqlite_file> --sqlite_file=<sqlite_file> Location of the sql_lite file. | ||
| -o <file> --output=<file> per-well aggregated data. | ||
| -t <sc_type> --sc_type=<sc_type> which sc_type to focus on [default: none]' -> doc | ||
|
|
||
| suppressWarnings(suppressMessages(library(docopt))) | ||
| suppressWarnings(suppressMessages(library(dplyr))) | ||
| suppressWarnings(suppressMessages(library(magrittr))) | ||
| suppressWarnings(suppressMessages(library(stringr))) | ||
|
|
||
| opts <- docopt(doc) | ||
|
|
||
| sql_file <- opts[["sqlite_file"]] | ||
| output_file <- opts[["output"]] | ||
| sc_type <- opts[["sc_type"]] | ||
|
|
||
| db <- src_sqlite(path = sql_file) | ||
|
|
||
| image <- tbl(src = db, "image") %>% | ||
| select(TableNumber, ImageNumber, Metadata_Plate, Metadata_Well) | ||
|
|
||
| aggregate_objects <- function(compartment, sc_type = "none") { | ||
| object <- tbl(src = db, compartment) | ||
|
|
||
| object %<>% inner_join(image, by = c("TableNumber", "ImageNumber")) | ||
|
|
||
| if (sc_type == "isolated") { | ||
| object <- object %>% dplyr::filter(Cells_Neighbors_NumberOfNeighbors_Adjacent == 0) | ||
|
There was a problem hiding this comment. @gwaygenomics This wouldn't work when There was a problem hiding this comment. I originally wrote this (see below), but then realized I am wrong. This may have been my original plan, but I do not remember 100% you are correct, this will fail as currently implemented! This is because we're aggregating compartment by compartment. If we merge everything first and then aggregate, it should work. It seems like some version of this must have worked at one point (because the call on L72 would have errored out for other non-cell compartments). Original thought that is wrongIf I'm remembering correctly, Since we're doing an inner join on L72, this should take care of filtering for other compartments |
||
| } else if (sc_type == "colony") { | ||
| object <- object %>% dplyr::filter(Cells_Neighbors_NumberOfNeighbors_Adjacent >= 4) | ||
| } | ||
|
|
||
| # compartment tag converts nuclei to ^Nuclei_ | ||
| compartment_tag <- | ||
| str_c("^", str_sub(compartment, 1, 1) %>% str_to_upper(), str_sub(compartment, 2), "_") | ||
|
|
||
| variables <- colnames(object) %>% stringr::str_subset(compartment_tag) | ||
| futile.logger::flog.info(str_c("Started aggregating ", compartment)) | ||
|
|
||
| object <- object %>% | ||
|
There was a problem hiding this comment. This bit here is also specific to aggregating single cell profiles It takes code from |
||
| dplyr::as_tibble() %>% | ||
| dplyr::group_by_(.dots = c("TableNumber", "ImageNumber", | ||
| "ObjectNumber", "Metadata_Plate", | ||
| "Metadata_Well")) %>% | ||
| dplyr::summarise_at(.funs = 'mean', .vars = variables) %>% | ||
| dplyr::ungroup() | ||
|
|
||
| return(object) | ||
| } | ||
|
|
||
| aggregate_cols <- c("TableNumber", "ImageNumber", "ObjectNumber", | ||
|
There was a problem hiding this comment. This is also specific to single cells and using it in other repositories will not always work |
||
| "Metadata_Plate", "Metadata_Well") | ||
| sc_objects <- | ||
| aggregate_objects("cells", | ||
| sc_type = sc_type) %>% | ||
| inner_join(aggregate_objects("cytoplasm"), | ||
| by = aggregate_cols) %>% | ||
| inner_join(aggregate_objects("nuclei"), | ||
| by = aggregate_cols) | ||
|
|
||
| futile.logger::flog.info(paste0("Now collapsing by well for ", output_file)) | ||
|
|
||
| variables <- colnames(sc_objects) %>% stringr::str_subset("^Cells|^Nuclei|^Cytoplasm") | ||
|
|
||
| sc_objects <- sc_objects %>% | ||
| dplyr::group_by_(.dots = c("Metadata_Plate", "Metadata_Well")) %>% | ||
| dplyr::summarize_at(.funs = 'mean', .vars = variables) %>% | ||
| dplyr::ungroup() | ||
|
|
||
| futile.logger::flog.info(paste0("Writing aggregated to ", output_file)) | ||
|
|
||
| sc_objects %>% readr::write_csv(output_file) | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -24,21 +24,15 @@ suppressWarnings(suppressMessages(library(readr))) | |
| opts <- docopt(doc) | ||
|
|
||
| batch_id <- opts[["batch_id"]] | ||
| external_metadata <- opts[["external_metadata"]] | ||
| cell_id <- opts[["cell_id"]] | ||
| format_broad_cmap <- opts[["format_broad_cmap"]] | ||
| plate_id <- opts[["plate_id"]] | ||
| perturbation_mode <- opts[["perturbation_mode"]] | ||
| metadata_dir <- paste("../..", "metadata", batch_id, sep = "/") | ||
| backend_dir <- paste("../..", "backend", batch_id, unlist(strsplit(plate_id, "_"))[1], sep = "/") | ||
|
There was a problem hiding this comment. @gwaygenomics Do you know why the string is being split? There was a problem hiding this comment. As we discussed on the brief chat, indeed, it looks like I am splitting the plate ID because I am inputing the plate ID with the This is a result of the output from the modified This will definitely break in situations where the plate ID has an underscore. Will need to address! |
||
|
|
||
| print(backend_dir) | ||
| # read profiles and rename column names | ||
|
|
||
| profiles <- suppressMessages(readr::read_csv(paste(backend_dir, paste0(plate_id, ".csv"), sep = "/"))) | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
after reviewing, I don't think this change is ready for prime time. A couple more conditional statements that only behave when
sc_typeis specified are required.