Answers.Body: <p>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.</p>
  <p>The other issues of variation in length - that depends on the how you've assembled the genes and exactly what project this is. <em>Sinorhizobium</em> 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.</p>
  <p>You could extract the amino acids from your <code>fna</code> file. Or ... use gff format to get the amino acids via <code>gffread -y protein.fa</code> that sort of thing.</p>
  <p>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' <code>tranalign</code>, amongst many other programs.</p>
  <p>The tranalign command is:</p>
  <pre><code>tranalign -asequence [nucleotide filename (fasta)] -bsequence [aligned protein filename (fasta)] -outseq 
  </code></pre>
  <p>tranalign is also available <a href="https://www.bioinformatics.nl/cgi-bin/emboss/tranalign" rel="nofollow noreferrer">online</a>, albeit the program is nothing special and its just as easy to write the code.</p>
  <p>You need a bit of code to translate the nucleotides to amino acids [I use Biopython <code>translate()</code>] or just use <code>seqotron</code> of <code>MEGAX</code> translate function (if you've a Mac) for a GUI.</p>
  <p><strong>Warning</strong> 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'.</p>
  <p>I hope this all sort of makes sense.</p>
  <hr />
  <p>HOWEVER, if you want to assess whether KOfamscan is giving the correct results (which is what you are questioning) you could use</p>
  <p><a href="https://academic.oup.com/bioinformatics/article/36/7/2251/5631907" rel="nofollow noreferrer">KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold available</a></p>
  <p>This is a webserver available <a href="https://www.genome.jp/tools/kofamkoala/" rel="nofollow noreferrer">here</a> that will do the same thing and should give the same results. The authors consider it an upgrade, but I've not used it.</p>
  <hr />
  <p>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 <em>per se</em> omits upwards of half of all genes. Nobody else needs to know that but as long as you are aware of this thats fine.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-01T21:21:07.827Z
Answers.Text: 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.
Answers.UserDisplayName: Lily Li
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-01T11:34:52.130Z
Answers.LastActivityDate: 2023-04-01T14:58:23.313Z
Answers.LastEditDate: 2023-04-01T14:58:23.313Z
Answers.UserDisplayName: M__
Body: <p>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.</p>
  <p>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.</p>
  <p>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
  <a href="https://i.sstatic.net/piLpj.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/piLpj.jpg" alt="" /></a></p>
  
Tags: <gff3><kegg>
Title: How to extract sequences of a specific gene from genomes, with genomes' .fna, .gff and results from KofamScan?
ThreadText: # How to extract sequences of a specific gene from genomes, with genomes' .fna, .gff and results from KofamScan?
  
  (Asked by: username_0 on 2023-04-01)
  
  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
  [![](https://i.sstatic.net/piLpj.jpg)](https://i.sstatic.net/piLpj.jpg)
  
  **Tags:** gff3, kegg
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-04-01. Score: 1 | Accepted answer
  
  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](https://www.bioinformatics.nl/cgi-bin/emboss/tranalign), 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](https://academic.oup.com/bioinformatics/article/36/7/2251/5631907)
  
  This is a webserver available [here](https://www.genome.jp/tools/kofamkoala/) 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.
  
  > **Comments:**
  > **username_0**: 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.  
  
  ---
  
  