From da77d2d6ea0114ea462d04254c78a6242999c9b2 Mon Sep 17 00:00:00 2001 From: Paula Ruiz Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Sun, 18 Oct 2026 21:11:33 +0200 Subject: [PATCH] Count the pseudogenes AMAP called in each sample in the report's chart The report's Pseudogene Detection chart counted the genes flagged by Mycobrowser's Is_Pseudogene or AMAP's pseudogene column, summed by lineage under a Mycobrowser title, so CDC1551 showed 129 where AMAP calls 123, a lineage with none dropped out and a lineage with more samples looked as if it had more. REPORT_INPUTS now gives each master table row the sample it is from, and the chart draws every sample at the number AMAP called, the QC table's, a sample with none at 0, with a table of the counts beside it. The master table explorer shows the note that says why a gene is a pseudogene, and tests/test_report_pseudogenes.py runs the chained steps and the section's own R code. --- CHANGELOG.md | 4 +- bin/report/demo/make_demo.py | 5 +- bin/report/report_AMAP.Rmd | 81 +++++++------- docs/about/changelog.html | 4 +- docs/assets/search-index.js | 2 +- docs/user-guide/report.html | 2 +- docs/user-guide/report.md | 4 +- modules/local/freport.nf | 10 +- tests/test_report_pseudogenes.py | 178 +++++++++++++++++++++++++++++++ 9 files changed, 242 insertions(+), 48 deletions(-) create mode 100644 tests/test_report_pseudogenes.py diff --git a/CHANGELOG.md b/CHANGELOG.md index cb7b4e7..08e8735 100644 --- a/CHANGELOG.md +++ b/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 ### ๐Ÿ› Pipeline @@ -107,10 +108,11 @@ directly comparable with 1.0.1. - **`tests/test_duplication_report.py`** runs GET_DUPLICATIONS and GET_ALLDUP as Nextflow renders them, chained as the pipeline chains them, with eggNOG giving one copy of a duplicated gene a category, another none and the last no row at all, and checks the columns and COG counts of `homologous_communities.tsv`; where R is installed, it knits the report's duplication summary on that table and on the one `GET_ALLDUP` writes when nothing is duplicated. It also fails when a chunk of the report prints markdown without passing it on as markdown - **`tests/test_profile_caps.py`** fails when a profile sets a cap, or when the order of `-profile` changes one: it reads `nextflow.config` and applies its profiles in every order, as both of Nextflow's parsers do, and with Nextflow on PATH also asks `nextflow config` - **`tests/test_qc_pseudogenes.py`** runs the steps after the correction as `tests/test_sample_ids.py` chains them, with FIX_CDS calling some genes pseudogenes and Mycobrowser marking another, with Mycobrowser on and off, and checks that `*-AMAP.qc.tsv`, REPORT_INPUTS and the report's own QC table count the pseudogenes AMAP called, as many as the GFF3 and the GenBank mark +- **`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 - **`tests/test_duplication_reproducible.py`** runs `get_duplications.py` and `network_genomes.py` under three hash seeds, on the Diamond tables of three genomes with a family whose average clustering changes with the order of its proteins, and fails unless every file they write is the same under all three; it also fails when either script opens a web browser, with hits or without, as a stand-in browser that records what it is asked to open tells - **`tests/test_custom_reference.py`** extracts a reference whose plasmid has a gene across its origin from a GenBank record and from GFF3 written either way, on either strand, compares Bakta's calls with it sequence by sequence, and runs `PREPARE_REFERENCE`'s own script, rendered as Nextflow renders it, with a stand-in for Bakta that names sequences as Bakta does - **`tests/test_amap_report.py`** writes the report from four genomes of a real run cut to two stretches of the reference (`tests/data/report/`, rebuilt by its `make_fixture.py`), where M. bovis lacks RD4 and a contig end cuts dnaE2, plus a gene across the origin of a circular plasmid, and checks that every gene, and every Bakta call AMAP did not keep, comes back out of the page as it went in, with Python and, where node is installed, with the page's own JavaScript, along with the states, the flags, the correction's steps and what the process hands the script -- **682 tests**, up from 92 +- **694 tests**, up from 92 ## [1.0.1] - 2026-07-25 diff --git a/bin/report/demo/make_demo.py b/bin/report/demo/make_demo.py index 9d87ef1..c27f1d0 100644 --- a/bin/report/demo/make_demo.py +++ b/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", diff --git a/bin/report/report_AMAP.Rmd b/bin/report/report_AMAP.Rmd index 55c0876..f7865aa 100644 --- a/bin/report/report_AMAP.Rmd +++ b/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)) diff --git a/docs/about/changelog.html b/docs/about/changelog.html index 9992ee1..05efef3 100644 --- a/docs/about/changelog.html +++ b/docs/about/changelog.html @@ -94,7 +94,7 @@

