Count the pseudogenes AMAP called in each sample in the report's chart

CHANGELOG.md

@@ -35,6 +35,7 @@ directly comparable with 1.0.1.
 - **Genes without a COG category are counted once**: `count_annot.py` filled in a missing category as `Unassigned` and split it into letters, as it splits `JK` into J and K, so every gene eggNOG gave no category was counted once in A, D, E, G, I and U and twice in N and S, in `*-AMAP.COG.tsv` and in the report's COG figures and statistics, where it also shrank every other category's share. A test-profile run with a stand-in for eggNOG that gave none of its 4036 genes a category reported 4036 genes in each of those six and 8072 in N and S. Such a gene now counts once, as `NaN`, which the report shows as None, and its percentage is of all the genome's genes like every other category's
 - **COG categories R, W, X, Y and Z in the report**: its COG figures took their categories from a list that ran from A to V without R, so a gene eggNOG put in R, W, X, Y or Z was plotted as NA in the overview and in none of the per-category figures, which stayed empty for those five, X - Mobilome (prophages and transposons) among them. The categories now come from the report's own list of figures, A to Z and None, so each is plotted in its own. The statistics table never used that list and is unchanged
 - **The QC table counts the pseudogenes AMAP called**: `annotation_qc.py` read the first of `Is_Pseudogene`, `pseudogene` and `is_pseudogene` that the master table has, and every master table has Mycobrowser's `Is_Pseudogene`, empty when Mycobrowser is off. So `pseudogene_count` in `*-AMAP.qc.tsv` and the report's QC table was 0 in every organism-agnostic run (E. coli Sakai against K-12, where AMAP called 129 pseudogenes), and in an MTBC run the number of the sample's genes named after one of the 13 that Mycobrowser marks as pseudogenes in H37Rv, whatever AMAP made of them: 13 for CDC1551, where AMAP calls 123 and finds 6 of those 13 intact. It now counts AMAP's own call, the master table's `pseudogene` column, the same genes the GFF3 marks `pseudo=true` and the GenBank `/pseudo`
+- **The report's pseudogene chart counts AMAP's calls, sample by sample**: its Pseudogene Detection section counted the master table's genes flagged in Mycobrowser's `Is_Pseudogene` or in AMAP's `pseudogene` column, lineage by lineage, as "Pseudogenes per Lineage (Mycobrowser annotation)". `Is_Pseudogene` says which genes are pseudogenes in H37Rv, not in the sample, so CDC1551 came out at 129 where AMAP calls 123 and the QC table says 123, the other 6 being genes Mycobrowser marks that CDC1551 has intact (Rv0947c, Rv1150, Rv2023A, Rv2098c, Rv3021c, Rv3324A), and without Mycobrowser the chart showed AMAP's calls under Mycobrowser's name. It also left out a lineage with none, and added up a lineage's samples, so a lineage with more samples looked as if it had more pseudogenes. REPORT_INPUTS now hands the report every master table row with the sample it is from, and the chart shows how many pseudogenes AMAP called in each sample, by lineage, a sample with none at 0, next to a table of them. The master table explorer shows the `note` that says why a gene is a pseudogene
 - **A sample with no gene in a category counts 0 in the report**: `count_annot.py` writes a row for each COG or Mycobrowser category a sample has genes in, and the report plotted the rows as they came, so a sample with no gene in a category was left out of that category's figures, of the Dunn tests between lineages run on them, of the overview and of the statistics table, instead of counting 0. Its lineage's box, mean and median were those of the samples that have some, a lineage none of whose samples has any dropped out of the comparison, and the comparison with the median never showed a sample that had lost a category. Each sample now counts 0 in every category another sample of the run has genes in, and the statistics table's Count is every sample of the lineage. A category no sample has keeps its empty figure, and `*-AMAP.COG.tsv` and `*-AMAP.MycoBrowser.tsv` are unchanged
 - **The clean proteome leaves out the pseudogenes AMAP called**: `export_clean_proteome.py` took them from the first of `Is_Pseudogene`, `pseudogene` and `is_pseudogene` that the master table has, which is Mycobrowser's, and looked each protein up by its whole FASTA ID, `L4_CDC1551_Rv0001;+;0-1524` where the master table has `L4_CDC1551_Rv0001`, so it never found one and its log said 0 pseudogenes. Only the proteins with an internal stop codon were left out, and the pseudogenes without one stayed in `*-AMAP.clean.faa`, which the docs say has none: 42 of the 123 AMAP calls in CDC1551 (a frameshift or lost stop codon, no stop codon within reach, no stop codon as in the reference) and 27 of the 129 in E. coli Sakai against K-12. It now leaves out the genes whose `pseudogene` is yes, found by the ID before the `;`, and on both genomes holds the CDS the GFF3 does not mark `pseudo=true`, 3891 and 5194, with the proteins the GenBank translates for them. A master table without that column stops it, where a warning let every pseudogene through
 
