<p>Extract desired gene sequences using <a href="https://blast.ncbi.nlm.nih.gov/Blast.cgi" rel="nofollow noreferrer">standalone blast</a> </p>

<p>Simply provide a reference database with your desired output. 
Set up your command and away you go. You can set the search up with a for loop for a batch of sequences. Command may look like </p>

<pre><code>for f in *.fasta; do
   f=$(basename $f .fasta)
   blastn \
   -outfmt "6 sseqid qseq %" \
   -query $f.fasta \
   -subject reference.fna \
   &gt; out/$f.fas
done
</code></pre>

<p>Watch the output as blast will spit out the gene as detected in + or - sense. If you want to only gather positive sense use <code>-strand</code>option. 
The default output I have here is tab output which requires a few sed commands to make into fasta.</p>

<pre><code>sed -i \
-e 's/\s*$//g' \
-e 's/^/&gt;/g'  \
-e 's/\s\+/\n/g' \
*.fas 
</code></pre>

<p>Online alignment servers are an easy way to align small datasets e.g. <a href="http://www.ebi.ac.uk/Tools/msa/" rel="nofollow noreferrer">EBI</a></p>

CC BY-SA 3.0
2017-07-05T15:17:32.013Z
It might be worth mentioning that this approach work for species with no splicing (so it should be fine for the OP), but isn't very useful for eukaryotic organisms. Also, why would you run multiple blasts in a loop like this? Why not give a multifasta input file to blast and run it once?
terdon
CC BY-SA 3.0
2017-07-05T15:29:19.183Z
In my experience the blast out uses the gene names from the reference. and does not append the original genome name. I go about a few QC/ parsing bits and then combine the sequences.
AudileF
CC BY-SA 3.0
2017-07-05T15:30:27.070Z
I'm not sure what my "desired gene sequences" are supposed to be. What is my query?
Daniel Harris
CC BY-SA 3.0
2017-07-05T15:34:37.623Z
Your query is the genome sequences you have $f will insert the name for each in a loop. your desired gene sequences will be in your reference file. This could be the 16s or MLST gene sequences. Whatever you want to build a phylogeny of.
AudileF
CC BY-SA 3.0
2017-07-05T15:36:44.963Z
My question is more how can I find those MLST gene sequences in the first place?
Daniel Harris
CC BY-SA 3.0
2017-07-05T15:40:04.210Z
A quick google brought this up http://bigsdb.pasteur.fr/klebsiella/klebsiella.html . Have a look about there is usually a file full of genes and alleles. In my experience one (allele) set of the 7 genes should suffice.
AudileF
CC BY-SA 3.0
2017-07-05T15:48:13.480Z
@DanielHarris take one MLST gene and use that as a query to find the rest.
terdon
CC BY-SA 3.0
2017-07-05T15:49:18.653Z
@AudileF blast can take a multifasta file as a query and the results for each sequence in the input file will be shown in a separate `Query=` section. The names are indeed the names taken from the input file, but that shouldn't be a problem as long as they're all unique.
terdon
CC BY-SA 3.0
2017-07-05T15:57:38.043Z
So use one of these MLST genes such as from this link http://bigsdb.pasteur.fr/perl/bigsdb/bigsdb.pl?db=pubmlst_klebsiella_seqdef_public&page=downloadAlleles as a query. Two questions now... Do I use the entire locus file, such as gapA has many >gapA_#'s or just one >gapA_1? Second, how does the output of blast then be used to build a phylogenetic tree? Blast output is far from multiple alignment format I think.
Daniel Harris
CC BY-SA 3.0
2017-07-05T20:01:56.630Z
@terdon must admit I never know. Does this work in the case of contig assemblies?
AudileF
CC BY-SA 3.0
2017-07-05T20:08:51.593Z
@DanielHarris 1) use one copy of each of the 7 genes. 2) you are correct the output can be quite jumbled. I usually use the `-strand`option to gather the + and - sense sequences. Split them by taxa name. Order sequences in each file by name. remove fasta headers and space. Create new fasta header of file name. Then concatenate them all. For simplicity you could use a online msa server like on ebi.  Long winded and confusing. Im not a master computational biologist terdon might have a better suggestion.
AudileF
CC BY-SA 3.0
2017-07-05T14:54:52.480Z
2017-07-05T20:13:16.533Z
2017-07-05T20:13:16.533Z
AudileF
<p>There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the <a href="https://sanger-pathogens.github.io/Roary/" rel="nofollow noreferrer">Roary</a> website:</p>

<pre><code>for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff
</code></pre>

<p>Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses</p>

CC BY-SA 3.0
2017-07-05T15:02:38.310Z
Roary is my go to but you have no say in what come out in the alignment. If you want to select your gene content use blast as I suggest. Concatenate the output and use an alignment software.
AudileF
CC BY-SA 3.0
2017-07-05T20:39:10.287Z
@terdon I added a code snippet to do the analysis. It can be used to isolate genes (plural). If the OP just wants to build a phylogenetic tree, there's no need to isolate specific genes. With almost 500 genomes from a single species there's unlikely to me much resolution from analyses of a single gene or from a handful of MLST genes.
heathobrien
CC BY-SA 3.0
2017-07-05T14:50:16.667Z
2017-07-05T20:39:21.993Z
2017-07-05T20:39:21.993Z
heathobrien
<p>I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA [or] all ribosomal-protein-coding genes" <a href="https://www.biostars.org/p/1930/" rel="noreferrer">https://www.biostars.org/p/1930/</a>. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).</p>