Changelog

[Unreleased]#

Annotating H37Rv against its own reference now gives back every CDS of the AL123456.3 GenBank record at exactly its coordinates, strand and name: 4031 of 4031, against 4020 before. The changes apply to every genome, so results are not directly comparable with 1.0.1.

๐Ÿ”ฌ Correctness#

- +

๐Ÿ› Pipeline#

๐Ÿ“Š Report#

@@ -104,7 +104,7 @@

๐Ÿงน Code๐Ÿ“š Documentation#

๐Ÿงช Testing#

- +

[1.0.1] - 2026-07-25#

Bug fixes, most of them found by running the pipeline on real genomes rather than by reading it. Several changed what the pipeline outputs, so results from 1.0.0 are not directly comparable.

๐Ÿ”ฌ Correctness#

diff --git a/docs/assets/search-index.js b/docs/assets/search-index.js index dc78ce6..118ef02 100644 --- a/docs/assets/search-index.js +++ b/docs/assets/search-index.js @@ -454,7 +454,7 @@ window.AMAP_DOCS_INDEX = [ "page": "The report", "section": "Functional annotation", "url": "user-guide/report.html#functional-annotation", -"text": "COG categories from eggNOG, Mycobrowser functional categories for MTBC, and a comparison against the reference distribution. Each has its own per-category breakdown in counts and percentages, with the statistical tests alongside. PE/PPE gene families get their own section, since they are both the hardest part of an MTBC genome to annotate and usually the most interesting. The pseudogenes of each lineage close it." +"text": "COG categories from eggNOG, Mycobrowser functional categories for MTBC, and a comparison against the reference distribution. Each has its own per-category breakdown in counts and percentages, with the statistical tests alongside. PE/PPE gene families get their own section, since they are both the hardest part of an MTBC genome to annotate and usually the most interesting. The pseudogenes close it: how many AMAP called in each sample, lineage by lineage, a sample with none at 0, the same genes the QC table counts and the GFF3 marks pseudo=true. The master table's note column says why each is a " }, { "page": "The report", diff --git a/docs/user-guide/report.html b/docs/user-guide/report.html index 31674eb..3e227b0 100644 --- a/docs/user-guide/report.html +++ b/docs/user-guide/report.html @@ -103,7 +103,7 @@

Annotation resultsFunctional annotation#

COG categories from eggNOG, Mycobrowser functional categories for MTBC, and a comparison against the reference distribution. Each has its own per-category breakdown in counts and percentages, with the statistical tests alongside.

-

PE/PPE gene families get their own section, since they are both the hardest part of an MTBC genome to annotate and usually the most interesting. The pseudogenes of each lineage close it.

+

PE/PPE gene families get their own section, since they are both the hardest part of an MTBC genome to annotate and usually the most interesting. The pseudogenes close it: how many AMAP called in each sample, lineage by lineage, a sample with none at 0, the same genes the QC table counts and the GFF3 marks pseudo=true. The master table's note column says why each is a pseudogene.

KEGG pathways#

The most represented metabolic pathways, and how they differ between lineages.

Duplications#

diff --git a/docs/user-guide/report.md b/docs/user-guide/report.md index c6e1f40..0b43828 100644 --- a/docs/user-guide/report.md +++ b/docs/user-guide/report.md @@ -47,7 +47,9 @@ breakdown in counts and percentages, with the statistical tests alongside. PE/PPE gene families get their own section, since they are both the hardest part of an MTBC genome to annotate and usually the most interesting. The -pseudogenes of each lineage close it. +pseudogenes close it: how many AMAP called in each sample, lineage by lineage, +a sample with none at 0, the same genes the QC table counts and the GFF3 marks +`pseudo=true`. The master table's `note` column says why each is a pseudogene. ### KEGG pathways diff --git a/modules/local/freport.nf b/modules/local/freport.nf index 06cd594..c3f033d 100644 --- a/modules/local/freport.nf +++ b/modules/local/freport.nf @@ -42,17 +42,19 @@ process REPORT_INPUTS { tail -n +2 "\$f" | awk -v s="\$sample" -F'\\t' '{print \$0 "\\t" s}' done >> input_rvmissing.txt - # Master tables โ€” use first NON-EMPTY file for header, strip rest + # Master tables: the first NON-EMPTY file's header, then every file's rows, + # each with the sample it is from (the name the file carries) in a first + # column, as the QC tables below, so FINAL_REPORT can count sample by sample touch input_mastertables.txt header_written=false find . -maxdepth 1 -name '*AMAPmastertable.tsv' -print0 | sort -z | while IFS= read -r -d '' f; do [ ! -s "\$f" ] && continue + sample=\$(basename "\$f" -AMAPmastertable.tsv) if ! \$header_written; then - cat "\$f" + head -1 "\$f" | awk '{print "sample\\t" \$0}' header_written=true - else - tail -n +2 "\$f" fi + tail -n +2 "\$f" | awk -v s="\$sample" '{print s "\\t" \$0}' done >> input_mastertables.txt # Pangene Rtab โ€” merge with outer join (handles different gene families) diff --git a/tests/test_report_pseudogenes.py b/tests/test_report_pseudogenes.py new file mode 100644 index 0000000..49b16cc --- /dev/null +++ b/tests/test_report_pseudogenes.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +The report's pseudogene chart counts the pseudogenes AMAP called, sample by +sample. + +Its Pseudogene Detection section kept the master table's genes flagged in +Mycobrowser's Is_Pseudogene or in AMAP's `pseudogene` column and counted them +by lineage, as "Pseudogenes per Lineage (Mycobrowser annotation)". Is_Pseudogene +says which genes are pseudogenes in H37Rv, not in the sample: CDC1551 came out +at 129, where AMAP calls 123 and its QC table says 123, the other 6 being genes +Mycobrowser marks that are intact in CDC1551. Without Mycobrowser the chart +showed AMAP's calls under Mycobrowser's name. Counting the flagged rows left +out a lineage with none, and adding up a lineage's samples made one with more +samples look as if it had more pseudogenes. REPORT_INPUTS now hands the report +the master tables with the sample each row is from, and the chart counts +AMAP's own call in each sample, a sample with none at 0. + +The tests run 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, as tests/test_qc_pseudogenes.py marks it, with Mycobrowser on and off; +then the report's own R code for the section on what REPORT_INPUTS made, where +R is installed. +""" + +import sys +from pathlib import Path + +import pytest + +ROOT = Path(__file__).resolve().parent.parent +sys.path.insert(0, str(ROOT / 'bin')) + +from tests import test_qc_pseudogenes as qc # noqa: E402 +from tests import test_sample_ids as chain # noqa: E402 +from tests.test_sample_ids import rscript, tools # noqa: E402,F401 (fixtures) + +# (the name a genome's files carry, its sample ID, its lineage). The k-th has +# Rv0001, Rv0002, Rv0003 and 2 + k copies of a gene of no reference. +SAMPLES = [('L2_Beijing_1', 'Beijing_1', 'L2'), + ('L2_Beijing_2', 'Beijing_2', 'L2'), + ('L4_H37Rv_test', 'H37Rv_test', 'L4'), + ('NA_K12_self', 'K12_self', 'NA')] # typing found no lineage + +# How many genes FIX_CDS calls pseudogenes in each. L2's two samples have 4 +# together, more than the 3 of lineage NA's one sample, which has more than +# either of them; L4's one sample has none. +PSEUDOGENES = {'L2_Beijing_1': 2, 'L2_Beijing_2': 2, 'L4_H37Rv_test': 0, 'NA_K12_self': 3} +REASONS = {'Rv0002': 'frameshift or lost stop codon: its frame runs on 40 residues ' + 'past the end of the reference protein', + 'Rv0003': 'no stop codon in its frame within 300 bp', + 'novel': '1 internal stop codon(s)'} + +# tests/test_sample_ids.py's writer, kept before the runs patch it +_write_sample = chain.write_sample + + +def write_sample(d, name, copies): + """tests/test_sample_ids.py's sample, with FIX_CDS calling as many of its + genes pseudogenes as PSEUDOGENES says, never Rv0001, which Mycobrowser + marks, nor the first copy of the gene of no reference, with the reason in + the note.""" + _write_sample(d, name, copies) + table = d / f'{name}.AMAP.tsv' + header, *rows = table.read_text().splitlines() + flagged = novel = 0 + for i, row in enumerate(rows): + cells = row.split('\t') + key = cells[4] or 'novel' + novel += key == 'novel' + if flagged < PSEUDOGENES[name] and key != 'Rv0001' and novel != 1: + cells[5], cells[8] = 'yes', REASONS[key] + rows[i] = '\t'.join(cells) + flagged += 1 + assert flagged == PSEUDOGENES[name] + table.write_text('\n'.join([header] + rows) + '\n') + return d + + +RUNS = ['mycobrowser', 'no_mycobrowser'] + + +@pytest.fixture(scope='module', params=RUNS) +def run(request, tmp_path_factory, tools): + """(the tasks of a run, whether Mycobrowser was on) for both runs.""" + mycobrowser = request.param == 'mycobrowser' + with pytest.MonkeyPatch.context() as mp: + mp.setattr(chain, 'write_sample', write_sample) + mp.setattr(chain, 'write_mycobrowser', qc.write_mycobrowser) + if not mycobrowser: + mp.setattr(chain, 'run_task', qc.skipping_mycobrowser(chain.run_task)) + tasks = chain.run_pipeline(tmp_path_factory.mktemp(request.param), SAMPLES, tools) + return tasks, mycobrowser + + +def section(tasks, rscript, tables, replace=None): + """Run the report's pseudogene section on a run's report inputs, as + tests/test_sample_ids.py runs the report's chunks, and read back `tables` + ({R expression: name}).""" + with pytest.MonkeyPatch.context() as mp: + mp.setattr(chain, 'R_CHUNKS', ['setup', 'pseudogenes']) + mp.setattr(chain, 'R_TABLES', tables) + return chain.report_tables(tasks, rscript, replace=replace) + + +@pytest.fixture(scope='module') +def chart(run, rscript): + """The count the section draws for each sample, and the chart's labels.""" + tasks, _ = run + return section(tasks, rscript, { + 'pseudo_count': 'counts', + 'data.frame(title = p$labels$title, x = p$labels$x, y = p$labels$y)': 'labels'}) + + +def per_sample(counts): + """{(lineage, sample): pseudogenes} from the counts the section draws.""" + assert {'Lineage', 'Sample', 'Pseudogenes'} <= set(counts.columns), \ + f'the chart counts {list(counts.columns)}, not each sample' + return dict(zip(zip(counts['Lineage'], counts['Sample']), counts['Pseudogenes'])) + + +def amap_calls(): + """{(lineage, sample): pseudogenes} as FIX_CDS called them.""" + return {(lineage, name): PSEUDOGENES[name] for name, _, lineage in SAMPLES} + + +class TestReportInputs: + + def test_each_master_table_row_says_its_sample(self, run): + # the name its file carries, as the QC table's rows say it + tasks, _ = run + header, *rows = (tasks['report'] / 'input_mastertables.txt').read_text().splitlines() + for name, _, _ in SAMPLES: + own = (tasks['outputs'][name] / f'{name}-AMAPmastertable.tsv').read_text().splitlines() + assert header == 'sample\t' + own[0] + assert [r.split('\t', 1)[1] for r in rows if r.split('\t', 1)[0] == name] == own[1:] + assert len(rows) == sum( + len((tasks['outputs'][name] / f'{name}-AMAPmastertable.tsv').read_text().splitlines()) - 1 + for name, _, _ in SAMPLES) + + +class TestChart: + + def test_counts_the_pseudogenes_amap_called_in_each_sample(self, chart): + # Main counted by lineage, Mycobrowser's Rv0001 in every sample with + # it on: L2 6, NA 4, L4 1; and without it L2 4, NA 3, and no L4 + assert per_sample(chart['counts']) == amap_calls() + + def test_a_sample_with_none_counts_zero(self, chart): + counts = per_sample(chart['counts']) + assert counts[('L4', 'L4_H37Rv_test')] == 0 + + def test_as_many_as_the_qc_table_says(self, run, chart): + tasks, _ = run + qc_table = chain.read_tsv(tasks['report'] / 'input_qc.txt') + assert {sample: n for (_, sample), n in per_sample(chart['counts']).items()} == \ + dict(zip(qc_table['sample'], qc_table['pseudogene_count'])) + + def test_the_title_names_amap_s_calls_per_sample(self, chart): + labels = chart['labels'].iloc[0] + assert 'Mycobrowser' not in labels['title'] + assert 'AMAP' in labels['title'] and 'per Sample' in labels['title'] + assert (labels['x'], labels['y']) == ('Lineage', 'Pseudogenes per sample') + + def test_every_sample_at_zero_when_none_has_any(self, run, rscript): + # no chart then, but the table still lists every sample, at 0 + tasks, _ = run + header, *rows = (tasks['report'] / 'input_mastertables.txt').read_text().splitlines() + col = header.split('\t').index('pseudogene') + intact = [] + for row in rows: + cells = row.split('\t') + cells[col] = 'no' + intact.append('\t'.join(cells)) + tables = section(tasks, rscript, {'pseudo_count': 'counts'}, + replace={'input_mastertables.txt': '\n'.join([header] + intact) + '\n'}) + assert per_sample(tables['counts']) == {key: 0 for key in amap_calls()}