Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:17:32.013Z
Answers.Text: 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?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:29:19.183Z
Answers.Text: 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.
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:30:27.070Z
Answers.Text: I'm not sure what my "desired gene sequences" are supposed to be. What is my query?
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:34:37.623Z
Answers.Text: 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.
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:36:44.963Z
Answers.Text: My question is more how can I find those MLST gene sequences in the first place?
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:40:04.210Z
Answers.Text: 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.
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:48:13.480Z
Answers.Text: @DanielHarris take one MLST gene and use that as a query to find the rest.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:49:18.653Z
Answers.Text: @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.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:57:38.043Z
Answers.Text: 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.
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T20:01:56.630Z
Answers.Text: @terdon must admit I never know. Does this work in the case of contig assemblies?
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T20:08:51.593Z
Answers.Text: @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.
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T14:54:52.480Z
Answers.LastActivityDate: 2017-07-05T20:13:16.533Z
Answers.LastEditDate: 2017-07-05T20:13:16.533Z
Answers.UserDisplayName: AudileF
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T15:02:38.310Z
Answers.Text: 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.
Answers.UserDisplayName: AudileF
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T20:39:10.287Z
Answers.Text: @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.
Answers.UserDisplayName: heathobrien
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T14:50:16.667Z
Answers.LastActivityDate: 2017-07-05T20:39:21.993Z
Answers.LastEditDate: 2017-07-05T20:39:21.993Z
Answers.UserDisplayName: heathobrien
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-05T14:47:29.650Z
Comments.Text: 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.
Comments.UserDisplayName: AudileF
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-05T15:15:42.953Z
Comments.Text: 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-05T21:44:35.447Z
Comments.Text: 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.
Comments.UserDisplayName: Joe Healey
Tags: <gene><phylogenetics><phylogeny>
Title: How to isolate genes from whole genomes for phylogenetic tree analysis?
ThreadText: # 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 username_0, 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](https://blast.ncbi.nlm.nih.gov/Blast.cgi)
  
  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 `-strand`option.
  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](http://www.ebi.ac.uk/Tools/msa/)
  
  > **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**: @username_0 take one MLST gene and use that as a query to find the rest.  
  > **username_2**: @username_1 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**: @username_2 must admit I never know. Does this work in the case of contig assemblies?  
  > **username_1**: @username_0 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 username_2 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](https://sanger-pathogens.github.io/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**: @username_2 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.  
  
  ---
  
  