library(LBDiscover)
#> Loading LBDiscover package

# Define the primary term of interest for our analysis
primary_term <- "migraine"

# Search for migraine-related articles
migraine_articles <- pubmed_search(
  query = paste0(primary_term, " pathophysiology"),
  max_results = 1000
)

# Search for treatment-related articles
drug_articles <- pubmed_search(
  query = "neurological drugs pain treatment OR migraine therapy OR headache medication",
  max_results = 1000
)

# Combine and remove duplicates
all_articles <- merge_results(migraine_articles, drug_articles)
cat("Retrieved", nrow(all_articles), "unique articles\n")
#> Retrieved 1929 unique articles

# Extract variations of our primary term using the utility function
primary_term_variations <- get_term_vars(all_articles, primary_term)
cat("Found", length(primary_term_variations), "variations of", primary_term, "in the corpus:\n")
#> Found 19 variations of migraine in the corpus:
print(head(primary_term_variations, 10))
#>  [1] "migraine"    "Migraine"    "MigrainE"    "MIGRAINE"    "migraines"  
#>  [6] "Migraines"   "migraineur"  "migraineurs" "Migraineurs" "AimMigraine"

# Preprocess text
preprocessed_articles <- preprocess_text(
  all_articles,
  text_column = "abstract",
  remove_stopwords = TRUE,
  min_word_length = 2  # Set min_word_length to capture short terms
)
#> Tokenizing text...

# Create a custom dictionary with all variations of our primary term
custom_dictionary <- data.frame(
  term = c(primary_term, primary_term_variations),
  type = rep("disease", length(primary_term_variations) + 1),
  id = paste0("CUSTOM_", 1:(length(primary_term_variations) + 1)),
  source = rep("custom", length(primary_term_variations) + 1),
  stringsAsFactors = FALSE
)

# Define additional MeSH queries for extended dictionaries
mesh_queries <- list(
  "disease" = paste0(primary_term, " disorders[MeSH] OR headache disorders[MeSH]"),
  "protein" = "receptors[MeSH] OR ion channels[MeSH]",
  "chemical" = "neurotransmitters[MeSH] OR vasoactive agents[MeSH]",
  "pathway" = "signal transduction[MeSH] OR pain[MeSH]",
  "drug" = "analgesics[MeSH] OR serotonin agonists[MeSH] OR anticonvulsants[MeSH]",
  "gene" = "genes[MeSH] OR channelopathy[MeSH]"
)

# Sanitize the custom dictionary
custom_dictionary <- sanitize_dictionary(
  custom_dictionary,
  term_column = "term",
  type_column = "type",
  validate_types = FALSE  # Don't validate custom terms as they're trusted
)
#> Sanitizing dictionary with 20 terms...
#> Sanitization complete. 20 terms remaining (100% of original)

# Extract entities using our custom dictionary
custom_entities <- extract_entities(
  preprocessed_articles,
  text_column = "abstract",
  dictionary = custom_dictionary,
  case_sensitive = FALSE,
  overlap_strategy = "priority",
  sanitize_dict = FALSE  # Already sanitized
)

# Extract entities using the standard workflow with improved entity validation
# Check if running in R CMD check environment
is_check <- !interactive() && 
            (!is.null(Sys.getenv("R_CHECK_RUNNING")) && 
             Sys.getenv("R_CHECK_RUNNING") == "true")
             
