Add per-population copy number (PPCN) normalization to anvi-script-gen-function-matrix-across-genomes

Adds a --per-population-copy-number flag that reports an additional matrix in
which each function's frequency in a metagenomic assembly is divided by the number
of populations estimated to be present in that assembly. Population counts are
estimated from single-copy core genes via anvi'o's NumGenomesEstimator (the mode
of SCG hits per domain, summed across Bacteria, Archaea, and Eukarya), the same
estimate anvi-display-contigs-stats reports. This normalization, introduced by
Iva Veseli (Veseli et al. 2025), makes functional capacity comparable across
metagenomic assemblies of very different community sizes.

The flag works only with an external-genomes file (-e), since the normalization is
meaningful for whole metagenomic assemblies (or long-read metagenomes) but not for
individual genomes or a genomes storage database. Assemblies whose SCGs are too
sparse to estimate a population count are reported as NA rather than dividing by
zero. All matrix outputs also lead with their descriptive columns (key, function,
accession) before the per-genome columns. Help page and CLI help included.

Author: jessika-fuessel

anvio/__init__.py

@@ -1329,6 +1329,23 @@ def TABULATE(table, header, numalign="right", max_width=0):
                      "anvi'o will only use the best hit for any gene that has multiple hits by default. But you can turn "
                      "that behavior off explicitly and show anvi'o who is the boss by using this flag."}
                 ),
+    'per-population-copy-number': (
+            ['--per-population-copy-number'],
+            {'default': False,
+             'action': 'store_true',
+             'help': "Request an additional output matrix that reports the per-population copy number (PPCN) of each "
+                     "function in every metagenomic assembly (or long-read metagenome). anvi'o divides the frequency of a "
+                     "function in a given assembly by the number of populations estimated to be present in that same "
+                     "assembly (estimated from single-copy core genes: the mode of the number of hits to the single-copy "
+                     "core genes of each domain, summed across Bacteria, Archaea, and Eukarya, just like "
+                     "`anvi-display-contigs-stats` does). This normalization is useful for metagenomic assemblies, where the "
+                     "raw frequency of a function scales with the number of populations in the assembly and is therefore "
+                     "not comparable across samples of vastly different community sizes. It follows the per-population copy "
+                     "number normalization introduced by Iva Veseli (Veseli et al. 2025). It requires that your "
+                     "contigs databases have been annotated with single-copy core genes via `anvi-run-hmms`, and it works "
+                     "only with an external-genomes file (`-e`), since population normalization is meaningful for whole "
+                     "metagenomic assemblies but not for individual genomes or a genomes storage database."}
+                ),
 
     'include-gc-identity-as-function': (
             ['--include-gc-identity-as-function'],

anvio/cli/gen_function_matrix_across_genomes.py

@@ -81,6 +81,15 @@ def get_args():
                                 "given function accross genomes. If you set a value, those functions that occur in less number "
                                 "of genomes will be excluded."), type=int)
 