CC BY-SA 3.0
2017-07-05T14:47:29.650Z
Hi Daniel, I would be cautious to build a phylogeny on 16s rRNA sequences as informative variability is essentially nil. I would suggest a good starting point being concatenated MLST genes, and work from there. You have the whole genome seems wrong to focus on the least variable gene.
AudileF
CC BY-SA 3.0
2017-07-05T15:15:42.953Z
Yes. Using 16S makes sense for distant species but not for closely related ones. it is extremely likely that there will be next to no variability at all.
terdon
CC BY-SA 3.0
2017-07-05T21:44:35.447Z
If you haven't yet annotated the genomes, do so with Prokka as suggested in some of the other answers. Once you've got that, you can easily pull gene features out of a genbank by name using biopython, once you've identified a few sequences you want to use as the basis for your sequence typing.
Joe Healey
<gene><phylogenetics><phylogeny>
How to isolate genes from whole genomes for phylogenetic tree analysis?
# How to isolate genes from whole genomes for phylogenetic tree analysis?

(Asked by: username_0 on 2017-07-05)

I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA [or] all ribosomal-protein-coding genes" <https://www.biostars.org/p/1930/>. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).

**Tags:** gene, phylogenetics, phylogeny
> **Comments:**
> **username_1**: Hi Daniel, I would be cautious to build a phylogeny on 16s rRNA sequences as informative variability is essentially nil. I would suggest a good starting point being concatenated MLST genes, and work from there. You have the whole genome seems wrong to focus on the least variable gene.  
> **username_2**: Yes. Using 16S makes sense for distant species but not for closely related ones. it is extremely likely that there will be next to no variability at all.  
> **username_3**: If you haven't yet annotated the genomes, do so with Prokka as suggested in some of the other answers. Once you've got that, you can easily pull gene features out of a genbank by name using biopython, once you've identified a few sequences you want to use as the basis for your sequence typing.  

---

## 2 answers

### Answer by username_1 on 2017-07-05. Score: 3

Extract desired gene sequences using standalone blast 
Simply provide a reference database with your desired output. 
Set up your command and away you go. You can set the search up with a for loop for a batch of sequences. Command may look like 
for f in *.fasta; do
   f=$(basename $f .fasta)
   blastn \
   -outfmt "6 sseqid qseq %" \
   -query $f.fasta \
   -subject reference.fna \
   > out/$f.fas
done

Watch the output as blast will spit out the gene as detected in + or - sense. If you want to only gather positive sense use -strandoption. 
The default output I have here is tab output which requires a few sed commands to make into fasta.
sed -i \
-e 's/\s*$//g' \
-e 's/^/>/g'  \
-e 's/\s\+/\n/g' \
*.fas 

Online alignment servers are an easy way to align small datasets e.g. EBI

> **Comments:**
> **username_2**: It might be worth mentioning that this approach work for species with no splicing (so it should be fine for the OP), but isn't very useful for eukaryotic organisms. Also, why would you run multiple blasts in a loop like this? Why not give a multifasta input file to blast and run it once?  
> **username_1**: In my experience the blast out uses the gene names from the reference. and does not append the original genome name. I go about a few QC/ parsing bits and then combine the sequences.  
> **username_0**: I'm not sure what my "desired gene sequences" are supposed to be. What is my query?  
> **username_1**: Your query is the genome sequences you have $f will insert the name for each in a loop. your desired gene sequences will be in your reference file. This could be the 16s or MLST gene sequences. Whatever you want to build a phylogeny of.  
> **username_0**: My question is more how can I find those MLST gene sequences in the first place?  
> **username_1**: A quick google brought this up http://bigsdb.pasteur.fr/klebsiella/klebsiella.html . Have a look about there is usually a file full of genes and alleles. In my experience one (allele) set of the 7 genes should suffice.  
> **username_2**: @DanielHarris take one MLST gene and use that as a query to find the rest.  
> **username_2**: @AudileF blast can take a multifasta file as a query and the results for each sequence in the input file will be shown in a separate `Query=` section. The names are indeed the names taken from the input file, but that shouldn't be a problem as long as they're all unique.  
> **username_0**: So use one of these MLST genes such as from this link http://bigsdb.pasteur.fr/perl/bigsdb/bigsdb.pl?db=pubmlst_klebsiella_seqdef_public&page=downloadAlleles as a query. Two questions now... Do I use the entire locus file, such as gapA has many >gapA_#'s or just one >gapA_1? Second, how does the output of blast then be used to build a phylogenetic tree? Blast output is far from multiple alignment format I think.  
> **username_1**: @terdon must admit I never know. Does this work in the case of contig assemblies?  
> **username_1**: @DanielHarris 1) use one copy of each of the 7 genes. 2) you are correct the output can be quite jumbled. I usually use the `-strand`option to gather the + and - sense sequences. Split them by taxa name. Order sequences in each file by name. remove fasta headers and space. Create new fasta header of file name. Then concatenate them all. For simplicity you could use a online msa server like on ebi.  Long winded and confusing. Im not a master computational biologist terdon might have a better suggestion.  

---

### Answer by username_4 on 2017-07-05. Score: 3

There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:
for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses

> **Comments:**
> **username_1**: Roary is my go to but you have no say in what come out in the alignment. If you want to select your gene content use blast as I suggest. Concatenate the output and use an alignment software.  
> **username_4**: @terdon I added a code snippet to do the analysis. It can be used to isolate genes (plural). If the OP just wants to build a phylogenetic tree, there's no need to isolate specific genes. With almost 500 genomes from a single species there's unlikely to me much resolution from analyses of a single gene or from a handful of MLST genes.  

---

