<p>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</p>
<hr />
<p>It called a combined alignment, or concatenated genes.</p>
<p>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.</p>
<p>You take your two alignments and test them for equal rates via IQTREE</p>
<pre><code>iqtree2 -s example.phy -p example.nex --symtest-only
</code></pre>
<p>Testing for a combined alignment is also described in my past answer <a href="https://bioinformatics.stackexchange.com/questions/19659/single-or-partitioned-model-for-supermatrix-tree-inference/19667#19667">here</a>: The previous question was can you combined an rRNA gene with a protein gene and the answer is ...</p>
<blockquote>
<p>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.</p>
</blockquote>
<hr />
<p>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 <a href="https://bioinformatics.stackexchange.com/questions/19814/is-there-an-equivalent-of-entrez-queries-for-local-blast">here</a>:</p>
<p>If you loads of fasta seqs you'll need to make a database...</p>
<pre><code>makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title &quot;Laura 700 - 1000&quot; -dbtype nucl
</code></pre>
<p>Then do this ...</p>
<pre><code>blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt
</code></pre>
<p>However if its just reference vs. query then you can just use the above command.</p>
<hr />
<p>The alignment is best done via a specialist alignment package Muscle Super5 and that is described in my answer <a href="https://bioinformatics.stackexchange.com/questions/18899/sequence-alignment-with-muscle">here</a>. This is done after your obtained the homologous regions in your reference genomes.</p>
<blockquote>
<p>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.</p>
</blockquote>

CC BY-SA 4.0
2022-10-26T08:53:11.783Z
Can you really get meaningful results using just two proteins, each aligned to its own reference genome, as input? Surely any tree built from such limited data will be unreliable, right? Why would you even build a tree if all you have are two proteins? Also, blast really isn't a good tool for this unless we know the OP is working on species with no introns, and in any case, they would need tblastn if they start from proteins.
terdon
CC BY-SA 4.0
2022-10-26T10:49:51.790Z
@terdon a tree needs 4-taxa ... so 3 reference genomes and the sequence of interest. Low reliability ... no ancient history, most tree methods moved on, its not affected by the sample size - absolutely not. Blast, despite being slow, is the work behind most published phylogenies.  .. introns .. it doesn't matter most quality tree building methods (IQTREE for example), ignore indels (introns are excluded) - why is a separate question -- thats a long explanation -- but based on their mutation matrix (indels are absent).
M__
CC BY-SA 4.0
2022-10-26T11:07:35.047Z
Edit comment above ... introns are excluded as phylogenetic information if they result in indels, things like Beast, IQTREE, RAxML all ignore indels. Thus even if they is a single indel for 10000 sequences that entire column is excluded for that particular nucleotide/amino acid.
M__
CC BY-SA 4.0
2022-10-26T11:38:59.867Z
OK, but how would you even build a tree using just two alignments, of two sequences each? As you said you need at least 4 taxa, and if (big if) I understand the question correctly, the OP has at best 2 species (no idea if they're the same taxon or not) and at worse just two genes from the same species.
terdon
CC BY-SA 4.0
2022-10-26T11:44:24.890Z
@terdon I edited the posted previously ...
M__
CC BY-SA 4.0
2022-10-26T11:47:45.767Z
Brilliant, thanks!
terdon
CC BY-SA 4.0
2022-10-25T22:19:25.220Z
2022-10-26T10:40:51.320Z
2022-10-26T10:40:51.320Z
M__
<p>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.</p>
<p>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 <a href="https://en.wikipedia.org/wiki/Molecular_clock" rel="nofollow noreferrer">molecular clock</a>). After that, you can do typical phylogenetic analysis.</p>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2022-10-26T17:14:25.147Z
Definitely agreed on 4 sequences. IMO an unrooted tree obviates the need to assign an outgroup. BTW mid-point rooting doesn't assume a global molecular clock, most tree building produce phylograms - non-clocks (Beast is an exception, options depending). UPGMA method assumes a clock.
M__
CC BY-SA 4.0
2022-10-27T04:03:40.810Z
@M__ I am repeating myself, so I don't think this discussion goes anywhere. I have deleted my comments and summarized them in the answer. Feel free to disagree but I am confident in my knowledge as well.
user172818
CC BY-SA 4.0
2022-10-27T15:47:23.160Z
Thanks @user172818 I've done likewise.
M__
CC BY-SA 4.0
2022-10-26T16:22:32.973Z
2022-10-27T03:57:24.797Z
2022-10-27T03:57:24.797Z
user172818
<p>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?</p>

CC BY-SA 4.0
2022-10-25T19:08:46.913Z
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Community
CC BY-SA 4.0
2022-10-26T08:55:01.317Z
What is the actual question you are trying to answer here? Are you trying to identify the gene coding for these proteins? Is the genome annotated? Are you working on a species that has introns? Alternative splicing?
terdon
<python><phylogenetics><phylogeny>
How can I align two novel sequences to a reference genome and build a phylogenetic tree?
# How can I align two novel sequences to a reference genome and build a phylogenetic tree?

(Asked by: username_0 on 2022-10-25)

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?

**Tags:** python, phylogenetics, phylogeny
> **Comments:**
> **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
> **username_2**: What is the actual question you are trying to answer here? Are you trying to identify the gene coding for these proteins? Is the genome annotated? Are you working on a species that has introns? Alternative splicing?  

---

## 2 answers

### Answer by username_3 on 2022-10-25. Score: 2

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.

> **Comments:**
> **username_2**: Can you really get meaningful results using just two proteins, each aligned to its own reference genome, as input? Surely any tree built from such limited data will be unreliable, right? Why would you even build a tree if all you have are two proteins? Also, blast really isn't a good tool for this unless we know the OP is working on species with no introns, and in any case, they would need tblastn if they start from proteins.  
> **username_3**: @terdon a tree needs 4-taxa ... so 3 reference genomes and the sequence of interest. Low reliability ... no ancient history, most tree methods moved on, its not affected by the sample size - absolutely not. Blast, despite being slow, is the work behind most published phylogenies.  .. introns .. it doesn't matter most quality tree building methods (IQTREE for example), ignore indels (introns are excluded) - why is a separate question -- thats a long explanation -- but based on their mutation matrix (indels are absent).  
> **username_3**: Edit comment above ... introns are excluded as phylogenetic information if they result in indels, things like Beast, IQTREE, RAxML all ignore indels. Thus even if they is a single indel for 10000 sequences that entire column is excluded for that particular nucleotide/amino acid.  
> **username_2**: OK, but how would you even build a tree using just two alignments, of two sequences each? As you said you need at least 4 taxa, and if (big if) I understand the question correctly, the OP has at best 2 species (no idea if they're the same taxon or not) and at worse just two genes from the same species.  
> **username_3**: @terdon I edited the posted previously ...  
> **username_2**: Brilliant, thanks!  

---

### Answer by username_4 on 2022-10-26. Score: 0

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](https://en.wikipedia.org/wiki/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.

> **Comments:**
> **username_3**: Definitely agreed on 4 sequences. IMO an unrooted tree obviates the need to assign an outgroup. BTW mid-point rooting doesn't assume a global molecular clock, most tree building produce phylograms - non-clocks (Beast is an exception, options depending). UPGMA method assumes a clock.  
> **username_4**: @M__ I am repeating myself, so I don't think this discussion goes anywhere. I have deleted my comments and summarized them in the answer. Feel free to disagree but I am confident in my knowledge as well.  
> **username_3**: Thanks @user172818 I've done likewise.  

---

