How to determine the phylogeny of a specific set of genes

Question

My goal is to understand the phylogeny of genes involved in bacterial restriction-modification systems, in psychrophilic bacteria.

To do so here is what I have done.

1. Downloaded the protein sequences of all RM genes in psychrophilic bacteria from REBASE. I also downloaded a set of genomes I'd like understand the phylogeny of from NCBI.
2. Clustered these genes using usearch's uclust method.
3. Aligned the resulting centroids using MUSCLE5.
4. Built an HMM from the aligned centroids using HMMER3.
5. Ran hmmsearch to get hits from my HMM on each genome (AA).

I am unsure how to proceed now. I now a phylogenetic tree requires an alignment, but I'm not sure where my HMM hits fit-in. Am I meant to use hmmalign on the output of hmmsearch and then feed the resulting alignment to IQ-TREE or similar?

Thank you for your help as I learn bioinformatics.

PS: Reading through the forum it seems like MOLE-BLAST does exactly this, however I have too many genes and their length is too long to use the MOLE-BLAST as hosted by NCBI at least.

I should mention I set the identity threshold of search to 50% for proteins. Also, I got this method from this paper which achieved something similar (https://pubmed.ncbi.nlm.nih.gov/30893937/).

Answer 1

Your strategy is fine but the order is wrong.

What I would do is:

• 2 -> checked by 4 and 6

Run usearch and check the output of that in HMMER3. You can then use this to kick out stuff from usearch, that are not real homologues. Thats what HMMER3, hmmeralign (I dunno hmmersearch), are doing (below). It's also possible the HMMER approaches could replace usearch.

Once you have identified your homologues with confidence you can then align each separate group via muscle5

Your way

What you are doing is using HMMER3 to assess the quality of the alignment. What could be a problem is feeding usearch straight into muscle5 ... a non-homologous protein will get you a bad alignment. The bad alignment will not give you a good tree. Both aligners and trees need nice homologues, otherwise they don't work. The errors between aligner (muscle5) and tree are also synergistic. Thus if that bad sequences isn't removed by hmmer stuff then alignment errors compound onto sequence errors in the bad sequence giving a really bad tree.

You could of course align, HMMER3, filter stuff then align again (thats much better a bit complicated).

I don't know how good usearch is - but you are relying on it alot, supposing it gives you orthologues at <25% identity that could be any protein. Basically you are feeding that stuff straight into muscle5 and if HMMER3 doesn't trap it that will go straight onto a tree. That could be a bad alignment and that would be a bad tree.

Summary

1. What you have to understand is phylogenetic trees are very sensitive to bad sequences, i.e. stuff that is not true homologues. The aligner will also be sensitive. There will also be bad synergistic relationships between the two.

2. HMMER3 and is associated methods (hmmeralign) are trying to clean stuff up before that stage, rather than simply relying on usearch. Prof Eddie (HMMER and pfam[?] author) has an excellent reputation in this area - not in tree building but identifying protein families and thats what you are doing here.

3. hmmeralign doesn't replace muscle5, but it will provide a very good assessment of the protein families around a given seed. It's going to give you the homologues. That output would then go to muscle5 to get homologous sites.

4. The tree needs both homologues and homologous sites.

The whole hmmer project has a really good reputation in the area of protein family identification and thats why its been recommended.

Probable simple answer I think the answer to the question is you're are confusing hmmeralign with muscle5. They are complimentary but distinct. Again hmmeralign identifies homologues and the muscle5 designates homologous sites.

---

Answering the comments

50% cut-off

It looks like usearch clusters according to a %age identity threshold.

There are three weaknesses:

1. Gene duplication can take place about 50% ... it depends whether thats cool, usually it isn't, but it's your call.
2. There can be complete change of function >50% identity and this is about the definition of a homologue. If you are using a functional definition its not a homologue
3. Homologues can have <<50% identity and thats very common certainly in viruses I don't know about bacteria.

hmmersearch is going to rank your data according to the hmmer algorithms. It's going to be like blast,but all hmmer algorithms are focused around defining protein families (it's a better approach than blast for this application). So it provides another option to assess the output of usearch. Thus if there a gene duplication for example you should be able to spot it (it will be a sharp break in the continuous output). This will help you assess you 50% cut-off.

hmmeralign would be quite cool, because you can take a representation of usearch and then search that via hmmeralign and see if you get a similar output. It is extremely important for assessing 50% cut-off, does the homology drop below 50%? hmmeralign will tell you, but it must see the complete original database (not the stuff extracted from usearch. You can cross check that with hmmersearch to assess where hmmeralign is making the cut-off. So hmmeralign will independently (or usearch) will find you homologues and hmmersearch will give you an insight about how it's done it.