pyani compare

[baileythegreen]

@widdowquinn I've been thinking about how this should work, so I can start writing it. Here are the things I think:

• it works by saying pyani compare --dbpath <db file> --runs <M>,<N> -o <outdir> [other stuff]
• M and N are two runs that either use different methods, or different parameters for the same method
• output includes a summary report with:
  • the commands used to generate each run, or a breakdown of same, indicating parameters
  • some plots, probably, either in the report, or as separate files
    • are the plots of the differences, or plots for both, just side-by-side?
  • things like ANI value, coverage, et cetera
  • metrics of similarity? I don't know which ones we may want to include

[widdowquinn]

it works by saying pyani compare --dbpath <db file> --runs <M>,<N> -o <outdir> [other stuff]

That's what I was also thinking. I don't have a strong preference for --runs <M>,<N> or --run1 <M> --run2 <N>. The first is more readily expandable to --runs <M>,<N>,<P>,... which I like the idea of (with more than two runs, each pairwise run comparison could get its own subdirectory, and there could be a summary subdirectory with the mean/std dev of each measure across all the runs being compared)

In the [other stuff] there's the usual -v -l logfile.txt etc. but we could use --nocoverage, --nosimerrs toggles to help the user control the output where those values are uninteresting or don't make sense.

M and N are two runs that either use different methods, or different parameters for the same method

I'm happy to let the user provide the same run twice if they like ;)

some plots, probably, either in the report, or as separate files

The outputs I'm thinking of are:

• tables/heatmaps of the difference in measurement for each of the runs, for each pairwise run comparison. Sometimes we will not have a comparator value (e.g. TETRA doesn't give coverage), so single cells would be NA and we'd have to handle that somehow.
• for pairwise run comparisons, a scatterplot for each metric of run1 value vs run2 value, with the 1:1 line shown - this will allow identification of bias/diagnostics
• a summary distribution of differences in scores between methods for each pairwise genome comparison; in "good" cases this will be a Normal distribution with mean zero (indicating essentially "random" differences) - in others this may indicate skew/bias/be useful for diagnostics
• a summary report as you describe: the pyani compare command used; what was being compared/command lines; a summary (possibly count) of which genome comparisons were in common and which not; general summary stats for each run - highest ANI, lowest coverage, that sort of thing

[baileythegreen]

@widdowquinn For the tables/heatmaps of differences in measures:

• should these be sided differences?
• just absolute differences?
• should the values be normalised?

I'm reasonably confident (subject to some testing) that I have the database logic of all of this sorted out. I'm now working on molding the inputs for the plotting functions.

[baileythegreen]

Pros/cons of the options, as I see them:

• Sided:
  • the choice of direction is inherently arbitrary and it makes the whole thing asymmetric.
• Absolute:
  • any information about which run scored 'higher' for a given comparison is lost.
• Normalisation:
  • looking at a single set of runs becomes less informative, but
  • more useful for looking at several pairs of runs.

[widdowquinn]

By "sided" do you mean "signed"?

I think that we need to report both absolute and relative differences. I also think there are a number of ways to avoid the disadvantages you mention.

For instance, suppose we compare two runs (A and B) and for argument's sake look at the percentage identity difference. The things I would care about are:

• what are the absolute differences in the identities calculated by each method?
• what are the relative differences in the identities calculated by each method (expressed as a percentage of the original value)?

Using a two-colour, divergent scale means we can represent both value and sign, relative to a reference point. Clearly, one of the runs must be chosen as the reference. I'd expect the first run to be the default reference.

With this scheme we have choices for the colour scales. We could, for instance:

• allow "positive" and "negative" sign colours to have maximum intensity at whatever the maximum value is (equivalent to a mode of normalisation, except in terms of the number that is reported)
• define arbitrary (maybe choosable?) "positive" and "negative" sign colour maximum identities.

I'd prefer the second, especially if by the time the colour reaches its maximum intensity it's already an indicator of a worrying difference. Making the limits parametrisable by the user would enable consistent colouring across multiple pairwise comparisons.

I'd also be interested in seeing not just an all-vs-all comparison heatmap, but the ability to choose a set of comparisons, such that a smaller heatmap could be produced, allowing us to focus in on some particular comparisons, e.g. (each [ ] being a cell in the heatmap):

method: 1 2 3 4 (implicitly against a reference)
A v B [ ] [ ] [ ] [ ]
B v A [ ] [ ] [ ] [ ]
A v C [ ] [ ] [ ] [ ]
C v A [ ] [ ] [ ] [ ]
B v C [ ] [ ] [ ] [ ]
C v B [ ] [ ] [ ] [ ]

[baileythegreen]

By "sided" do you mean "signed"?

I mean being consistent about the side of the mathematical operator on which the runs in the various subtraction operations reside (always doing run_A_<measure> - run_B_<measure>; which I suppose is not different from 'signed'.

I'm going to put work on colour choices and subsetting comparisons off until after I've gotten it happily feeding stuff into the plotting commands. The subsetting, in particular, will take some thinking.

[widdowquinn]

We could provide (or allow the user to select) additional output that is essentially abs(a - b), so order doesn't matter and there's just a magnitude of difference and not a sign - but that wouldn't work for a / b as a percentage or ratio.

[baileythegreen]

@widdowquinn I have been thinking, (and trying to name things unambiguously in the code for this subcommand), and the words 'compare' and 'comparison' are pretty overloaded in pyani, generally.

I propose this feature instead be called pyani contrast.

[widdowquinn]

I don't know. We refer to "sequence comparison" and "comparative genomics" rather than "sequence contrast" and "contrasting genomics"… what's the positive case for "contrast" beyond "too much 'compare'"?

[baileythegreen]

(Reminder that there is a difference between adding a new comment on a discussion and replying to an existing one.)

I don't know that I have a positive case for 'contrast', except that as a term it's relevant and easily understood.

However, I wouldn't consider the terms you cite as particularly relevant here. The purpose of the subcommand in question is not to run comparative genomics or to do sequence comparisons; it isn't dealing with the sequences at all, beyond their IDs in the database. It is supposed to show the similarities and differences that exist between two different pyani run results. To compare and contrast them. I see those two actions as equal parts of what the subcommand should do, and I don't see a compelling reason for preferring 'compare'.

[widdowquinn]

When I'm asked to "compare" things, I take it to mean that I'm looking for similarities and differences; when I'm asked to "contrast" I assume that I'm looking primarily for differences.

I think that's how these terms are generally understood and that this is why we have "sequence comparison" and not "sequence contrast" as the terminology that the field has settled on. I think it's also why "contrast" is used in the technical sense, in statistics, too. I wasn't implying that pyani compare does sequence comparisons.