+    groupD = parser.add_argument_group('PER-POPULATION COPY NUMBER', "If you are working with metagenomic assemblies, the raw frequency of a "
+                                "function will scale with the number of populations in each assembly, which makes it difficult to "
+                                "compare functions across samples of different community sizes. If you are simply looking at "
+                                "presence/absence, mostly everything will be present everywhere when looking at metagenomic "
+                                "assemblies as a whole. To account for this, anvi'o can "
+                                "report an additional matrix in which function frequencies are normalized by the estimated number "
+                                "of populations in each metagenome (an approach known as per-population copy number, or PPCN).")
+    groupD.add_argument(*anvio.A('per-population-copy-number'), **anvio.K('per-population-copy-number'))
+
     groupE = parser.add_argument_group('GENES', "By default, anvi'o will look for genes in contigs databases that are identified "
                                 "by `pyrodigal-gv`. But if you have generated your contigs databse with external gene calls, or have "
                                 "otherwise used another gene caller than the default, you can explicitly ask anvi'o to use that "

anvio/docs/artifacts/functions-across-genomes-txt.md

@@ -2,58 +2,64 @@ A TAB-delimited output file that describes the distribution of functions across
 
 ## General format
 
-With this file type anvi'o may describe the presence-absence of individual functions across genomes, or the frequency of functions. The header of the last column in a %(functions-across-genomes-txt)s file will be the name of the function annotation source from which the gene functions were recovered.
+With this file type anvi'o may describe the presence-absence of individual functions across genomes, or the frequency of functions. The first column contains a unique `key` for each function, and the second column the annotation generated from the annotation source the user chose. The third column is optional and includes the accession ids. The remaining columns hold one value per (meta)genome.
 
 ### Frequency of function hits across genomes
 
 Here is an example output file that shows the frequency of functions across 9 _Bifidobacterium_ genomes:
 
-|**key**|**B_adolescentis_1_11**|**B_adolescentis_22L**|**B_adolescentis_6**|**B_lactis_Bl_04**|**B_lactis_CNCM_I_2494**|**B_lactis_DSM_10140**|**B_longum_JDM301**|**B_longum_KACC_91563**|**B_longum_NCIMB8809**|**COG20_FUNCTION**|
+|**key**|**COG20_FUNCTION**|**B_adolescentis_1_11**|**B_adolescentis_22L**|**B_adolescentis_6**|**B_lactis_Bl_04**|**B_lactis_CNCM_I_2494**|**B_lactis_DSM_10140**|**B_longum_JDM301**|**B_longum_KACC_91563**|**B_longum_NCIMB8809**|
 |:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|
-|`func_3874a4219ffa`|1|1|1|1|1|1|1|1|1|Chromosomal replication initiation ATPase DnaA (DnaA) (PDB:1L8Q)|
-|`func_1530552f6b61`|1|1|1|1|1|1|1|1|1|DNA polymerase III sliding clamp (beta) subunit, PCNA homolog (DnaN) (PDB:1JQJ)|
-|`func_060988dd6d3a`|1|1|1|1|1|1|1|1|1|Recombinational DNA repair ATPase RecF (RecF) (PDB:5Z67)|
-|`func_46477a763d38`|1|1|1|1|1|1|1|1|1|Predicted  nucleic acid-binding protein, contains Zn-ribbon domain (includes truncated derivatives)|
-|`func_184ca83fe067`|2|2|2|2|2|2|2|2|2|DNA gyrase/topoisomerase IV, subunit B (GyrB) (PDB:1EI1)|
-|`func_2ee67347f04b`|2|2|2|2|2|2|2|2|2|DNA gyrase/topoisomerase IV, subunit A (GyrA) (PDB:1SUU)|
-|`func_ff160914972f`|1|0|0|0|0|0|3|1|1|Predicted ATPase, archaeal AAA+ ATPase superfamily|
-|`func_598325c18ceb`|1|0|0|0|0|0|1|1|1|Molybdopterin or thiamine biosynthesis adenylyltransferase (ThiF) (PDB:1ZUD) (PUBMED:32239579)|
-|`func_8ce5ff84aa42`|14|13|13|10|10|10|22|17|16|Predicted arabinose efflux permease AraJ, MFS family (AraJ) (PDB:4LDS)|
-|`func_300e2c6e37e4`|1|1|1|1|1|1|1|1|1|Glutamate dehydrogenase/leucine dehydrogenase (GdhA) (PDB:1B3B) (PUBMED:24391520)|
-|`func_a666bc87a03f`|1|1|1|1|1|1|1|1|1|Large-conductance mechanosensitive channel (MscL) (PDB:2OAR)|
-|`func_71379d89f0c6`|1|1|1|1|1|1|1|1|1|UTP pyrophosphatase, metal-dependent hydrolase family (YgjP) (PDB:4JIU) (PUBMED:27941785)|
-|`func_772efbb7cb2e`|2|2|2|2|2|2|2|2|2|Hemolysin-related protein, contains CBS domains, UPF0053 family (TlyC) (PDB:2NQW)|
-|`func_5a5bf9735c8a`|1|1|1|1|1|1|1|1|1|Carbonic anhydrase (CynT) (PDB:1EKJ) (PUBMED:22081392)|
-|`func_1ffc26d0034f`|4|4|1|3|3|3|12|13|2|Transposase (or an inactivated derivative) (IS285)|
-|`func_0631d9fecae6`|2|1|1|1|1|1|2|1|1|Phosphoenolpyruvate carboxylase (Ppc) (PDB:1FIY)|
-|`func_62068e09f9f3`|2|1|1|0|0|0|1|1|1|Chromosome segregation ATPase Smc (Smc) (PDB:5XG3)|
-|`func_f4b8779f6ad8`|2|0|0|0|0|0|0|0|0|Peptidoglycan-binding (PGRP) domain of peptidoglycan hydrolases (PGRP) (PDB:4FET)|
-|`func_700652d0a2dd`|1|1|1|1|1|1|1|1|1|Uncharacterized membrane protein YjjP, DUF1212 family (YjjP)|
+|`func_3874a4219ffa`|Chromosomal replication initiation ATPase DnaA (DnaA) (PDB:1L8Q)|1|1|1|1|1|1|1|1|1|
+|`func_1530552f6b61`|DNA polymerase III sliding clamp (beta) subunit, PCNA homolog (DnaN) (PDB:1JQJ)|1|1|1|1|1|1|1|1|1|
+|`func_060988dd6d3a`|Recombinational DNA repair ATPase RecF (RecF) (PDB:5Z67)|1|1|1|1|1|1|1|1|1|
+|`func_46477a763d38`|Predicted  nucleic acid-binding protein, contains Zn-ribbon domain (includes truncated derivatives)|1|1|1|1|1|1|1|1|1|
+|`func_184ca83fe067`|DNA gyrase/topoisomerase IV, subunit B (GyrB) (PDB:1EI1)|2|2|2|2|2|2|2|2|2|
+|`func_2ee67347f04b`|DNA gyrase/topoisomerase IV, subunit A (GyrA) (PDB:1SUU)|2|2|2|2|2|2|2|2|2|
+|`func_ff160914972f`|Predicted ATPase, archaeal AAA+ ATPase superfamily|1|0|0|0|0|0|3|1|1|
+|`func_598325c18ceb`|Molybdopterin or thiamine biosynthesis adenylyltransferase (ThiF) (PDB:1ZUD) (PUBMED:32239579)|1|0|0|0|0|0|1|1|1|
+|`func_8ce5ff84aa42`|Predicted arabinose efflux permease AraJ, MFS family (AraJ) (PDB:4LDS)|14|13|13|10|10|10|22|17|16|
+|`func_300e2c6e37e4`|Glutamate dehydrogenase/leucine dehydrogenase (GdhA) (PDB:1B3B) (PUBMED:24391520)|1|1|1|1|1|1|1|1|1|
+|`func_a666bc87a03f`|Large-conductance mechanosensitive channel (MscL) (PDB:2OAR)|1|1|1|1|1|1|1|1|1|
+|`func_71379d89f0c6`|UTP pyrophosphatase, metal-dependent hydrolase family (YgjP) (PDB:4JIU) (PUBMED:27941785)|1|1|1|1|1|1|1|1|1|
+|`func_772efbb7cb2e`|Hemolysin-related protein, contains CBS domains, UPF0053 family (TlyC) (PDB:2NQW)|2|2|2|2|2|2|2|2|2|
+|`func_5a5bf9735c8a`|Carbonic anhydrase (CynT) (PDB:1EKJ) (PUBMED:22081392)|1|1|1|1|1|1|1|1|1|
+|`func_1ffc26d0034f`|Transposase (or an inactivated derivative) (IS285)|4|4|1|3|3|3|12|13|2|
+|`func_0631d9fecae6`|Phosphoenolpyruvate carboxylase (Ppc) (PDB:1FIY)|2|1|1|1|1|1|2|1|1|
+|`func_62068e09f9f3`|Chromosome segregation ATPase Smc (Smc) (PDB:5XG3)|2|1|1|0|0|0|1|1|1|
+|`func_f4b8779f6ad8`|Peptidoglycan-binding (PGRP) domain of peptidoglycan hydrolases (PGRP) (PDB:4FET)|2|0|0|0|0|0|0|0|0|
+|`func_700652d0a2dd`|Uncharacterized membrane protein YjjP, DUF1212 family (YjjP)|1|1|1|1|1|1|1|1|1|
 |(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|
 
 ### Presence/absence of function hits across genomes
 
 In contrast, here is the presence/absence report for the same data:
 
-|**key**|**B_adolescentis_1_11**|**B_adolescentis_22L**|**B_adolescentis_6**|**B_lactis_Bl_04**|**B_lactis_CNCM_I_2494**|**B_lactis_DSM_10140**|**B_longum_JDM301**|**B_longum_KACC_91563**|**B_longum_NCIMB8809**|**COG20_FUNCTION**|
+|**key**|**COG20_FUNCTION**|**B_adolescentis_1_11**|**B_adolescentis_22L**|**B_adolescentis_6**|**B_lactis_Bl_04**|**B_lactis_CNCM_I_2494**|**B_lactis_DSM_10140**|**B_longum_JDM301**|**B_longum_KACC_91563**|**B_longum_NCIMB8809**|
 |:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|
-|`func_3874a4219ffa`|1|1|1|1|1|1|1|1|1|Chromosomal replication initiation ATPase DnaA (DnaA) (PDB:1L8Q)|
-|`func_1530552f6b61`|1|1|1|1|1|1|1|1|1|DNA polymerase III sliding clamp (beta) subunit, PCNA homolog (DnaN) (PDB:1JQJ)|
-|`func_060988dd6d3a`|1|1|1|1|1|1|1|1|1|Recombinational DNA repair ATPase RecF (RecF) (PDB:5Z67)|
-|`func_46477a763d38`|1|1|1|1|1|1|1|1|1|Predicted  nucleic acid-binding protein, contains Zn-ribbon domain (includes truncated derivatives)|
-|`func_184ca83fe067`|1|1|1|1|1|1|1|1|1|DNA gyrase/topoisomerase IV, subunit B (GyrB) (PDB:1EI1)|
-|`func_2ee67347f04b`|1|1|1|1|1|1|1|1|1|DNA gyrase/topoisomerase IV, subunit A (GyrA) (PDB:1SUU)|
-|`func_ff160914972f`|1|0|0|0|0|0|1|1|1|Predicted ATPase, archaeal AAA+ ATPase superfamily|
-|`func_598325c18ceb`|1|0|0|0|0|0|1|1|1|Molybdopterin or thiamine biosynthesis adenylyltransferase (ThiF) (PDB:1ZUD) (PUBMED:32239579)|
-|`func_8ce5ff84aa42`|1|1|1|1|1|1|1|1|1|Predicted arabinose efflux permease AraJ, MFS family (AraJ) (PDB:4LDS)|
-|`func_300e2c6e37e4`|1|1|1|1|1|1|1|1|1|Glutamate dehydrogenase/leucine dehydrogenase (GdhA) (PDB:1B3B) (PUBMED:24391520)|
-|`func_a666bc87a03f`|1|1|1|1|1|1|1|1|1|Large-conductance mechanosensitive channel (MscL) (PDB:2OAR)|
-|`func_71379d89f0c6`|1|1|1|1|1|1|1|1|1|UTP pyrophosphatase, metal-dependent hydrolase family (YgjP) (PDB:4JIU) (PUBMED:27941785)|
-|`func_772efbb7cb2e`|1|1|1|1|1|1|1|1|1|Hemolysin-related protein, contains CBS domains, UPF0053 family (TlyC) (PDB:2NQW)|
-|`func_5a5bf9735c8a`|1|1|1|1|1|1|1|1|1|Carbonic anhydrase (CynT) (PDB:1EKJ) (PUBMED:22081392)|
-|`func_1ffc26d0034f`|1|1|1|1|1|1|1|1|1|Transposase (or an inactivated derivative) (IS285)|
-|`func_0631d9fecae6`|1|1|1|1|1|1|1|1|1|Phosphoenolpyruvate carboxylase (Ppc) (PDB:1FIY)|
-|`func_62068e09f9f3`|1|1|1|0|0|0|1|1|1|Chromosome segregation ATPase Smc (Smc) (PDB:5XG3)|
-|`func_f4b8779f6ad8`|1|0|0|0|0|0|0|0|0|Peptidoglycan-binding (PGRP) domain of peptidoglycan hydrolases (PGRP) (PDB:4FET)|
-|`func_700652d0a2dd`|1|1|1|1|1|1|1|1|1|Uncharacterized membrane protein YjjP, DUF1212 family (YjjP)|
+|`func_3874a4219ffa`|Chromosomal replication initiation ATPase DnaA (DnaA) (PDB:1L8Q)|1|1|1|1|1|1|1|1|1|
+|`func_1530552f6b61`|DNA polymerase III sliding clamp (beta) subunit, PCNA homolog (DnaN) (PDB:1JQJ)|1|1|1|1|1|1|1|1|1|
+|`func_060988dd6d3a`|Recombinational DNA repair ATPase RecF (RecF) (PDB:5Z67)|1|1|1|1|1|1|1|1|1|
+|`func_46477a763d38`|Predicted  nucleic acid-binding protein, contains Zn-ribbon domain (includes truncated derivatives)|1|1|1|1|1|1|1|1|1|
+|`func_184ca83fe067`|DNA gyrase/topoisomerase IV, subunit B (GyrB) (PDB:1EI1)|1|1|1|1|1|1|1|1|1|
+|`func_2ee67347f04b`|DNA gyrase/topoisomerase IV, subunit A (GyrA) (PDB:1SUU)|1|1|1|1|1|1|1|1|1|
+|`func_ff160914972f`|Predicted ATPase, archaeal AAA+ ATPase superfamily|1|0|0|0|0|0|1|1|1|
+|`func_598325c18ceb`|Molybdopterin or thiamine biosynthesis adenylyltransferase (ThiF) (PDB:1ZUD) (PUBMED:32239579)|1|0|0|0|0|0|1|1|1|
+|`func_8ce5ff84aa42`|Predicted arabinose efflux permease AraJ, MFS family (AraJ) (PDB:4LDS)|1|1|1|1|1|1|1|1|1|
+|`func_300e2c6e37e4`|Glutamate dehydrogenase/leucine dehydrogenase (GdhA) (PDB:1B3B) (PUBMED:24391520)|1|1|1|1|1|1|1|1|1|
+|`func_a666bc87a03f`|Large-conductance mechanosensitive channel (MscL) (PDB:2OAR)|1|1|1|1|1|1|1|1|1|
+|`func_71379d89f0c6`|UTP pyrophosphatase, metal-dependent hydrolase family (YgjP) (PDB:4JIU) (PUBMED:27941785)|1|1|1|1|1|1|1|1|1|
+|`func_772efbb7cb2e`|Hemolysin-related protein, contains CBS domains, UPF0053 family (TlyC) (PDB:2NQW)|1|1|1|1|1|1|1|1|1|
+|`func_5a5bf9735c8a`|Carbonic anhydrase (CynT) (PDB:1EKJ) (PUBMED:22081392)|1|1|1|1|1|1|1|1|1|
+|`func_1ffc26d0034f`|Transposase (or an inactivated derivative) (IS285)|1|1|1|1|1|1|1|1|1|
+|`func_0631d9fecae6`|Phosphoenolpyruvate carboxylase (Ppc) (PDB:1FIY)|1|1|1|1|1|1|1|1|1|
+|`func_62068e09f9f3`|Chromosome segregation ATPase Smc (Smc) (PDB:5XG3)|1|1|1|0|0|0|1|1|1|
+|`func_f4b8779f6ad8`|Peptidoglycan-binding (PGRP) domain of peptidoglycan hydrolases (PGRP) (PDB:4FET)|1|0|0|0|0|0|0|0|0|
+|`func_700652d0a2dd`|Uncharacterized membrane protein YjjP, DUF1212 family (YjjP)|1|1|1|1|1|1|1|1|1|
 |(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|(...)|
+
+### Per-population copy number of functions across metagenomes
+
+When the program that generates this file is run with the `--per-population-copy-number` flag, anvi'o additionally reports the per-population copy number (PPCN) of each function. In this file, the frequency of a function in a given metagenome is divided by the number of populations estimated to be present in that metagenome (which anvi'o estimates from single-copy core genes). This normalization makes the values comparable across metagenomes that describe communities of different sizes.
+
+The layout is identical to the frequency file, except that the values are per-population copy numbers rather than raw counts. If anvi'o was unable to estimate the number of populations for a given metagenome, its values in this file will be reported as `NA`.

anvio/docs/programs/anvi-script-gen-function-matrix-across-genomes.md

@@ -20,10 +20,33 @@ The command above will result in two files in your work directory, both of which
 * MY-GENOMES-FREQUENCY.txt
 * MY-GENOMES-PRESENCE-ABSENCE.txt
 
+In each of these files, the first columns describe each function (a unique `key`, the function name, and optionally the function accession ids), and the remaining columns hold one value per genome.
+
 {:.notice}
 You can always learn about which functions are in a given %(contigs-db)s using the program %(anvi-db-info)s.
 
-Alternatively you can run it with a %(groups-txt)s that associates sets of genomes with distinct groups,
+## Per-population copy number normalization for metagenomic assemblies
+
+If we want to get an idea of differences in functional capacity across different metagenomic assemblies (or long-read sequence metagenomes) but, for some good reasons, we do not have MAGs from these assemblies, or we want to make sure that we make use of all the sequence data we have and not only those reads that are used to reconstruct genomes, we can't just look at the distribution of functions across genomes, because all functions and metabolic pathways will most likely occur nearly everywhere, at least in one population, and population numbers may differ dramatically across samples, and just counting the occurrence of a given function would not provide ecologically meaningful insights. To overcome this, Iva Veseli introduced the [per-population copy number](https://elifesciences.org/reviewed-preprints/89862) normalization. Based on the same principle, adding the flag `--per-population-copy-number` allows you to normalize individual functional annotations within a metagenomic assembly using the SCG-based estimate of population numbers within the sample.
+
+{{ codestart }}
+anvi-script-gen-function-matrix-across-genomes -e %(external-genomes)s \
+                                               --annotation-source COG20_FUNCTION \
+                                               --output-file-prefix MY-METAGENOMES \
+                                               --per-population-copy-number
+{{ codestop }}
+
+Adding the flag generates an additional output file of type %(functions-across-genomes-txt)s:
+
+* MY-METAGENOMES-PER-POPULATION-COPY-NUMBER.txt
+
+By dividing the frequency of each function in a given metagenomic assembly by the number of populations estimated to be present in that same assembly based on counts of single-copy core genes (SCGs) in each %(contigs-db)s: for each domain-specific SCG set, anvi'o takes the mode of the number of hits across all SCGs, and sums these per-domain estimates across Bacteria, Archaea, and Eukarya.
+
+The script relies on SCG annotations generated by %(anvi-run-hmms)s, and you need an %(external-genomes)s file to let it know which assemblies to compare.
+
+A per-population copy number is an *average number of copies of a function per population*: a value near 1 means the function tends to occur once in every population, a value below 1 means it is present in only a fraction of the populations, and a value above 1 means it tends to occur in multiple copies per population (as is common for, e.g., transposases). If a given metagenome has single-copy core genes that are too sparse to yield a reliable estimate, anvi'o cannot compute its number of populations, and the per-population copy number values for that assembly will be reported as `NA`.
+
+Alternatively, you can run it with a %(groups-txt)s that associates sets of genomes with distinct groups,
 
 {{ codestart }}
 anvi-script-gen-function-matrix-across-genomes -i %(internal-genomes)s \