# How to quantify similarity of genomes and find differences in set of S aureus genomes?

I have around 500 annotated proteomes of different bacterial strains and would like to quantify their similarity (or difference). I found gt genomediff from genometools gives me some scores that I can use to generate nice clusters, but I am not sure whether that tools really works. The fasta-files that I use contain multiple sequences.

I ran some tests and it looks ok.

1_reference.fna:
>1
TAAGTTACT
>2
TAAGTTACA

2_eq_to_ref.fna:
>1
TAAGTTACT
>2
TAAGTTACA

3_tags_diff.fna:
TAAGTTACT
>2asdfasdfffa
TAAGTTACA

4_orde_diff.fna:
>2
TAAGTTACA
>1
TAAGTTACT

>1
TAAGTTACT
>2
TAAGTTACA
>2
TTACA

6_point_mut.fna:
>1
AAAGTTACT
>2
TAAGTTACA

7_different.fna:
>1
TAAGTTACT
ATTACCTAA
>2
AAAAAAAAA


Then:

gt genomediff --indexname test *fna
7
1_reference.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
2_eq_to_ref.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
3_tags_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
4_orde_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
5_add_subse.fna 0.206969        0.206969        0.206969        0.206969        0.000000        0.212596        1.349206
6_point_mut.fna 0.199527        0.199527        0.199527        0.199527        0.212596        0.000000        0.569180
7_different.fna 0.794496        0.794496        0.794496        0.794496        1.349206        0.569180        0.000000


I am running the clustering on the difference matrix calculated with genomediff: "These distances are Jukes-Cantor corrected divergence between the pairs of genomes, that is, the number of mutations per base between them."

Currently, we are studying S aureus. We genomes are assembled genomes (three different methods). My guess is that the sequences from the plasmids are present. Furthermore, we do have drug-resistance measured in culture. So, we will be able to compare the genomes and the resistances.

• Can help here (I think), but need more info. Could you explain why you talk about "proteins" but present nucleotide data. What do the fasta sequences represent and what is your goal of analysis? Also which bacteria, pathogens?
– M__
Feb 13 '19 at 4:23
• Arguably, depending on what the goal of using said differences is for, it would be far more useful to infer a phylogeny of the bacterial strains you're analyzing and consider differences within the context of their evolutionary relationships. Feb 13 '19 at 14:34
• Thats what I was thinking too. In pathogens its quite common to track microbial drug-resistance and that tends to be presence/absence stuff.
– M__
Feb 13 '19 at 15:54
• The talk about proteins was an error, but I also have annotated proteomes generated from the dna sequence files. From here I would do hierarchical clustering. What tool would you use to generate a phylogenetic tree? First of all I want to quantify the difference between the strains to see whether drug-resistant or other classes of strains cluster together. Feb 13 '19 at 17:02
• Ah! right now we have a question. Its a while since I've done anything on microbial drug resistance so I'll need to collect my thoughts.
– M__
Feb 13 '19 at 22:08

... stuff of original post deleted. On second thoughts what you might be doing is a template assembly of your genomes. It is a possible interpretation of the 1.2. fasta sequences above (i.e. 1 is template).

A microbial genomics professional would advise also performing a de novo assembly, particularly if you are interested in presence and absence. The reason is that if a gene is present in your query and absent in your template, you will miss it. Again it all depends on what bacteria are being assessed, some are more prone to "genetic islands" than others.

You need a collaborator beyond that.

For the onlookers here bacterial genetic behaviour is very different to eukaryotes and what they get up, switching DNA etc.., to would appear bizzare from a higher eukaryotic world.

You mentioned a comparative analysis of 500 strains: I have worked on this in MRSA. Note, the bacteria is important considering the approach you will adopt.

Anyway, you want a single aligned file to produce a phylogeny, Bayes or likelihood. This is a model of point mutations. Bootstrapping doesn't really help because of SNP differences between isolates.

There is a complex problem. Generic phylogeny can hit the buffers in my opinion because often:

• low numbers of SNP differences between isolates.
• The error across the genome is poorly defined
• Multi-clonal infection ignored and your isolates not cloned.
• The other issue is the QC of the genome, I've been fairly astonished by the variation (it affects the tree)

Put it all together and you can get a nice tree but the topology of a given MLST is in my opinion unlikely to be correct. A referee may not bother with this, but they might.

You then map the phenotype against the phylogeny (the tips of the tree) and look for clusters. HOWEVER, there are problems:

• Drug-resistant e.g. mec cassette (methicillin-resistant ) are often on plasmids and these are lost during isolation. Albeit they can integrate onto the genome. So an isolate can be drug resistance but you fail to find the gene.
• Drug-resistance transmission is unlikely to cluster, its too quick, so you don't see nice tight clusters against a theoretically perfect tree
• The best approach is drug-resistance in culture.

Hierarchical clustering, has many meanings. It is used in presence/absence data, I demonstrated the method failed (on bacteria) based on cluster analysis. It is currently revised and used in unsupervised deep (or machine) learning prior a training method. The question then is what is what are you modelling? Presence/absence (of genes), point mutations, epidemiological data, drug-resistance/non-drug resistance?

My assessment is that given the potential complexity of these biological scenarios and the quality of the question you should seek formal collaboration and this is before considering what bacteria are being assessed. Some of epidemic drug-resistent bacteria rip up this rule book.

• I am running the clustering on the difference matrix calculated with genomediff. "These distances are Jukes-Cantor corrected divergence between the pairs of genomes, that is, the number of mutations per base between them." Currently, we are studying S aureus. We do have drug-resistance measured in culture. And we do have assembled genomes. My guess is that the sequences from the plasmids are present. Feb 14 '19 at 16:58
• JC is a primitive algorithm, particularly using cluster analysis ... a definite no no. Acually its quite rubbish and if I reviewed the ms I'd politely return it. On other issue is that I was authoritively informed that plasmids for MRSA (i.e. staph, i.e. S aureus) could not be guarrenteed and can be lost in isolation. E.g. the mec cassette may be present, but could be lost. I'm not a microbial biologist so wouldn't have a clue about plasmid retention, but it was the clear advise I got from exactly the same question, i.e. from people who are supposed to know.
– M__
Feb 14 '19 at 17:04
• Anyway I think we're done here ... at least get a negative rep for this thread (at least so far :-O ) ... Oh I see "S aureus" is mentioned in the title, my bad, ... the good thing is that you can use phylogeny on this bacteria, i.e. the theoretical model (which is what I do) matches the biology of the bacterial genetics ...
– M__
Feb 14 '19 at 17:07
• Michael, do you have a publication for your SA project? It seems to be very relevant for us. Feb 14 '19 at 17:15
• I do, but its not an area I work on and it didn't float my boat, so I prefer to remain anonymous. There are lots of publications on MRSA you can follow. I personally don't think the area is very well developed in terms of modeling brilliance, but one of Peacock's studies stands out.
– M__
Feb 14 '19 at 17:23