How can I align two novel sequences to a reference genome and build a phylogenetic tree?

Question

I have two new amino acid gene sequences that I want to align to a reference genome. Is it a good idea to have a combined phylogenetic tree? what will this type phylogenetic analysis called?

Answer 1

Edit: just to point out a tree needs 4 taxa (sequences), so rather than one reference genome, 3 reference genomes. Thus for Blast you might need to produce a database (using 3 reference genomes) before Blasting. There are other Blast approaches to obtaining a large meaningful sample size ... thats a separate question

---

It called a combined alignment, or concatenated genes.

The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a single alignment and make a tree from them. The theory goes back a long but the modern test is in IQTREE.

You take your two alignments and test them for equal rates via IQTREE

iqtree2 -s example.phy -p example.nex --symtest-only

Testing for a combined alignment is also described in my past answer here: The previous question was can you combined an rRNA gene with a protein gene and the answer is ...

Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.

---

The first part of the question was mapping the gene to a reference genome. Anything will do that ... local Blast works and using that is described in my answer to this question here:

If you loads of fasta seqs you'll need to make a database...

makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

Then do this ...

blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

However if its just reference vs. query then you can just use the above command.

---

The alignment is best done via a specialist alignment package Muscle Super5 and that is described in my answer here. This is done after your obtained the homologous regions in your reference genomes.

HMM Super5 is a hidden Markov model (HMM) where gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

Answer 2

I would strongly recommend to find an outgroup or at least a fourth sequence. Building a phylogenetic tree with 3 sequences is doable but is tricky and error prone.

If you still want to do a 3-sequence analysis, you can align the two homologous proteins to the reference genome with exonerate, spaln or miniprot, extract the reference subsequence in the alignment, do a proper multiple-sequence alignment (MSA) with muscle or MAFFT, build a tree and put the root at the middle point of the longest branch (this assumes a molecular clock). After that, you can do typical phylogenetic analysis.

Whether the above works depends on the divergence of your sequences. If they are closely related and one sequence is a clear outgroup to the other two, the procedure could give you something, though having more sequences is still recommended.

EDIT in response to the comments: midpoint rooting assumes an approximate molecular clock across the whole tree. Severe violation of the assumption may lead to a wrong root. This approach is not very reliable due to the clock assumption.