@@ -123,7 +124,8 @@ directly comparable with 1.0.1.
 - **`tests/test_environment_pins.py`** fails when `environment.yml` pins a tool at a version `docker/Dockerfile` does not install
 - **`tests/test_citation.py`** fails when anything cites AMAP with a DOI other than the concept DOI of the manifest
 - **`tests/test_clean_proteome.py`** runs the steps after the correction as `tests/test_sample_ids.py` chains them, with FIX_CDS calling some genes pseudogenes, one with a stop codon in its protein and the others without, and Mycobrowser marking a gene AMAP finds intact, with Mycobrowser on and off, and checks that `*-AMAP.clean.faa` holds every protein FIX_CDS wrote but the pseudogenes', as it wrote them: the CDS the GFF3 does not mark `pseudo=true`, with the GenBank's translations
-- **732 tests**, up from 92
+- **`tests/test_report_pseudogenes.py`** runs the steps after the correction as `tests/test_sample_ids.py` chains them, with FIX_CDS calling a different number of pseudogenes in each sample, none in one, and Mycobrowser marking a gene that every sample has intact, with Mycobrowser on and off, and checks that REPORT_INPUTS hands over each master table row with its sample; where R is installed, it runs the report's own pseudogene section on what REPORT_INPUTS made and checks that it draws each sample at the number AMAP called, the QC table's, a sample with none at 0, under a title that does not name Mycobrowser
+- **744 tests**, up from 92
 
 ## [1.0.1] - 2026-07-25
 

bin/report/demo/make_demo.py

@@ -64,7 +64,8 @@ def w(name, lines):
     "pseudogene_count", "hypothetical_pct", "mean_gene_length",
     "median_gene_length", "coding_density_pct", "eggnog_annotation_pct"]))
 
