Repository: widdowquinn/pyani (language: python)

Issue:
ANIm should not be symmetric
#### Summary:
PyANI treats ANIm as symmetric, but  the order of genomes impacts this calculation.

#### Description:
I have two test genomes. The ANIm results produced by PyANI change depending on the order these genomes are processed. That is, the results depend on which genome is considered the "reference" and which the "query" when running nucmer.

#### Reproducible Steps:
1. Process the attached genomes with PyANI using `average_nucleotide_identity.py -m ANIm --workers 1`.
2. Swap the names of the genomes (i.e. genome1.fna becomes genome2.fna and genome2.fna becomes genome1.fna).
3. Run PyANI again as in step 1, but changing the output directory.
4. Compare the results and observe that the ANIm values are now different. 

As a concrete results, the number of similarity errors changes from 512 to 561 depending on the order of the genomes.

### Expected Output:
[test_genomes.tar.gz](https://github.com/widdowquinn/pyani/files/7883433/test_genomes.tar.gz)


The ANI-percentage_identity.tab, ANIm_similarity_errors.tab, and ANIm_alignment_lengths.tab files are non-symmetric matrices that reflect which genome was the reference and which genome was the query.

#### pyani Version:
v0.2.11

#### Operating System:
Ubuntu

Alignment coverage >1.0
#### Summary:
I am comparing a large number of genomes and occasionally the ANIm_alignment_coverage.tab file indicates a value >1.0. I would have though this value was bounded between [0, 1.0]. How does one interpret an alignment coverage > 1.0?


#### Description:
I am using PyANI v0.2.11 as follows:
`average_nucleotide_identity.py -i genomes/ -o anim -v -m ANIm`

I can narrow the issue down to a pair of viral genomes which are in the `genomes` directory.


#### Reproducible Steps:
I can provide the pair of viral genomes which cause this issue.

#### Current Output:
The `ANIm_alignment_coverage.tab` file indicates the alignment coverage between these pairs of genomes is 1.49 and 1.51, respectively.


#### Expected Output:
I would have though the alignment coverage would have a maximum value of 1.0. That is, any region of the query genome would align to at most one region in the target genome and thus it would be impossible to cover more than 100% of a genome.

#### pyani Version:
0.2.11



Discussion / hints:
Thanks @dparks1134 - apologies for the bug, and we'll deal with this as a priority.

Quick notes

- although it stems from the same misunderstanding of/assumptions about `mummer`'s operation, this is not the same issue as #340 and will probably not be fixed by the proposed changes to resolve that issue.
- the resolution currently appears to be to run `mummer` for the reciprocal pairwise comparisons.

L.
NOTE: Handling under branch `issue_373` for version 0.3.
I can add that ANIb produces the expect result for these two genomes with an alignment coverage of 0.98 and 0.99.
Thanks @dparks1134 - that's weird. If you could provide the input genomes that would be very helpful.

Cheers,

L.
FWIW I found (and fixed) a related bug while trying this out on a couple of (unrelated) virus genomes. I don't think it will fix your issue, though.
Zip file with the two genomes and PyANI ANIm results are attached.

[high_align_cov.zip](https://github.com/widdowquinn/pyani/files/7261637/high_align_cov.zip)

Thanks - I can confirm your ANIm output. Looking into *why*, just now.
Well - the issue appears to be that MUMmer is reporting two overlapping alignments. One alignment runs from position 85 to 37713 in the query; the other alignment runs from 17709 to 39253 in 0264574. Inspecting the alignment length output shows that the aligned sequence length is longer than either genome.

We use the `--mum` option to ensure that there are unique anchors for the NUCmer alignment, so this was not the behaviour I was expecting. I errored in my understanding of how MUMmer generates alignment output.

I think that what we should have been doing all along is *also* to use the `--noextend` option, by default. With this option we still get alignments at {85, 27720} and {27702, 39253} in the query and the ANI results are otherwise in line with ANIb. I suspect, though don't yet know for certain, that the default `--extend` option in NUCmer allows the anchors to extend into each other - and using `--mum --noextend` is needed to ensure that alignments are unique to both query and subject sequence (though in this example they *still* overlap by 18 bases!).

I feel quite bad about missing that. I'm not going to defend it by claiming that we're reimplementing the Rossello-Mora et al. implementation - by default we use `--mum` (not indicated in that paper) to try to ensure we don't align repeats. It's just a mistake.

I'm going to add a command-line option `--noextend` to pyani v0.2, so that we retain backwards compatibility by default. @baileythegreen - I think we should maybe make `--noextend` the default NUCmer behaviour for v0.3, and implement an `--extend` option for `pyani anim` that allows us to use `pyani compare` to see what the damage is likely to have been under the assumption above.

I'll make the changes and push a v0.2.12 release with an apology and explanation.

Many thanks for spotting this and bringing this to our attention @dparks1134 - this is a very important improvement that will impact directly on the next iteration of `pyani`. I appreciate the heads-up, and hope it doesn't put you off using the tool.

Cheers,

L.






@all-contributors please add @dparks1134 for bug
@widdowquinn 

I've put up [a pull request](https://github.com/widdowquinn/pyani/pull/341) to add @dparks1134! :tada:
Thanks for the detailed explanation. Nice to know what is going on. 
@widdowquinn I've created a new issue for this in version 3, and will prioritise that, since it probably qualifies as a bug.
Hi. I implemented the proposed `noextend` fix in v0.2.11. This does reduce the number of pairs with an alignment coverage >1, but this still occur in a few instances. The values are just above one in all cases I have observed so perhaps this is fine, but it would be good to ensure things aren't being double counted. I've attached an email pair of genomes.

[AF_bug_v2.zip](https://github.com/widdowquinn/pyani/files/7267941/AF_bug_v2.zip)
Thanks @dparks1134 - it does look like even `--noextend` cannot guarantee that no bases are double-counted in the MUMmer alignment.

I've been talking with @baileythegreen about potentially using interval graphs to ensure we don't double-count and, if you've found further examples of double-counting with `--noextend` already, it looks like we're going to have to do that *instead*.

I expect this hasn't been noticed already because the overlap size probably gets lost in the noise for bacterial genomes, and the values look plausible.
Hi @widdowquinn, do you know if it is sensible to run `dnadiff` and look at the generated 1-to-1 coordinates and/or alignment file? I know at least one person that has taken this approach, but I don't know enough about the MUMmer package to know if this is actually correct. 
Thanks for the suggestion - I've taken a look at `dnadiff` and it doesn't appear to do anything differently to what `pyani` already does - it runs a default `nucmer` search and then uses `delta-filter` to generate a 1:1 alignment (extension `.1coords` with `dnadiff`). It does some things _in addition_ like identifying SNPs and reporting M-to-M alignments, but they don't contribute to ANI calculations.

`pyani`'s current ANIm `nucmer` output and the `.1coords` output are identical. There are no useful options to `dnadiff` that I see; we can't, for instance, pass through either `--mum` or `--noextend` to `nucmer`.

`nucmer` default, followed by `delta-filter -1`:

```text
/Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0264574.fna /Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0266457.fna
NUCMER
>MGV_MGV-GENOME-0264574 MGV_MGV-GENOME-0266457 39253 39594
85 37713 1 37636 215 215 0
-3013
-24624
-1
-1
-1
-1
-1
0
17709 39253 17626 39176 7 7 0
-9994
-1
-1
-1
-1
-1
0
```

`dnadiff`'s `.1delta` output:

```text
/Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0264574.fna /Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0266457.fna
NUCMER
>MGV_MGV-GENOME-0264574 MGV_MGV-GENOME-0266457 39253 39594
85 37713 1 37636 215 215 0
-3013
-24624
-1
-1
-1
-1
-1
0
17709 39253 17626 39176 7 7 0
-9994
-1
-1
-1
-1
-1
0
```

`dnadiff's` `.report` - notice that it claims a total alignment of 59kbp for the 39kbp genomes:

```text
$ mkdir -p dnadiff_default && dnadiff -p dnadiff_default/high_align_test high_align_cov/MGV-GENOME-0264574.fna high_align_cov/MGV-GENOME-0266457.fna
[...]
/Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0264574.fna /Users/lpritc/Development/GitHub/pyani/issue_340/high_align_cov/MGV-GENOME-0266457.fna
NUCMER

                               [REF]                [QRY]
[Sequences]
TotalSeqs                          1                    1
AlignedSeqs               1(100.00%)           1(100.00%)
UnalignedSeqs               0(0.00%)             0(0.00%)

[Bases]
TotalBases                     39253                39594
AlignedBases           39169(99.79%)        39176(98.94%)
UnalignedBases             84(0.21%)           418(1.06%)

[Alignments]
1-to-1                             2                    2
TotalLength                    59174                59187
AvgLength                   29587.00             29593.50
AvgIdentity                    99.63                99.63
[...]
```

I'm surprised and disappointed that even the `--noextend` option in `nucmer` gives overlapping alignments. We _can_ account for this, but it will take some additional output filtering of the `nucmer` output prior to calculation. We're investigating how much of a difference using `--noextend` makes to our analyses (the tests seem to be nearly all within tolerance, so far), but it's clearly not going to be the actual solution, if we want to ensure that there are no overlaps between alignments.
Dose it make sense to follow the `AlignedBases` approach used by dnadiff? The answer certainly seems sensible in this case, but it is unclear to me exactly how this number is established.
Yes, the `AlignedBases` value looks to be a good place to start for coverage. For instance, `dnadiff` reports 39169 `AlignedBases` for the reference, which corresponds to the 39253 - 85 interval we'd expect when not double-counting the overlap. That helps towards improving the coverage/aligned fraction calculation.

Currently, the ANI identity calculation will be double-counting the overlap, and `AlignedBases` doesn't look to help with that. I think we may still need to deduplicate overlapping regions in order to give an accurate ANI value.
Largely a note for me and @baileythegreen after today's discussion...

It looks like we might be able to get everything we need from a combination of the `.1coords` and `.snps` files that `dnadiff` produces, so we can get the "right" answers for identity, coverage, etc.

However, for very large comparisons on a cluster we probably don't want `dnadiff` to generate that suite of multiple output files, so it might be best to do something like run `show-snps` and `show-coords` directly (catching each STDOUT) in a wrapper, and pass the results of the calculation to the database, in v0.3+ - though there may be a better solution.
My initial optimism may have been misplaced. There appear to be some disturbing differences in the _Caulobacter_ test set - especially wrt coverage.

[ANIm_alignment_coverage_noextend.pdf](https://github.com/widdowquinn/pyani/files/7514027/ANIm_alignment_coverage_noextend.pdf)
[ANIm_alignment_coverage.pdf](https://github.com/widdowquinn/pyani/files/7514028/ANIm_alignment_coverage.pdf)


Further investigation of the `dnadiff` and `nucmer` parsing is making me optimistic again. It's quite clear that `dnadiff` is *not* using `--noextend`, but is using something like a tweaked `--maxmatch`, and is then using the output of M-to-M matches (which we avoid), and processing them to remove overlaps to get the value in `AlignedBases`.

The difference in output between what we do and what `dnadiff` does turns out to be quite small for a bacterial test set.

Where the `dnadiff` `.report` output gives `AlignedBases = 1414384` (reference) or `1424179` (query), the 1:1 alignments are actually `1408533` and `1409770` in total length on reference and query genomes; the M:M alignments are `1442608` and `1443869` respectively.

I'd argue that we should be only using 1:1 alignments, so `AlignedBases` is not the right number to use.

However, we were not originally accounting for overlapping alignments. When we *do* account for them, `dnadiff` gives (look at `ref_aligned_size` and `query_aligned_size`):

```text
infname=PosixPath('2021-11-10/test_dnadiff_output/NC_002696_vs_NC_010338.1coords')
aln_ref=1408533, aln_query=1409770, cov_ref=35.050000000000125, cov_query=25.850000000000044, num_rows=1206, overlaps=11, overlap_len=17664
len(ref_tree)=1195, len(query_tree)=1201
ref_aligned_size=1401112, query_aligned_size=1408639

infname=PosixPath('2021-11-10/test_dnadiff_output/NC_002696_vs_NC_010338.mcoords')
aln_ref=1442608, aln_query=1443869, cov_ref=35.82000000000012, cov_query=26.40000000000006, num_rows=1275, overlaps=53, overlap_len=50359
len(ref_tree)=1222, len(query_tree)=1230
ref_aligned_size=1414420, query_aligned_size=1424194
```

and our standard `pyani` approach gives:

```text
infname=PosixPath('2021-11-10/test_nucmer_output/NC_002696_vs_NC_010338.fcoords')
aln_ref=1388340, aln_query=1389491, cov_ref=34.59000000000013, cov_query=25.51000000000005, num_rows=1173, overlaps=5, overlap_len=7845
len(ref_tree)=1168, len(query_tree)=1172
ref_aligned_size=1384814, query_aligned_size=1389479
```

With our alignment choices (including 1:1 match filtering) we align about 1.38Mbp where `dnadiff` aligns 1.41Mbp.

However, if we use `--maxmatch` rather than `--mum`, we get:

```text
infname=PosixPath('2021-11-10/test_nucmer_maxmatch_output/NC_002696_vs_NC_010338.fcoords')
aln_ref=1408533, aln_query=1409770, cov_ref=35.050000000000125, cov_query=25.850000000000044, num_rows=1206, overlaps=11, overlap_len=17664
len(ref_tree)=1195, len(query_tree)=1201
ref_aligned_size=1401112, query_aligned_size=1408639
```

and almost exactly recover the `dnadiff` output.

So, while our output differs, it appears to do so *mostly* because we use `--mum` rather than `--maxmatch`. As `--maxmatch` is an option intended to find inexact repeats and uses all anchors, whether unique or not, `--mum` uses only anchors that are unique in both reference and query. We do offer the `--maxmatch` option to users already if they don't agree with our default, so I think we might have this covered.

However, that's not the problem with issue #340 ;)

The issue with #340 is that we were not emulating the overlap removal that `dnadiff` employs when generating its `.report` file. Note that the raw alignment lengths, and the `.1coord`/`.mcoord` outputs have the same problem we did - this from the original virus comparison:

```text
[Bases]
TotalBases                     39253                39594
AlignedBases           39169(99.79%)        39176(98.94%)
UnalignedBases             84(0.21%)           418(1.06%)

[Alignments]
1-to-1                             2                    2
TotalLength                    59174                59187
AvgLength                   29587.00             29593.50
AvgIdentity                    99.63                99.63

M-to-M                             2                    2
TotalLength                    59174                59187
AvgLength                   29587.00             29593.50
AvgIdentity                    99.63                99.63
```

`AlignedBases` is great, but the alignment lengths are not.

We can fix this by running an `IntervalTree` over the alignment co-ordinates. I did that for the virus comparisons to get:

```text
infname=PosixPath('2021-11-10/virus_dnadiff_output/virus_dnadiff.1coords')
aln_ref=59174, aln_query=59187, cov_ref=150.75, cov_query=149.48, num_rows=2, overlaps=1, overlap_len=39169
len(ref_tree)=1, len(query_tree)=1
ref_aligned_size=39169, query_aligned_size=39176

infname=PosixPath('2021-11-10/virus_dnadiff_output/virus_dnadiff.mcoords')
aln_ref=59174, aln_query=59187, cov_ref=150.75, cov_query=149.48, num_rows=2, overlaps=1, overlap_len=39169
len(ref_tree)=1, len(query_tree)=1
ref_aligned_size=39169, query_aligned_size=39176

infname=PosixPath('2021-11-10/virus_nucmer_output/virus_nucmer.coords')
aln_ref=59174, aln_query=59187, cov_ref=150.75, cov_query=149.48, num_rows=2, overlaps=1, overlap_len=39169
len(ref_tree)=1, len(query_tree)=1
ref_aligned_size=39169, query_aligned_size=39176
```

and recovers `AlignedBases` exactly, whether we use `dnadiff` or our existing approach.
So, I think we should be generating a `.coords` file output and parsing it through an `IntervalTree` to remove/account for overlaps. That's the main change.

I think we already allow users to choose between `--maxmatch`/`--mum` so no change there.

I don't understand the process by which `AlignedBases` is generated exactly, but as it appears to be based on M:M alignments, I think it's possibly the wrong value for us to use, if we're trying to identify uniquely-homologous regions.

Once we correct for the overlaps, I think that our agreement with `dnadiff`/other methods should be pretty good and, except in unusual cases like the MAG viruses, I don't currently think it will affect `pyani` outputs to a large degree - though it *will* be an improvement.

Fix patch:
diff --git a/Makefile b/Makefile
index 90678b94..629ff1f0 100644
--- a/Makefile
+++ b/Makefile
@@ -47,7 +47,7 @@ walkthrough: clean_walkthrough
 	pyani report --runs -o C_blochmannia_ANIm/ --formats html excel stdout
 	pyani report --run_results 1 --formats html excel stdout -o C_blochmannia_ANIm/
 	pyani report --run_matrices 1 --formats html excel stdout -o C_blochmannia_ANIm/
-	pyani plot --formats png pdf --method seaborn -o C_blochmannia_ANIm --run_ids 1
+	pyani plot --formats png pdf --method seaborn -o C_blochmannia_ANIm --run_id 1
 	# pyani anib C_blochmannia C_blochmannia_ANIb \
     #     --name "C. blochmannia run 2" \
     #     --labels C_blochmannia/labels.txt --classes C_blochmannia/classes.txt
diff --git a/pyani/pyani_graphics/sns/__init__.py b/pyani/pyani_graphics/sns/__init__.py
index c8bcb284..fdcddea5 100644
--- a/pyani/pyani_graphics/sns/__init__.py
+++ b/pyani/pyani_graphics/sns/__init__.py
@@ -37,6 +37,7 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 # THE SOFTWARE.
 """Code to implement Seaborn graphics output for ANI analyses."""
+import warnings
 
 import matplotlib  # pylint: disable=C0411
 import pandas as pd
@@ -106,17 +107,37 @@ def get_clustermap(dfr, params, title=None, annot=True):
     :param title:  str, plot title
     :param annot:  Boolean, add text for cell values?
     """
-    fig = sns.clustermap(
-        dfr,
-        cmap=params.cmap,
-        vmin=params.vmin,
-        vmax=params.vmax,
-        col_colors=params.colorbar,
-        row_colors=params.colorbar,
-        figsize=(params.figsize, params.figsize),
-        linewidths=params.linewidths,
-        annot=annot,
-    )
+
+    # If we do not catch warnings here, then we often get the following warning:
+    #   ClusterWarning: scipy.cluster: The symmetric non-negative hollow
+    #   observation matrix looks suspiciously like an uncondensed distance matrix
+    # The usual solution would be to convert the array with
+    # scipy.spatial.distance.squareform(), but this requires that all values in
+    # the main diagonal are zero, which is not the case for ANI.
+    # As we know this is a (1-distance) matrix, we could just set the diagonal
+    # to zero and fudge it, but this is not a good solution. Instead, we suppress
+    # the warning in a context manager for this function call only, because we
+    # know the warning is not relevant.
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore",
+            message=(
+                "scipy.cluster: The symmetric non-negative "
+                "hollow observation matrix looks suspiciously like an "
+                "uncondensed distance matrix"
+            ),
+        )
+        fig = sns.clustermap(
+            dfr,
+            cmap=params.cmap,
+            vmin=params.vmin,
+            vmax=params.vmax,
+            col_colors=params.colorbar,
+            row_colors=params.colorbar,
+            figsize=(params.figsize, params.figsize),
+            linewidths=params.linewidths,
+            annot=annot,
+        )
 
     # add labels for each of the input genomes
     add_labels(fig, params)
diff --git a/pyani/scripts/average_nucleotide_identity.py b/pyani/scripts/average_nucleotide_identity.py
index 10fc7767..0fc0dcbe 100755
--- a/pyani/scripts/average_nucleotide_identity.py
+++ b/pyani/scripts/average_nucleotide_identity.py
@@ -812,7 +812,6 @@ def draw(args: Namespace, filestems: List[str], gformat: str) -> None:
         infilename = fullstem.with_suffix(".tab")
         dfm = pd.read_csv(infilename, index_col=0, sep="\t")
         logger.info("Writing heatmap to %s", outfilename)
-        print(args.labels, args.classes)
         params = pyani_graphics.Params(
             params_mpl(dfm)[filestem],
             pyani_tools.get_labels(args.labels, logger=logger),
diff --git a/pyani/scripts/parsers/plot_parser.py b/pyani/scripts/parsers/plot_parser.py
index 0aac8031..b35f8e01 100644
--- a/pyani/scripts/parsers/plot_parser.py
+++ b/pyani/scripts/parsers/plot_parser.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2016-2019
-# (c) University of Strathclyde 2019-2020
+# (c) University of Strathclyde 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -17,7 +17,7 @@
 # The MIT License
 #
 # Copyright (c) 2016-2019 The James Hutton Institute
-# Copyright (c) 2019-2020 University of Strathclyde
+# Copyright (c) 2019-2022 University of Strathclyde
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
diff --git a/pyani/scripts/parsers/report_parser.py b/pyani/scripts/parsers/report_parser.py
index 22767a79..c1ca01e1 100644
--- a/pyani/scripts/parsers/report_parser.py
+++ b/pyani/scripts/parsers/report_parser.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2016-2019
-# (c) University of Strathclyde 2019-2020
+# (c) University of Strathclyde 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -17,7 +17,7 @@
 # The MIT License
 #
 # Copyright (c) 2016-2019 The James Hutton Institute
-# Copyright (c) 2019-2020 University of Strathclyde
+# Copyright (c) 2019-2022 University of Strathclyde
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
diff --git a/pyani/scripts/subcommands/subcmd_plot.py b/pyani/scripts/subcommands/subcmd_plot.py
index d3682d1d..4e3f1ab0 100644
--- a/pyani/scripts/subcommands/subcmd_plot.py
+++ b/pyani/scripts/subcommands/subcmd_plot.py
@@ -194,6 +194,9 @@ def write_distribution(
             title=f"matrix_{matdata.name}_run{run_id}",
         )
 
+    # Be tidy with matplotlib caches
+    plt.close("all")
+
 
 def write_heatmap(
     run_id: int,
diff --git a/pyani/scripts/subcommands/subcmd_report.py b/pyani/scripts/subcommands/subcmd_report.py
index 82423893..453b29cc 100644
--- a/pyani/scripts/subcommands/subcmd_report.py
+++ b/pyani/scripts/subcommands/subcmd_report.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2017-2019
-# (c) University of Strathclyde 2019-2020
+# (c) University of Strathclyde 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -18,7 +18,7 @@
 # The MIT License
 #
 # Copyright (c) 2017-2019 The James Hutton Institute
-# Copyright (c) 2019-2020 University of Strathclyde
+# Copyright (c) 2019-2022 University of Strathclyde
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -40,6 +40,7 @@
 """Provides the report subcommand for pyani."""
 
 import logging
+import warnings
 
 from argparse import Namespace
 from typing import List, NamedTuple
@@ -154,12 +155,7 @@ def subcmd_report(args: Namespace) -> int:
             "genome label",
             "genome class",
         ]
-        report(
-            args,
-            session,
-            formats,
-            ReportParams("runs_genomes", statement, headers),
-        )
+        report(args, session, formats, ReportParams("runs_genomes", statement, headers))
 
     # Report table of all runs in which a genome is involved
     if args.show_genomes_runs:
@@ -196,12 +192,7 @@ def subcmd_report(args: Namespace) -> int:
             "method",
             "date run",
         ]
-        report(
-            args,
-            session,
-            formats,
-            ReportParams("genomes_runs", statement, headers),
-        )
+        report(args, session, formats, ReportParams("genomes_runs", statement, headers))
 
     # Report table of comparison results for the indicated runs
     if args.run_results:
@@ -234,6 +225,7 @@ def subcmd_report(args: Namespace) -> int:
                 .filter(Run.run_id == run_id)
                 .statement
             )
+            logger.debug("Results query: %s", statement)
             headers = [
                 "Comparison ID",
                 "Query ID",
@@ -300,12 +292,7 @@ def subcmd_report(args: Namespace) -> int:
     return 0
 
 
-def report(
-    args: Namespace,
-    session,
-    formats: List[str],
-    params: ReportParams,
-) -> None:
+def report(args: Namespace, session, formats: List[str], params: ReportParams) -> None:
     """Write tabular report of pyani runs from database.
 
     :param args:  Namespace of command-line arguments
@@ -319,7 +306,47 @@ def report(
     logger.debug(
         "Writing table of pyani %s from the database to %s.*", params.name, outfname
     )
-    data = pd.read_sql(params.statement, session.bind)
+
+    # With newer versions of SQLAlchemy, Pandas may throw a warning due to the composition
+    # of our statement including a Cartesian product, even though this is what we want:
+    #   SAWarning: SELECT statement has a cartesian product between FROM element(s)
+    #   "runs" and FROM element "genome_query".  Apply join condition(s) between each
+    #   element to resolve.
+    # We could use SQLAlchemy's true() function to force the join condition, but this has
+    # to be done from within Pandas, and is an issue for them to fix.
+    # We suppress the warning, instead.
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore",
+            message=(
+                "SELECT statement has a cartesian product between FROM "
+                'element\\(s\\) "runs" and FROM element "genome_query"'
+            ),
+        )
+        warnings.filterwarnings(
+            "ignore",
+            message=(
+                "SELECT statement has a cartesian product between FROM "
+                'element\\(s\\) "runs" and FROM element "comparisons"'
+            ),
+        )
+        warnings.filterwarnings(
+            "ignore",
+            message=(
+                "SELECT statement has a cartesian product between FROM "
+                'element\\(s\\) "runs" and FROM element "genome_subject"'
+            ),
+        )
+        warnings.filterwarnings(
+            "ignore",
+            message=(
+                "SELECT statement has a cartesian product between FROM "
+                'element\\(s\\) "genome_query", "genome_subject", "comparisons" '
+                'and FROM element "runs"'
+            ),
+        )
+        data = pd.read_sql(params.statement, session.bind)
+
     data.columns = params.headers
     pyani_report.write_dbtable(data, outfname, formats)
 
diff --git a/requirements.txt b/requirements.txt
index 15d97e59..62841ab0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,5 +10,5 @@ Pillow
 PyQt5
 scipy
 seaborn
-sqlalchemy==1.3.10
+sqlalchemy
 tqdm


Test patch:
diff --git a/tests/test_anib.py b/tests/test_anib.py
index 2f8d2380..f522a107 100644
--- a/tests/test_anib.py
+++ b/tests/test_anib.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2017-2019
-# (c) University of Strathclyde 2019-2020
+# (c) University of Strathclyde 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact: leighton.pritchard@strath.ac.uk
@@ -17,7 +17,7 @@
 # The MIT License
 #
 # Copyright (c) 2017-2019 The James Hutton Institute
-# Copyright (c) 2019-2020 University of Strathclyde
+# Copyright (c) 2019-2022 University of Strathclyde
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -49,7 +49,6 @@
 import pandas as pd
 import pytest  # noqa: F401  # pylint: disable=unused-import
 import unittest
-import unittest.mock as mock
 
 from pandas.testing import assert_frame_equal
 
@@ -367,8 +366,8 @@ def test_parse_legacy_blastdir(anib_output_dir):
         anib_output_dir.legacyblastdir, orglengths, fraglengths, mode="ANIblastall"
     )
     assert_frame_equal(
-        result.percentage_identity.sort_index(1).sort_index(),
-        anib_output_dir.legacyblastresult.sort_index(1).sort_index(),
+        result.percentage_identity.sort_index(axis=1).sort_index(),
+        anib_output_dir.legacyblastresult.sort_index(axis=1).sort_index(),
     )
 
 
@@ -380,8 +379,8 @@ def test_parse_blastdir(anib_output_dir):
         anib_output_dir.blastdir, orglengths, fraglengths, mode="ANIb"
     )
     assert_frame_equal(
-        result.percentage_identity.sort_index(1).sort_index(),
-        anib_output_dir.blastresult.sort_index(1).sort_index(),
+        result.percentage_identity.sort_index(axis=1).sort_index(),
+        anib_output_dir.blastresult.sort_index(axis=1).sort_index(),
     )
 
 
diff --git a/tests/test_anim.py b/tests/test_anim.py
index d3c40621..baf5b5b9 100644
--- a/tests/test_anim.py
+++ b/tests/test_anim.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2017-2019
-# (c) University of Strathclyde 2019-2021
+# (c) University of Strathclyde 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -18,7 +18,7 @@
 # The MIT License
 #
 # Copyright (c) 2017-2019 The James Hutton Institute
-# Copyright (c) 2019-2021 University of Strathclyde
+# Copyright (c) 2019-2022 University of Strathclyde
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -43,7 +43,7 @@
 
 pytest -v
 """
-import sys
+
 from pathlib import Path
 from typing import List, NamedTuple, Tuple
 
@@ -51,7 +51,7 @@
 import pytest
 import unittest
 
-from pandas.util.testing import assert_frame_equal
+from pandas.testing import assert_frame_equal
 
 from pyani import anim, pyani_files
 
@@ -208,8 +208,8 @@ def test_deltadir_parsing(delta_output_dir):
     orglengths = pyani_files.get_sequence_lengths(seqfiles)
     result = anim.process_deltadir(delta_output_dir.deltadir, orglengths)
     assert_frame_equal(
-        result.percentage_identity.sort_index(1).sort_index(),
-        delta_output_dir.deltaresult.sort_index(1).sort_index(),
+        result.percentage_identity.sort_index(axis=1).sort_index(),
+        delta_output_dir.deltaresult.sort_index(axis=1).sort_index(),
     )
 
 
@@ -251,17 +251,14 @@ def test_mummer_single(tmp_path, path_file_two):
         path_file_two[0], path_file_two[1], outdir=tmp_path
     )
     dir_nucmer = tmp_path / "nucmer_output"
+    outprefix = (
+        dir_nucmer
+        / str(path_file_two[0].stem)
+        / str(path_file_two[0].stem + "_vs_" + path_file_two[1].stem)
+    )
     expected = (
-        (
-            "nucmer --mum -p "
-            f"{dir_nucmer / str(path_file_two[0].stem) / str(path_file_two[0].stem + '_vs_' + path_file_two[1].stem)} "
-            f"{path_file_two[0]} {path_file_two[1]}"
-        ),
-        (
-            "delta_filter_wrapper.py delta-filter -1 "
-            f"{dir_nucmer / str(path_file_two[0].stem ) / str(path_file_two[0].stem + '_vs_' + path_file_two[1].stem + '.delta')} "
-            f"{dir_nucmer / str(path_file_two[0].stem ) / str(path_file_two[0].stem + '_vs_' + path_file_two[1].stem + '.filter')}"
-        ),
+        f"nucmer --mum -p {outprefix} {path_file_two[0]} {path_file_two[1]}",
+        f"delta_filter_wrapper.py delta-filter -1 {outprefix}.delta {outprefix}.filter",
     )
     assert cmds == expected
 
diff --git a/tests/test_subcmd_05_report.py b/tests/test_subcmd_05_report.py
index 49bf0419..d30096ce 100644
--- a/tests/test_subcmd_05_report.py
+++ b/tests/test_subcmd_05_report.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2017-2019
-# (c) The University of Strathclude 2019
+# (c) The University of Strathclude 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -19,7 +19,7 @@
 # The MIT License
 #
 # Copyright (c) 2017-2018 The James Hutton Institute
-# (c) The University of Strathclude 2019
+# (c) The University of Strathclude 2019-2022
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
diff --git a/tests/test_tetra.py b/tests/test_tetra.py
index 752dcd47..13507c9d 100644
--- a/tests/test_tetra.py
+++ b/tests/test_tetra.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # (c) The James Hutton Institute 2017-2019
-# (c) The University of Strathclude 2019-2020
+# (c) The University of Strathclude 2019-2022
 # Author: Leighton Pritchard
 #
 # Contact:
@@ -19,7 +19,7 @@
 # The MIT License
 #
 # Copyright (c) 2017-2019 The James Hutton Institute
-# (c) The University of Strathclude 2019-2020
+# (c) The University of Strathclude 2019-2022
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -46,13 +46,10 @@
 """
 
 import json
-import unittest
-
-from pathlib import Path
 
 import pandas as pd
 
-from pandas.util.testing import assert_frame_equal
+from pandas.testing import assert_frame_equal
 
 from pyani.tetra import (
     calculate_correlations,