README.Rmd

---
output: github_document
---

```{r, echo = FALSE, message = FALSE}
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  collapse = TRUE,
  echo = FALSE,
  comment = "#>",
  fig.path = "man/figures/",
  knitr::opts_chunk$set(dpi = 300, fig.width = 8)
)
```

# someta <img src="man/figures/logo.png" align="right" width="34%">
<!-- badges: start -->
[![R build status](https://github.com/rnabioco/someta/workflows/Query/badge.svg)](https://github.com/rnabioco/someta/actions)
[![Last Commit on GitHub](`r paste0("https://img.shields.io/badge/Last%20Run-", format(format = "%m--%d--%Y", Sys.Date()), "-brightgreen")`)](https://rnabioco.github.io/someta/articles/get_geo.html)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4686053.svg)](https://doi.org/10.5281/zenodo.4686053)
<!-- badges: end -->

Cell-type annotations are frequently excluded from public single cell datasets. This hinders single cell sequencing analysis reproducibility and accessibility. To better describe the issue, we monitor GEO entries monthly (currently set to auto-update at 1AM UTC, 1st of the month), and programmatically determine the fraction of entries with (potentially, likely overestimated) usable cell metadata. Descriptions of the issue and suggestions (in short: PLEASE PLEASE deposit some metadata at cell level for scRNA-seq data) are now published on [PLoS Biology](https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001077). Supplemental data files related to the manuscript are located in the `inst/manuscript/` directory.

Additional thoughtful guidelines for organizing scRNA-seq sample and cell metadata are discussed here by [Füllgrabe et al](https://www.nature.com/articles/s41587-020-00744-z).

The latest archive of GEO scRNA-seq records with other associated data can be directly downloaded here [current_geo.rds](https://github.com/rnabioco/someta/raw/master/inst/extdata/current_geo.rds). GEO filtering, preview first lines of text files, and content of .tar, and spot check feedback can be accessed from last tab on the `clustifyr` [web app](https://raysinensis.shinyapps.io/clustifyr-web-app/?tab=someta).

***

As of the initial presentation of this issue (10--17--2020), the number is a frustratingly low **0.122**.

```{r load}
library(tidyverse)
library(GEOquery)
library(here)
date <- format(format = "%m%d%y", Sys.Date())
date_1off <- format(format = "%m%d%y", Sys.Date() - 1)
if (file.exists(here("inst", "extdata", date, paste0("geo_", date, ".rds")))) {
  gds <- readRDS(here("inst", "extdata", date, paste0("geo_", date, ".rds"))) %>%
    mutate(year = as.numeric(year))
} else {
  # gds <- readRDS(here("inst", "extdata", date_1off, paste0("geo_", date_1off, ".rds"))) %>%
  #   mutate(year = as.numeric(year))
  gds <- readRDS(here("inst", "extdata", "current_geo.rds"))
}
gds_orig <- gds
```

```{r filteragain}
checks <- suppressWarnings(map(
  gds$geo,
  function(x) {
    str_detect(
      str_to_lower(x),
      "single nuclei|single-nucleus|single cell|single-cell|scrnaseq|scrna-seq|snrnaseq|snrna-seq"
    ) %>% any()
  }
))
gds <- gds[unlist(checks), ]
```

```{r superseries}
parse_subseries <- function(geolist) {
  tryCatch(
    {
      rel <- geolist$relation
      rel[str_detect(rel, "^SuperSeries of: ")] %>%
        str_remove("^SuperSeries of: ")
    },
    error = function(e) {
      return(NA)
    }
  )
}

res <- map(
  gds$geo,
  parse_subseries
)
names(res) <- gds$id
res2 <- res[lengths(res) != 0]
res2 <- map(
  res2,
  function(x) ifelse(x %in% gds$id, x, NA)
)
res3 <- data.frame(superseries = names(res2))
res3$id <- res2
dict_series <- res3 %>%
  unnest(id) %>%
  na.omit()
gds2 <- gds %>%
  left_join(dict_series) %>%
  mutate(superseries = ifelse(is.na(superseries),
    id,
    superseries
  )) %>%
  group_by(superseries) %>%
  summarize(
    has_meta1 = ifelse(any(has_meta1 == "yes"), "yes", "no"),
    has_meta2 = ifelse(any(has_meta2 == "yes"), "yes", "no"),
    has_object1 = ifelse(any(has_object1 == "yes"), "yes", "no"),
    has_object2 = ifelse(any(has_object2 == "yes"), "yes", "no")
  ) %>%
  mutate(usable = ifelse(has_meta1 == "yes" |
    has_meta2 == "yes" |
    has_object1 == "yes" |
    has_object2 == "yes",
  "yes", "no"
  )) %>%
  mutate(usable2 = ifelse(has_meta1 == "yes" | has_meta2 == "yes", "meta",
    ifelse(has_object1 == "yes" | has_object2 == "yes", "object", "no")
  ))

gds_fixed <- gds %>% inner_join(gds2 %>% select(id = superseries, fixed_usable = usable, usable2))
gds <- gds_fixed %>%
  select(-usable) %>%
  rename(usable = fixed_usable)
```

```{r fixfalsep}
gds_old <- gds
fp <- read_tsv(here("inst", "extdata", "falsepos.txt"), col_names = F) %>% pull(1)
gds <- gds %>% mutate(usable = ifelse(id %in% fp, "no", usable))
```

```{r ae}
sets <- read_tsv(here("inst", "extdata", "arrayexpress_sets_proc_add.tsv.gz"))

# remove GEO imported
sets2 <- sets %>% filter(!str_detect(ID, "GEO"))

sets_10x <- sets2 %>% filter(str_to_lower(str_sub(samplemetainfo_2, 1, 3)) == "10x")

sets_10x_withmeta <- sets_10x %>% 
  filter(
    str_detect(str_to_lower(samplemetainfo_1), "meta|annot|type|clustering|coldata") | 
    str_detect(str_to_lower(additionalfiles), "meta|annot|type|clustering|coldata")
  )

sets_final <- sets_10x %>% 
  mutate(has_meta = ifelse(
    str_detect(str_to_lower(samplemetainfo_1), "meta|annot|type|clustering|coldata") | 
    str_detect(str_to_lower(additionalfiles), "meta|annot|type|clustering|coldata"),
    "yes",
    "no")) %>% 
  mutate(has_meta = ifelse(
    is.na(has_meta), "no", has_meta
  )) %>% 
  select(ID, ReleaseDate, PubmedID, Species, 
         files = samplemetainfo_1, type = samplemetainfo_2,
         has_meta) %>% 
  filter(str_sub(ReleaseDate, 1, 4) != "2021")

# write_tsv(sets_final, "arrayexpres_analysis_2020.tsv.gz")
```

Current fraction in GEO with metadata: **`r (gds$usable == "yes") %>% mean(na.rm = TRUE)`**. In comparison, for ArrayExpress 10x datasets, the fraction is **`r nrow(sets_10x_withmeta) / nrow(sets_10x)`**).

```{r comp}
gds_comp <- readRDS(here("inst", "extdata", "101520", paste0("geo_", "101520", ".rds"))) %>%
  mutate(year = as.numeric(year))
id_comp <- gds_comp %>%
  filter(usable == "no") %>%
  pull(id)
id_gds <- gds_comp %>%
  filter(usable == "yes") %>%
  pull(id)
```

Number of depositions with updated metadata records since description of the issue here: **`r length(intersect(id_gds, id_comp))`**.

Please also see [full report page](https://rnabioco.github.io/someta/articles/get_geo.html).

```{r frac}
frac <- gds %>%
  na.omit() %>%
  group_by(year, usable) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  pivot_wider(names_from = usable, values_from = n) %>%
  replace_na(list(no = 0, yes = 0)) %>%
  mutate(fraction = yes / (yes + no), total = yes + no) %>%
  mutate(fraction = round(fraction, digits = 3))

frac2 <- gds %>%
  na.omit() %>%
  group_by(year, usable2) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  pivot_wider(names_from = usable2, values_from = n) %>%
  replace_na(list(no = 0, meta = 0, object = 0)) %>%
  mutate(fraction = (meta + object) / (meta + object + no), total = meta + object + no) %>%
  mutate(fraction = round(fraction, digits = 3))

g <- ggplot(
  gds %>% na.omit() %>%
    filter(year >= 2012) %>%
    mutate(usable2 = factor(usable2, levels = c("no", "meta", "object"))),
  aes(x = year)
) +
  geom_bar(aes(fill = usable2)) +
  scale_fill_brewer(name = "usable") +
  geom_text(data = frac2 %>% filter(year >= 2012), color = "black", aes(x = year, y = total + max(frac$total) * 0.03, label = fraction), size = 3) +
  cowplot::theme_cowplot() +
  scale_x_continuous(breaks = c(2012:format(Sys.time(), "%Y") %>% as.numeric())) +
  ylab("total studies") +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 14),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)
  )
g + ggtitle("scRNAseq datasets by year")

last12 <- gds %>%
  na.omit() %>%
  select(year, month) %>%
  distinct() %>%
  mutate(month = factor(month, levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))) %>%
  arrange(desc(year), desc(as.numeric(month))) %>%
  dplyr::slice(1:12) %>%
  mutate(full = str_c(year, month, sep = "-")) %>%
  pull(full)

last12_df <- gds %>%
  na.omit() %>%
  mutate(full = str_c(year, month, sep = "-")) %>%
  filter(full %in% last12) %>%
  mutate(month = factor(month, levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))) %>%
  arrange(year, as.numeric(month))

last12_frac <- last12_df %>%
  na.omit() %>%
  group_by(year, month, usable2) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  pivot_wider(names_from = usable2, values_from = n) %>%
  replace_na(list(no = 0, meta = 0, object = 0)) %>%
  mutate(fraction = (meta + object) / (meta + object + no), total = meta + object + no) %>%
  mutate(fraction = round(fraction, digits = 3)) %>%
  mutate(full = str_c(year, month, sep = "-"))

g2 <- ggplot(
  last12_df %>% mutate(full = factor(full, levels = unique(last12_df$full))) %>%
    mutate(usable2 = factor(usable2, levels = c("no", "meta", "object"))),
  aes(x = full)
) +
  geom_bar(aes(fill = usable2)) +
  geom_text(
    data = last12_frac %>% mutate(full = factor(full, levels = unique(last12_df$full))),
    color = "black",
    aes(x = full, y = total + max(last12_frac$total) * 0.03, label = fraction), size = 3
  ) +
  scale_fill_brewer(name = "usable") +
  cowplot::theme_cowplot() +
  ylab("total studies") +
  xlab("last 12 months") +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 14),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)
  )
g2 + ggtitle("scRNAseq datasets by month")
```

***