-master_header = ["simple_id", "lineage", "start", "end", "strand", "rvs",
+# as REPORT_INPUTS writes it: the sample each gene is from, then its master table
+master_header = ["sample", "simple_id", "lineage", "start", "end", "strand", "rvs",
                  "pseudogene", "COG_category", "Description", "Preferred_name",
                  "KEGG_Pathway", "Functional_Category"]
 master_rows.append("\t".join(master_header))
@@ -126,7 +127,7 @@ def w(name, lines):
         cat = rng.choices([m[0] for m in MYCO], [m[1] for m in MYCO])[0]
         s = rng.randint(1, 4_400_000)
         master_rows.append("\t".join([
-            f"{sample}_{i:05d}", lin, str(s), str(s + rng.randrange(150, 3000, 3)),
+            sample, f"{sample}_{i:05d}", lin, str(s), str(s + rng.randrange(150, 3000, 3)),
             rng.choice("+-"), f"Rv{rng.randint(1, 3924):04d}",
             rng.choice(["no"] * 24 + ["yes"]), rng.choice(COGS),
             "hypothetical protein" if cat == "conserved hypotheticals" else "putative enzyme",

bin/report/report_AMAP.Rmd

@@ -403,7 +403,8 @@ if (has_rvmissing) {
   } else { has_rvmissing <- FALSE }
 }
 
-# Master table
+# Master table: every sample's genes, each row with the sample it is from in
+# the column REPORT_INPUTS adds, the name the sample's files carry
 has_mastertable <- file.exists(params$mastertable_file) && file.info(params$mastertable_file)$size > 0
 if (has_mastertable) {
   mastertable <- tryCatch(
@@ -909,44 +910,52 @@ if (has_mastertable && exists("kegg_data")) {
 
 ## Pseudogene Detection {.tabset}
 
-Genes without a stop codon in their annotated coordinates — potential pseudogenes.
+The genes AMAP calls pseudogenes in each sample, by lineage: internal stop codons, a frameshift or lost stop that makes the frame run on past the end of the reference protein, no stop codon within reach, or a stop the reference itself lacks. They are the genes the GFF3 marks `pseudo=true` and the QC table counts, and the master table's `note` column says why each one is a pseudogene. A sample with none counts 0. Mycobrowser's `Is_Pseudogene` says whether a gene is a pseudogene in H37Rv, not in the sample, and is not counted.
 
-### Per Lineage
+### By Lineage
 ```{r pseudogenes, echo=FALSE, fig.align='center', fig.width=12, fig.height=6, message=FALSE, warning=FALSE}
-if (has_mastertable) {
-  # Pseudogenes can be detected from genes with internal stops or no terminal stop
-  # Using the rvs column: genes with known Rv ID that were flagged as pseudogenes
-  # Check both Mycobrowser's Is_Pseudogene and AMAP's pseudogene columns
-  has_myco_pseudo <- "Is_Pseudogene" %in% colnames(mastertable)
-  has_amap_pseudo <- "pseudogene" %in% colnames(mastertable)
-  if (has_myco_pseudo || has_amap_pseudo) {
-    pseudo_data <- mastertable %>%
-
-      filter(
-        (if (has_myco_pseudo) (Is_Pseudogene %in% c("yes","Yes","TRUE",TRUE)) else FALSE) |
-        (if (has_amap_pseudo) (pseudogene %in% c("yes","Yes","TRUE",TRUE)) else FALSE)
-      )
-
-    if (nrow(pseudo_data) > 0) {
-      pseudo_count <- pseudo_data %>%
-        count(lineage) %>%
-        arrange(desc(n))
-
-      p <- ggplot(pseudo_count, aes(x=reorder(lineage, -n), y=n, fill=lineage)) +
-        geom_col(show.legend=FALSE) +
-        geom_text(aes(label=n), vjust=-0.3, size=3.5) +
-        scale_fill_manual(values=colors) +
-        labs(title="Pseudogenes per Lineage (Mycobrowser annotation)", x="Lineage", y="Count") +
-        theme_pubr()
-      print(p)
-    } else {
-      cat("No pseudogenes detected in the dataset.\n")
-    }
+if (!has_mastertable) {
+  cat("Master table not available.\n")
+} else if (!all(c("sample", "pseudogene") %in% colnames(mastertable))) {
+  cat("The master table has no sample or no pseudogene column.\n")
+} else {
+  # AMAP's own call, as the QC table counts it. Mycobrowser's Is_Pseudogene
+  # was counted too, which says which genes are pseudogenes in H37Rv: CDC1551
+  # came out at 129, where AMAP calls 123, the other 6 being genes Mycobrowser
+  # marks that are intact in CDC1551; and without Mycobrowser the chart showed
+  # AMAP's calls under its name. Counted sample by sample, a sample with none
+  # at 0: counting the flagged rows dropped a lineage with none, and adding up
+  # a lineage's samples made one with more samples look as if it had more.
+  pseudo_count <- mastertable %>%
+    mutate(lineage = ifelse(is.na(lineage) | lineage == "", "NA", as.character(lineage))) %>%
+    group_by(Lineage = lineage, Sample = sample) %>%
+    summarise(Pseudogenes = sum(pseudogene %in% "yes"), .groups = 'drop') %>%
+    arrange(Lineage, Sample)
+
+  if (sum(pseudo_count$Pseudogenes) == 0) {
+    cat("AMAP called no pseudogene in any sample.\n")
   } else {
-    cat("No pseudogene columns found in master table (neither Is_Pseudogene nor pseudogene).\n")
+    p <- ggplot(pseudo_count, aes(x=Lineage, y=Pseudogenes, color=Lineage)) +
+      geom_boxplot(outlier.shape=NA, fill=NA, alpha=0.6) +
+      # side to side only, so each sample is drawn at its count
+      geom_jitter(alpha=0.7, width=0.25, height=0) +
+      scale_color_manual(values=colors) +
+      scale_y_continuous(breaks=function(x) unique(floor(pretty(x)))) +
+      expand_limits(y=0) +
+      labs(title="Pseudogenes AMAP Called per Sample", x="Lineage", y="Pseudogenes per sample") +
+      theme_pubr() +
+      theme(legend.position="none")
+    print(p)
   }
+}
+```
+
+### Per Sample
+```{r pseudogenes_table, echo=FALSE, message=FALSE, warning=FALSE}
+if (exists("pseudo_count")) {
+  make_dt(pseudo_count, "pseudogenes_per_sample")
 } else {
-  cat("Master table not available.\n")
+  cat("No pseudogene counts available.\n")
 }
 ```
 
@@ -1214,9 +1223,9 @@ Interactive explorer for the full annotation master table with all functional an
 
 ```{r master_explorer, echo=FALSE, message=FALSE, warning=FALSE}
 if (has_mastertable) {
-  # Show a subset of the most useful columns
+  # Show a subset of the most useful columns, a pseudogene's reason in its note
   display_cols <- intersect(
-    c("id","strand","start","end","rvs","pseudogene","COG_category","Description","Preferred_name",
+    c("id","strand","start","end","rvs","pseudogene","note","COG_category","Description","Preferred_name",
       "Gene","Product_x","KEGG_Pathway","KEGG_ko","GOs","EC","PFAMs","CAZy",
       "Functional_Category","Is_Pseudogene","Locus"),
     colnames(mastertable))