I am trying to do alignment among orthologs and paralogs of a gene but I encountered some problems. I have came up with a workflow, but the sequences I extracted are not correct. They are of hugely varying length, and the alignment is very poor.

I have attached the workflow below, could anyone suggest if I am missing out on anything or if I am interpreting the files wrongly? Really appreciate your help.

My query is: with a list of target genes and bacteria genomes (.fna), for each gene, extract all sequences of it from all genomes and put together into one file for alignment


1 Answer 1


The answer is very simple, just align from amino acids. Amino acids are not subject to saturation - which is what you are encountering and is smashing up the nucleotide alignment. Amino acid will result in a high quality alignment. Thus you think the sequences are wrong, I think the sequences are simply divergent.

The other issues of variation in length - that depends on the how you've assembled the genes and exactly what project this is. Sinorhizobium is a pretty specific bacteria, so it's pretty clear what the biology is - its quite possibly legume agriculture. However, the precise molecular study is not so clear. Basically, what I'm saying is 'variations' in length are not necessarily 'errors' if this is metagenomics.

You could extract the amino acids from your fna file. Or ... use gff format to get the amino acids via gffread -y protein.fa that sort of thing.

The following might not make much sense ... but hope it helps ... to obtain the DNA alignment ... you convert the amino acid alignment back to a nucleotide alignment by transferring the alignment from amino acids on to the nucleotide data using Emboss' tranalign, amongst many other programs.

The tranalign command is:

tranalign -asequence [nucleotide filename (fasta)] -bsequence [aligned protein filename (fasta)] -outseq 

tranalign is also available online, albeit the program is nothing special and its just as easy to write the code.

You need a bit of code to translate the nucleotides to amino acids [I use Biopython translate()] or just use seqotron of MEGAX translate function (if you've a Mac) for a GUI.

Warning If you use this approach the taxa names (strain names) KOs etc ... need to be identical between amino acid and nucleotide files AND the order needs to be identical too. This can be a headache - but if you encounter errors here simply post back as a separate question and I'll forward the code to get around this 'glitch'.

I hope this all sort of makes sense.

HOWEVER, if you want to assess whether KOfamscan is giving the correct results (which is what you are questioning) you could use

KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold available

This is a webserver available here that will do the same thing and should give the same results. The authors consider it an upgrade, but I've not used it.

There are ALOT of other approaches to this problem BTW and all have over-reliance on KEGG, whilst it is the most comprehensive database available, is not perfect. KEGG per se omits upwards of half of all genes. Nobody else needs to know that but as long as you are aware of this thats fine.

  • 1
    $\begingroup$ Hi, thank you so much for your detailed answer. Sorry that I missed out another part of the workflow in my original question, but my plan was to translate the .fna to .faa first, do alignment on .faa, and do codon alignment to get the aligned nucleotides. I suspect that the start and end position on the .gff file is not correct, such that I cannot extract the correct nucleotide sequence. I have to get the nucleotide sequence for selection pressure analysis. $\endgroup$
    – Lily Li
    Commented Apr 1, 2023 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.