# More robust check for testing environment
if (!is_check && !is.null(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) {
  is_check <- TRUE
}

# Set number of cores based on environment
num_cores_to_use <- if(is_check) 1 else 4

standard_entities <- extract_entities_workflow(
  preprocessed_articles,
  text_column = "abstract",
  entity_types = c("disease", "drug", "gene"),
  parallel = !is_check,           # Disable parallel in check environment
  num_cores = num_cores_to_use,   # Use 1 core in check environment
  batch_size = 500                # Process 500 documents per batch
)
#> Warning in extract_entities_workflow(preprocessed_articles, text_column =
#> "abstract", : UMLS source requested but no API key provided. Skipping UMLS.

# Uncomment to include UMLS entities
# standard_entities <- extract_entities_workflow(
#   preprocessed_articles,
#   text_column = "abstract",
#   entity_types = c("disease", "drug", "gene", "protein", "pathway", "chemical"),
#   dictionary_sources = c("local", "mesh", "umls"),  # Including UMLS
#   additional_mesh_queries = mesh_queries,
#   sanitize = TRUE,
#   api_key = "your-umls-api-key",  # Your UMLS API key here
#   parallel = TRUE,
#   num_cores = 4
# )

# Combine entity datasets using our utility function
entities <- merge_entities(
  custom_entities,
  standard_entities,
  primary_term
)
#> Combined 6295 custom entities with 6101 standard entities.

# Filter entities to ensure only relevant biomedical terms are included
filtered_entities <- valid_entities(
  entities,
  primary_term,
  primary_term_variations,
  validation_function = is_valid_biomedical_entity
)
#> Filtered from 6358 to 6358 validated entities

# View the first few extracted entities
head(filtered_entities)
#>         entity entity_type doc_id start_pos end_pos
#> 1  AimMigraine     disease   1805         1      11
#> 2 antimigraine     disease    207       596     607
#> 3 antimigraine     disease    207       985     996
#> 4 antimigraine     disease    767      1488    1499
#> 5 antimigraine     disease    774      1671    1682
#> 6 antimigraine     disease   1505       948     959
#>                                                                                                                                                                                                            sentence
#> 1                               AimMigraine substantially affects daily functioning, yet its impact on unpaid domestic labor, an essential and sex-structured component of home life, remains poorly characterized.
#> 2            The aim of this review is to highlight the neurovascular mechanisms underlying poor response or nonresponse to antimigraine drugs, with particular focus on the phenomena related to the CGRP pathway.
#> 3                                                               RESULTS: Non-response to antimigraine treatment is closely linked to clinical and neurobiological factors that increase the overall disease burden.
#> 4 CONCLUSION AND IMPLICATIONS: DOP-mediated antimigraine effect is stronger in female than male rats and appears to be mediated in chronic conditions through peripheral DOPs in females and central DOPs in males.
#> 5                                            Our results suggest that the antagonism on the alpha-CGRP isoform plays a relevant role in the antimigraine action of mABs but not in the development of constipation.
#> 6                                                                   This review discusses the efficacy of nanocarriers loaded with antimigraine agents and delivered by the IN route in the management of migraine.
#>   frequency
#> 1         1
#> 2         7
#> 3         7
#> 4         7
#> 5         7
#> 6         7

# Create co-occurrence matrix with validated entities
co_matrix <- create_comat(
  filtered_entities,
  doc_id_col = "doc_id",
  entity_col = "entity",
  type_col = "entity_type",
  normalize = TRUE,
  normalization_method = "cosine"
)
#> Building entity-document matrix...
#> Calculating co-occurrence matrix...
#> Normalizing co-occurrence matrix using cosine method...

# Find our primary term in the co-occurrence matrix
a_term <- find_term(co_matrix, primary_term)
#> Found primary term in co-occurrence matrix

# Check matrix dimensions
dim(co_matrix)
#> [1] 15 15

# Apply the improved ABC model with enhanced term filtering and type validation
abc_results <- abc_model(
  co_matrix,
  a_term = a_term,
  c_term = NULL,  # Allow all potential C terms
  min_score = 0.001,  # Lower threshold to capture more potential connections
  n_results = 500,    # Increase to get more candidates before filtering
  scoring_method = "combined",
  # Focus on biomedically relevant entity types
  b_term_types = c("protein", "gene", "pathway", "chemical"),
  c_term_types = c("drug", "chemical", "protein", "gene"),
  exclude_general_terms = TRUE,  # Enable enhanced term filtering
  filter_similar_terms = TRUE,   # Remove terms too similar to migraine
  similarity_threshold = 0.7,    # Relatively strict similarity threshold
  enforce_strict_typing = TRUE   # Enable strict entity type validation
)
#> Filtered 13 B terms (92.9%) that weren't valid biomedical entities
#> Warning: Term filtering removed most B terms (92.9% filtered). Results may be limited.
#> Filtered 1 B terms that didn't match specified entity types: protein, gene, pathway, chemical
#> No suitable B terms found with association score > 0.001 after filtering

# If we don't have enough results, try with less stringent criteria
min_desired_results <- 10
if (nrow(abc_results) < min_desired_results) {
  cat("Not enough results with strict filtering. Trying with less stringent criteria...\n")
  abc_results <- abc_model(
    co_matrix,
    a_term = a_term,
    c_term = NULL,
    min_score = 0.0005,  # Even lower threshold
    n_results = 500,
    scoring_method = "combined",
    b_term_types = NULL,  # No type constraints
    c_term_types = NULL,  # No type constraints
    exclude_general_terms = TRUE,
    filter_similar_terms = TRUE,
    similarity_threshold = 0.8,    # More lenient similarity threshold
    enforce_strict_typing = FALSE  # Disable strict type validation as fallback
  )
}
#> Not enough results with strict filtering. Trying with less stringent criteria...
#> Filtered 13 B terms (92.9%) that weren't valid biomedical entities
#> Warning: Term filtering removed most B terms (92.9% filtered). Results may be limited.
#> Filtered out 0 B terms that were too similar to A term (similarity threshold: 0.8)
#> Filtered 13 potential C terms that weren't valid biomedical entities
#> Identifying potential C terms via 1 B terms...
#>   |                                                                              |                                                                      |   0%  |                                                                              |======================================================================| 100%
#> No ABC connections found

# Check if we have results and what columns are available
if (nrow(abc_results) > 0) {
  cat("Available columns in abc_results:", paste(colnames(abc_results), collapse = ", "), "\n")
  
  # Safely display results with available columns
  display_cols <- intersect(c("a_term", "b_term", "c_term", "abc_score", "b_type", "c_type"), 
                           colnames(abc_results))
  
  if (length(display_cols) > 0) {
    # View top results
    head(abc_results[, display_cols, drop = FALSE])
  } else {
    # If specific columns aren't available, show first few columns
    head(abc_results[, 1:min(6, ncol(abc_results)), drop = FALSE])
  }
} else {
  cat("No ABC results found\n")
}
#> No ABC results found

# Apply statistical validation to the results
validated_results <- tryCatch({
  validate_abc(
    abc_results,
    co_matrix,
    alpha = 0.1,  # More lenient significance threshold
    correction = "BH",  # Benjamini-Hochberg correction for multiple testing
    filter_by_significance = FALSE  # Keep all results but mark significant ones
  )
}, error = function(e) {
  cat("Error in statistical validation:", e$message, "\n")
  cat("Using original results without validation...\n")
  # Add dummy p-values based on ABC scores
  abc_results$p_value <- 1 - abc_results$abc_score / max(abc_results$abc_score, na.rm = TRUE)
  abc_results$significant <- abc_results$p_value < 0.1
  return(abc_results)
})
#> ABC results are empty

# Sort by ABC score and take top results
validated_results <- validated_results[order(-validated_results$abc_score), ]
top_n <- min(100, nrow(validated_results))  # Larger top N for diversification
top_results <- head(validated_results, top_n)

# Check available columns and display safely
if (nrow(top_results) > 0) {
  display_cols <- intersect(c("a_term", "b_term", "c_term", "abc_score", "p_value", "significant"), 
                           colnames(top_results))
  
  if (length(display_cols) > 0) {
    # View top validated results
    head(top_results[, display_cols, drop = FALSE])
  } else {
    # Show available columns
    head(top_results[, 1:min(6, ncol(top_results)), drop = FALSE])
  }
} else {
  cat("No validated results available\n")
}
#> No validated results available

# Diversify results using our utility function
diverse_results <- safe_diversify(
  top_results,
  diversity_method = "both",
  max_per_group = 5,
  min_score = 0.0001,
  min_results = 5
)
#> ABC results are empty

# Ensure we have enough results for visualization
diverse_results <- min_results(
  diverse_results,
  top_results,
  a_term,
  min_results = 3
)
#> No results found. Creating a placeholder result for demonstration.

# Create heatmap visualization
plot_heatmap(
  diverse_results,
  output_file = "migraine_heatmap.png",
  width = 1200,
  height = 900,
  top_n = 15,
  min_score = 0.0001,
  color_palette = "blues",
  show_entity_types = TRUE
)
#> Warning in vis_heatmap(results, top_n = min(top_n, nrow(results)), min_score =
#> min_score, : Entity types not found in results. Setting show_entity_types =
#> FALSE
#> No connections are statistically significant (p < 0.05)
#> Created heatmap visualization: migraine_heatmap.png

# Create network visualization
plot_network(
  diverse_results,
  output_file = "migraine_network.png",
  width = 1200,
  height = 900,
  top_n = 15,
  min_score = 0.0001,
  node_size_factor = 5,
  color_by = "type",
  title = "Migraine Treatment Network",
  show_entity_types = TRUE,
  label_size = 1.0
)
#> Warning in vis_network(results, top_n = min(top_n, nrow(results)), min_score =
#> min_score, : Entity types not found in results. Setting show_entity_types =
#> FALSE
#> Created network visualization: migraine_network.png

# Create interactive HTML network visualization
export_network(
  diverse_results,
  output_file = "migraine_network.html",
  top_n = min(30, nrow(diverse_results)),
  min_score = 0.0001,
  open = FALSE  # Don't automatically open in browser
)

# Create interactive chord diagram
export_chord(
  diverse_results,
  output_file = "migraine_chord.html",
  top_n = min(30, nrow(diverse_results)),
  min_score = 0.0001,
  open = FALSE
)
#> Number of unique terms: 7
#> First few terms: migraine, serotonin, CGRP, cortical spreading depression, sumatriptan
#> Role assignments: A=1, B=3, C=3

# Evaluate literature support for top connections
evaluation <- eval_evidence(
  diverse_results,
  max_results = 5,
  base_term = "migraine",
  max_articles = 5
)
#> 
#> === Evaluation of Top Results ===
#> 
#> Evaluating potential treatment: sumatriptan  
#> ABC score: 0.04 
#> P-value: 0.1 - Not statistically significant 
#> Connection through intermediary: serotonin
#> Found 5 articles directly linking migraine and sumatriptan 
#> Most recent article: Ibogalogs reduce acute and chronic headache-like behaviours in rodents by different co-modulatory mechanisms involving 5-HT1 and 5-HT2 receptor subtypes. 
#> 
#> Evaluating potential treatment: topiramate  
#> ABC score: 0.03 
#> P-value: 0.2 - Not statistically significant 
#> Connection through intermediary: CGRP
#> Found 5 articles directly linking migraine and topiramate 
#> Most recent article: More Than Meets the Eye: Topiramate-Induced Bilateral Secondary Angle Closure Unmasking Secondary Antiphospholipid Syndrome. 
#> 
#> Evaluating potential treatment: propranolol  
#> ABC score: 0.02 
#> P-value: 0.3 - Not statistically significant 
#> Connection through intermediary: cortical spreading depression
#> Found 5 articles directly linking migraine and propranolol 
#> Most recent article: Comparison table: Some drugs for migraine prevention.

# Prepare articles for report generation
articles_with_years <- prep_articles(all_articles)
#> Found 1929 articles with valid publication years

# Store results for report
results_list <- list(abc = diverse_results)

# Store visualization paths
visualizations <- list(
  heatmap = "migraine_heatmap.png",
  network = "migraine_network.html",
  chord = "migraine_chord.html"
)

# Create comprehensive report
gen_report(
  results_list = results_list,
  visualizations = visualizations,
  articles = articles_with_years,
  output_file = "migraine_discoveries.html"
)
#> Generated comprehensive report: migraine_discoveries.html

cat("\nDiscovery analysis complete!\n")
#> 
#> Discovery analysis complete!