Refactor custom probe annotation process

- Update Agilent catalog design handling to pull only ProbeID to Ensembl
  Transcript ID mapping. Map to Ensembl Gene IDs using BioMart and then
  proceed to standard GeneLab Annotation process
- Update outdated links
- Clean up formatting and typos in all documentation

Author: bnovak32

Microarray/Agilent_1-channel/Pipeline_GL-DPPD-7112_Versions/GL-DPPD-7112-A.md

@@ -1117,9 +1117,7 @@ if (organism %in% c("athaliana")) {
     }
   }
 
-  if (use_custom_annot) {
-    unloadNamespace("biomaRt")
-  } else {
+  if (!use_custom_annot) {
     print(ensembl)
 
     probe_ids <- unique(norm_data$genes$ProbeName)
@@ -1148,66 +1146,85 @@ if (organism %in% c("athaliana")) {
 
       Sys.sleep(10) # Slight break between requests to prevent back-to-back requests
     }
-  }
-}
-
-# At this point, we have df_mapping from either the BioMart live service or the ensembl genomes ftp archive depending on the organism
-# If no df_mapping obtained (e.g., not supported in BioMart), use custom annotations; otherwise, merge in-house annotations to df_mapping
-
-if (use_custom_annot) {
-  expected_attribute_name <- 'ProbeName'
-
-  annot_type <- 'NO_CUSTOM_ANNOT'
-  if (!is.null(local_annotation_dir) && !is.null(annotation_config_path)) {
-    config_df <- read.csv(annotation_config_path, row.names=1)
-    if (unique(df_rs$`biomart_attribute`) %in% row.names(config_df)) {
-      annot_config <- config_df[unique(df_rs$`biomart_attribute`), ]
-      annot_type <- annot_config$annot_type[[1]]
+  } else {
+    print(glue::glue("Using custom probe annotation for organism: {unique(df_rs$organism)}"))
+    expected_attribute_name <- 'ProbeName'
+
+    annot_type <- 'NO_CUSTOM_ANNOT'
+    if (!is.null(local_annotation_dir) && !is.null(annotation_config_path)) {
+      config_df <- read.csv(annotation_config_path, row.names=1)
+      if (unique(df_rs$`biomart_attribute`) %in% row.names(config_df)) {
+        annot_config <- config_df[unique(df_rs$`biomart_attribute`), ]
+        annot_type <- annot_config$annot_type[[1]]
+      } else {
+        warning(paste0("No entry for '", unique(df_rs$`biomart_attribute`), "' in provided config file: ", annotation_config_path))
+      }
     } else {
-      warning(paste0("No entry for '", unique(df_rs$`biomart_attribute`), "' in provided config file: ", annotation_config_path))
+      warning("Need to provide both local_annotation_dir and annotation_config_path to use custom annotation.")
     }
-  } else {
-    warning("Need to provide both local_annotation_dir and annotation_config_path to use custom annotation.")
-  }
-
-  if (annot_type == 'agilent') {
-    unique_probe_ids <- read.delim(
-      file.path(local_annotation_dir, annot_config$annot_filename[[1]]),
-      header = TRUE, na.strings = c('NA', '')
-    )[c('ProbeID', 'EnsemblID', 'GeneSymbol', 'GeneName', 'RefSeqAccession', 'EntrezGeneID', 'GO')]
 
-    stopifnot(nrow(unique_probe_ids) == length(unique(unique_probe_ids$ProbeID)))
+    if (annot_type == 'agilent') {
+      print(glue::glue("Using Agilent Ensembl transcript annotation to retrieve Ensembl Gene IDs from BioMart"))
+      print(ensembl)
+
+      # read in AllAnnotations file and strip version off of Ensembl Transcript ID
+      agilent_table <- read.delim(
+          file.path(local_annotation_dir, annot_config$annot_filename[[1]]),
+          header = TRUE, na.strings = c('NA', '')
+        )[c('ProbeID', 'EnsemblID')] %>%
+        dplyr::mutate(EnsemblID = stringr::str_split_i(EnsemblID, "\\.", 1))
+
+      stopifnot(nrow(agilent_table) == length(unique(agilent_table$ProbeID)))
+
+      # Map Ensembl Transcript IDs to Ensembl Gene IDs using BioMart
+      # Run BioMart Queries in chunks to prevent request timeouts
+      #   Note: If timeout is occurring (possibly due to larger load on BioMart), reduce chunk size
+      CHUNK_SIZE= 1500
+      transcript_id_chunks <- split(agilent_table$EnsemblID, ceiling(seq_along(agilent_table$EnsemblID) / CHUNK_SIZE))
+      id_map <- data.frame()
+      for (i in seq_along(transcript_id_chunks)) {
+        transcript_id_chunk <- transcript_id_chunks[[i]]
+        print(glue::glue("Running BioMart query chunk {i} of {length(transcript_id_chunks)}. Total transcript IDS in query ({length(transcript_id_chunk)})"))
+        chunk_results <- biomaRt::getBM(
+            attributes = c(
+                "ensembl_transcript_id",
+                "ensembl_gene_id"
+                ), 
+                filters = "ensembl_transcript_id", 
+                values = transcript_id_chunk, 
+                mart = ensembl)
+
+        if (nrow(chunk_results) > 0) {
+          id_map <- id_map %>% dplyr::bind_rows(chunk_results)
+        }
 
-    # Clean columns
-    unique_probe_ids$GO <- purrr::map_chr(stringr::str_extract_all(unique_probe_ids$GO, 'GO:\\d{7}'), ~paste0(unique(.), collapse = '|')) %>% stringr::str_replace('^$', NA_character_)
+        Sys.sleep(10) # Slight break between requests to prevent back-to-back requests
+      }
 
-    names(unique_probe_ids) <- c('ProbeName', 'ENSEMBL', 'SYMBOL', 'GENENAME', 'REFSEQ', 'ENTREZID', 'GOSLIM_IDS')
+      df_mapping <- 
+        dplyr::left_join(agilent_table, id_map, dplyr::join_by(EnsemblID == ensembl_transcript_id)) %>%
+        dplyr::select(ProbeID, ensembl_gene_id)
 
-    unique_probe_ids$STRING_id <- NA_character_
+      names(df_mapping) <- c('ProbeName', 'ensembl_gene_id')
 
-    gene_col <- 'ENSEMBL'
-    if (sum(!is.na(unique_probe_ids$ENTREZID)) > sum(!is.na(unique_probe_ids$ENSEMBL))) {
-      gene_col <- 'ENTREZID'
-    }
-    if (sum(!is.na(unique_probe_ids$SYMBOL)) > max(sum(!is.na(unique_probe_ids$ENTREZID)), sum(!is.na(unique_probe_ids$ENSEMBL)))) {
-      gene_col <- 'SYMBOL'
+    } else if (annot_type == 'custom') {
+      unique_probe_ids <- read.csv(
+        file.path(local_annotation_dir, annot_config$annot_filename[[1]]),
+        header = TRUE, na.strings = c('NA', '')
+      )
+    } else {
+      annot_cols <- c('ProbeName', 'ENTREZID', 'SYMBOL', 'GENENAME', 'ENSEMBL', 'REFSEQ', 'GOSLIM_IDS', 'STRING_id', 'count_gene_mappings', 'gene_mapping_source')
+      unique_probe_ids <- setNames(data.frame(matrix(NA_character_, nrow = 1, ncol = length(annot_cols))), annot_cols)
     }
-
-    unique_probe_ids <- unique_probe_ids %>%
-                        dplyr::mutate( 
-                          count_gene_mappings = 1 + stringr::str_count(get(gene_col), stringr::fixed("|")),
-                          gene_mapping_source = gene_col
-                        )
-  } else if (annot_type == 'custom') {
-    unique_probe_ids <- read.csv(
-      file.path(local_annotation_dir, annot_config$annot_filename[[1]]),
-      header = TRUE, na.strings = c('NA', '')
-    )
-  } else {
-    annot_cols <- c('ProbeName', 'ENTREZID', 'SYMBOL', 'GENENAME', 'ENSEMBL', 'REFSEQ', 'GOSLIM_IDS', 'STRING_id', 'count_gene_mappings', 'gene_mapping_source')
-    unique_probe_ids <- setNames(data.frame(matrix(NA_character_, nrow = 1, ncol = length(annot_cols))), annot_cols)
   }
-} else {
+}
+
+# At this point, we have df_mapping from either the BioMart live service directly or via Agilent probe 
+# annotation or the ensembl genomes ftp archive depending on the organism and presence of BioMart probe 
+# annotation. If no df_mapping obtained (e.g., not supported in BioMart or lacking Ensembl transcript 
+# IDs), use custom annotations; otherwise, merge in-house annotations to df_mapping
+
+if (!use_custom_annot || annot_type == 'agilent') {
   annot <- read.table(
       as.character(annotation_file_path),
       sep = "\t",
@@ -1217,19 +1234,20 @@ if (use_custom_annot) {
   )
 
   unique_probe_ids <- df_mapping %>% 
-                        dplyr::mutate(dplyr::across(!!sym(expected_attribute_name), as.character)) %>% # Ensure probe ids treated as character type
-                        dplyr::group_by(!!sym(expected_attribute_name)) %>% 
-                        dplyr::summarise(
-                          ENSEMBL = list_to_unique_piped_string(ensembl_gene_id)
-                          ) %>%
-                        # Count number of ensembl IDS mapped
-                        dplyr::mutate( 
-                          count_gene_mappings = 1 + stringr::str_count(ENSEMBL, stringr::fixed("|")),
-                          gene_mapping_source = annot_key
-                        ) %>%
-                        dplyr::left_join(annot, by = c("ENSEMBL" = annot_key))
+    dplyr::mutate(dplyr::across(!!sym(expected_attribute_name), as.character)) %>% # Ensure probe ids treated as character type
+    dplyr::group_by(!!sym(expected_attribute_name)) %>% 
+    dplyr::summarise(
+      ENSEMBL = list_to_unique_piped_string(ensembl_gene_id)
+      ) %>%
+    # Count number of ensembl IDS mapped
+    dplyr::mutate( 
+      count_gene_mappings = 1 + stringr::str_count(ENSEMBL, stringr::fixed("|")),
+      gene_mapping_source = annot_key
+    ) %>%
+    dplyr::left_join(annot, by = c("ENSEMBL" = annot_key))
 }
 
+# reformat merged probe annotations
 norm_data$genes <- norm_data$genes %>% 
   dplyr::left_join(unique_probe_ids, by = c("ProbeName" = expected_attribute_name ) ) %>%
   dplyr::mutate( count_gene_mappings := ifelse(is.na(count_gene_mappings), 0, count_gene_mappings) ) %>%

Microarray/Agilent_1-channel/README.md

@@ -26,5 +26,6 @@
 Jonathan Oribello  
 
 **Maintained by:**  
-Crystal Han (crystal.han@nasa.gov)  
+Jihan Yehia (jihan.yehia@nasa.gov)
 Alexis Torres (alexis.torres@nasa.gov)
+Crystal Han (crystal.han@nasa.gov)  

Microarray/Agilent_1-channel/Workflow_Documentation/NF_MAAgilent1ch/README.md

@@ -4,7 +4,7 @@
 
 ### Implementation Tools <!-- omit in toc -->
 
-The current GeneLab Agilent 1 Channel Microarray consensus processing pipeline (NF_MAAgilent1ch), [GL-DPPD-7112-A](../../Pipeline_GL-DPPD-7112_Versions/GL-DPPD-7112-A.md), is implemented as a [Nextflow](https://nextflow.io/) DSL2 workflow and utilizes [Singularity](https://docs.sylabs.io/guides/3.10/user-guide/introduction.html) to run all tools in containers. This workflow (NF_MAAgilent1ch) is run using the command line interface (CLI) of any unix-based system.  While knowledge of creating workflows in Nextflow is not required to run the workflow as is, [the Nextflow documentation](https://nextflow.io/docs/latest/index.html) is a useful resource for users who want to modify and/or extend this workflow.   
+The current GeneLab Agilent 1 Channel Microarray consensus processing pipeline (NF_MAAgilent1ch), [GL-DPPD-7112-A](../../Pipeline_GL-DPPD-7112_Versions/GL-DPPD-7112-A.md), is implemented as a [Nextflow](https://nextflow.io/) DSL2 workflow and utilizes [Singularity](https://docs.sylabs.io/guides/3.10/user-guide/introduction.html) or [Docker](https://docs.docker.com/get-started/) containers to run all tools. This workflow (NF_MAAgilent1ch) is run using the command line interface (CLI) of any unix-based system.  While knowledge of creating workflows in Nextflow is not required to run the workflow as is, [the Nextflow documentation](https://docs.seqera.io/nextflow/) is a useful resource for users who want to modify and/or extend this workflow.   
 
 ### Workflow & Subworkflows <!-- omit in toc -->
 
@@ -33,13 +33,13 @@ Below is a description of each subworkflow and the additional output files gener
      - This subworkflow performs validation and verification (V&V) on the raw and processed data files.  It performs a series of checks on the output files generated and flags the results, using the flag codes indicated in the table below, which are outputted into a log file.  
        **V&V Flags**:
 
-       |Flag Codes|Flag Name|Interpretation|
-       |:---------|:--------|:-------------|
-       | 2    | MANUAL   | Special flag that indicates a manual check that is advised. Often used to advise what assess in QA plots. |
-       | 20    | GREEN   | Indicates the check passed all validation conditions |
-       | 30    | YELLOW  | Indicates the check was flagged for minor issues (e.g. slight outliers) |
-       | 50    | RED     | Indicates the check was flagged for moderate issues (e.g. major outliers) |
-       | 80    | HALT    | Indicates the check was flagged for severe issues that trigger a processing halt (e.g. missing data) |
+       | Flag Codes | Flag Name | Interpretation                                                                                            |
+       | :--------- | :-------- | :-------------------------------------------------------------------------------------------------------- |
+       | 2          | MANUAL    | Special flag that indicates a manual check that is advised. Often used to advise what assess in QA plots. |
+       | 20         | GREEN     | Indicates the check passed all validation conditions                                                      |
+       | 30         | YELLOW    | Indicates the check was flagged for minor issues (e.g. slight outliers)                                   |
+       | 50         | RED       | Indicates the check was flagged for moderate issues (e.g. major outliers)                                 |
+       | 80         | HALT      | Indicates the check was flagged for severe issues that trigger a processing halt (e.g. missing data)      |
 
 <br>
 
@@ -64,9 +64,11 @@ Below is a description of each subworkflow and the additional output files gener
 
 #### 1a. Install Nextflow
 
-Nextflow can be installed either through [Anaconda](https://anaconda.org/bioconda/nextflow) or as documented on the [Nextflow documentation page](https://www.nextflow.io/docs/latest/getstarted.html).
+Nextflow can be installed either through the [Anaconda bioconda channel](https://anaconda.org/bioconda/nextflow) or as documented in the [Nextflow installation documentation](https://docs.seqera.io/nextflow/install).
 
-> Note: If you want to install Anaconda, we recommend installing a Miniconda, Python3 version appropriate for your system, as instructed by [Happy Belly Bioinformatics](https://astrobiomike.github.io/unix/conda-intro#getting-and-installing-conda).  
+
+> [!TIP]
+> If you wish to install Conda, we recommend installing a Miniforge version appropriate for your system, as documented on the [conda-forge website](https://conda-forge.org/download/), where you can find basic binaries for most systems. More detailed miniforge documentation is available in the [miniforge github repository](https://github.com/conda-forge/miniforge).
 > 
 > Once conda is installed on your system, you can install the latest version of Nextflow by running the following commands:
 > 
@@ -83,7 +85,9 @@ Singularity is a container platform that allows usage of containerized software.
 
 We recommend installing Singularity on a system wide level as per the associated [documentation](https://docs.sylabs.io/guides/3.10/admin-guide/admin_quickstart.html).
 
-> Note: Singularity is also available through [Anaconda](https://anaconda.org/conda-forge/singularity).
+> [!TIP]
+> - Singularity is also available through the [Anaconda conda-forge channel](https://anaconda.org/conda-forge/singularity). 
+> - Alternatively, Docker can be used in place of Singularity. To get started with Docker, see the [Docker CE installation documentation](https://docs.docker.com/engine/install/).
 
 <br>
 
@@ -106,7 +110,8 @@ unzip NF_MAAgilent1ch_1.0.5.zip
 ### 3. Run the Workflow
 
 While in the location containing the `NF_MAAgilent1ch_1.0.5` directory that was downloaded in [step 2](#2-download-the-workflow-files), you are now able to run the workflow. Below are three examples of how to run the NF_MAAgilent1ch workflow:
-> Note: Nextflow commands use both single hyphen arguments (e.g. -help) that denote general nextflow arguments and double hyphen arguments (e.g. --ensemblVersion) that denote workflow specific parameters.  Take care to use the proper number of hyphens for each argument.
+> [!NOTE]
+> Nextflow commands use both single hyphen arguments (e.g. -help) that denote general nextflow arguments and double hyphen arguments (e.g. --ensemblVersion) that denote workflow specific parameters.  Take care to use the proper number of hyphens for each argument.
 
 <br>
 
@@ -196,7 +201,7 @@ All parameters listed above and additional optional arguments for the NF_MAAgile
 nextflow run NF_MAAgilent1ch_1.0.5/main.nf --help
 ```
 
-See `nextflow run -h` and [Nextflow's CLI run command documentation](https://nextflow.io/docs/latest/cli.html#run) for more options and details common to all nextflow workflows.
+See `nextflow run -h` and [Nextflow's CLI run command documentation](https://docs.seqera.io/nextflow/cli) for more options and details common to all nextflow workflows.
 
 <br>
 
@@ -211,7 +216,7 @@ All R code steps and output are rendered within a Quarto document yielding the f
   
 
 The outputs from the Analysis Staging and V&V Pipeline Subworkflows are described below:
-> Note: The outputs from the Agilent 1 Channel Microarray Processing Subworkflow are documented in the [GL-DPPD-7112-A.md](../../../Pipeline_GL-DPPD-7112_Versions/GL-DPPD-7112-A.md) processing protocol.
+> Note: The outputs from the Agilent 1 Channel Microarray Processing Subworkflow are documented in the [GL-DPPD-7112-A.md](../../Pipeline_GL-DPPD-7112_Versions/GL-DPPD-7112-A.md) processing protocol.
 
 **Analysis Staging Subworkflow**
 > Note: only applicable for [Approach 1](#3a-approach-1-run-the-workflow-on-a-genelab-agilent-1-channel-microarray-dataset) and [Approach 3](#3c-approach-3-run-the-workflow-using-an-isa-archive)
@@ -229,7 +234,7 @@ The outputs from the Analysis Staging and V&V Pipeline Subworkflows are describe
 <br>
 
 Standard Nextflow resource usage logs are also produced as follows:
-> Further details about these logs can also found within [this Nextflow documentation page](https://www.nextflow.io/docs/latest/tracing.html#execution-report).
+> Further details about these logs can also found in the [Nextflow Report Documentation](https://docs.seqera.io/nextflow/reports).
 
 **Nextflow Resource Usage Logs**
    - Output:

Microarray/Agilent_1-channel/Workflow_Documentation/NF_MAAgilent1ch/examples/annotations/README.md

@@ -16,8 +16,8 @@
 | Column Name | Type | Description | Example |
 |:------------|:-----|:------------|:--------|
 | array_design | string | A bioMart attribute identifier denoting the microarray probe/probeset attribute used for annotation mapping. | AGILENT SurePrint G3 GE 8x60k v3 |
-| annot_type | string | Used to determine how the custom annotations are parsed before merging to the data. Currently, only the below are supported: <ul><li>`agilent`: Annotations file is expected to be in the AA (All Annotations) format by [Agilent](https://earray.chem.agilent.com/earray/)</li><li>`custom`: Annotations file is merged as is, expected to have the following columns: `ProbesetID`, `ENTREZID`, `SYMBOL`, `GENENAME`, `ENSEMBL`, `REFSEQ`, `GOSLIM_IDS`, `STRING_id`, `count_gene_mappings`, `gene_mapping_source`</li></ul> | agilent |
-| annot_filename | string | Name of the custom annotations file. This is the AllAnnotations file downloaded from Agilent's eArray web portal. | 072363_D_AA_20240521.txt |
+| annot_type | string | Used to determine how the custom annotations are parsed before merging to the data. Currently, only the below are supported: <ul><li>`agilent`: Annotations file is expected to be in the AA (All Annotations) format created by [Agilent](https://earray.chem.agilent.com/earray/) for catalog arrays. The Ensembl Transcript Stable IDs in the file are used to map to Ensembl Gene Stable IDs after which the standard GeneLab annotation process is used.</li><li>`custom`: Annotations file is merged as-is and expected to have the following columns: `ProbeName`, `ENTREZID`, `SYMBOL`, `GENENAME`, `ENSEMBL`, `REFSEQ`, `GOSLIM_IDS`, `STRING_id`, `count_gene_mappings`, `gene_mapping_source`, where `ENSEMBL` is the Ensembl Gene Stable ID, not the Transcript ID</li></ul> | agilent |
+| annot_filename | string | Name of the custom annotations file. | 072363_D_AA_20240521.txt |
 
 ## Optional columns 
 If the file was downloaded from a website, provide the download link used and date