dropped one pf8 study that overlaps with PMID

analysis/10_read_pf8k.R

@@ -76,6 +76,11 @@ mdf <- mdf |>
   mutate(country = ifelse(country == "Côte d'Ivoire", "Cote d'Ivoire", country)) |>
   filter(country %in% pf8_Africa)
 
+# drop one study that is duplicated with published paper (PMID 27899115
+# corresponding to MalariaGen study 1248-PF-MW-MCCANN-SM)
+mdf <- mdf |>
+  filter(study_id != "MalariaGEN_1248_PF_MW_MCCANN_SM")
+
 # count studies
 length(unique(mdf$study_id))
 

analysis/11_merge_stave.R

@@ -51,7 +51,7 @@ s$append_data(studies_dataframe = s_pf8$get_studies(),
               counts_dataframe = s_pf8$get_counts())
 
 # save combined object to file
-saveRDS(s, file = here("analysis", "data-out", "stave_data_2026.03.13.rds"))
+saveRDS(s, file = here("analysis", "data-out", "stave_data_2026.03.17.rds"))
 
 # ------------------------------------------------------------------
 # Some basic stats on the final object

analysis/13_deduplicate_pf8.R

@@ -0,0 +1,95 @@
+
+# ------------------------------------------------------------
+# get abstracts from PMID for all non-pf8 studies
+# (NB, if complains because too many queries, break PMID vector into multiple
+# sections and stick back together)
+
+library(rentrez)
+library(xml2)
+
+s_geoff <- readRDS(file = here("analysis", "data-derived", "geoff_STAVE.rds"))
+s_wwarn <- readRDS(file = here("analysis", "data-derived", "WWARN_STAVE.rds"))
+s_who <- readRDS(file = here("analysis", "data-derived", "WHO_STAVE.rds"))
+
+# Input PMIDs
+pmids <- c(s_geoff$get_studies()$PMID,
+           s_wwarn$get_studies()$PMID,
+           s_who$get_studies()$PMID)
+
+pmids <- pmids[!is.na(pmids)]
+
+# Fetch PubMed records as XML
+xml_text_raw <- entrez_fetch(
+  db = "pubmed",
+  id = pmids,
+  rettype = "xml",
+  retmode = "xml"
+)
+
+# Parse XML
+doc <- read_xml(xml_text_raw)
+
+# Find all PubMed articles
+articles <- xml_find_all(doc, ".//PubmedArticle")
+
+# Extract fields
+results <- lapply(articles, function(article) {
+
+  pmid_node <- xml_find_first(article, ".//PMID")
+  title_node <- xml_find_first(article, ".//ArticleTitle")
+  abstract_nodes <- xml_find_all(article, ".//Abstract/AbstractText")
+  journal_node <- xml_find_first(article, ".//Journal/Title")
+  year_node <- xml_find_first(article, ".//PubDate/Year")
+
+  pmid <- if (inherits(pmid_node, "xml_missing")) NA_character_ else xml_text(pmid_node)
+  title <- if (inherits(title_node, "xml_missing")) NA_character_ else xml_text(title_node)
+  journal <- if (inherits(journal_node, "xml_missing")) NA_character_ else xml_text(journal_node)
+  year <- if (inherits(year_node, "xml_missing")) NA_character_ else xml_text(year_node)
+
+  abstract <- if (length(abstract_nodes) == 0) {
+    NA_character_
+  } else {
+    paste(xml_text(abstract_nodes), collapse = " ")
+  }
+
+  data.frame(
+    PMID = pmid,
+    Title = title,
+    Abstract = abstract,
+    Journal = journal,
+    Year = year,
+    stringsAsFactors = FALSE
+  )
+})
+
+# Combine into one data frame
+pubmed_data <- do.call(rbind, results)
+
+# ------------------------------------------------------------
+# search for WGS data 
+
+# option to read in pre-computed from file
+#pubmed <- readRDS("/Users/rverity/Desktop/pubmed.rds")
+
+# Patterns related to whole genome sequencing
+wgs_patterns <- c(
+  "\\bwhole genomes? sequenc\\w*\\b",
+  "\\bwhole-genome sequenc\\w*\\b",
+  "\\bsequencing of the whole genome\\b",
+  "\\bwhole genome data\\b",
+  "\\bwgs\\b",
+  "\\bgenome[- ]wide sequencing\\b"
+)
+
+# match abstracts against patterns
+pubmed$WGS <- FALSE
+for (i in 1:nrow(pubmed)) {
+  pubmed$WGS[i] <- any(str_detect(tolower(pubmed$Abstract[[i]]), wgs_patterns))
+}
+
+# manually fix papers where cannot get the abstract automatically
+pubmed$WGS[346] <- FALSE
+
+pubmed |>
+  filter(WGS) |>
+  View()