Answers.Body: <p>With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of </p>
  
  <pre>
  grep ^'>' fasta_file | wc -l
  </pre>
  
  <p>? This would show how many contigs you have in the file.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T18:35:10.313Z
Answers.Text: Over at http://github.com/dib-lab/khmer we'd be happy to help you work up a script that does what you need. Just post an issue.
Answers.UserDisplayName: Titus Brown
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T19:32:25.940Z
Answers.Text: I searched 3135 different 28-mers and about half of them actually matched. What I'm concerned with is why they only happen in the first 100 kpb of the genome. Using grep ^'>' fasta_file | wc -l , I found 32 contigs in the original fasta file.
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T19:35:45.283Z
Answers.Text: Link to blast results (searched 3135 unique k-mers on one genome): http://ix.io/xqN
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T20:52:39.003Z
Answers.Text: Hi Daniel, so it looks like you have multiple contigs. 32 * 100kb is about 3.2 Mbp so many of those contigs could be less than 100kb. It looks like maybe all the accession numbers are the same, from the BLAST results. Weird. Not sure what's going on, would have to take a look at the genome file to be sure.
Answers.UserDisplayName: Titus Brown
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T14:09:26.443Z
Answers.Text: I see! Here is the genome file I am using. ftp://ftp.patricbrc.org/patric2/current_release/AMR_genome_sets/Staphylococcus/methicillin/Resistant/1413323.3.fna
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T14:10:20.777Z
Answers.Text: I have tested a bunch of these genomes from that ftp directory though with the same blast command and they all search 1 contig only but it varies which contig.
Answers.UserDisplayName: Daniel Harris
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T18:22:13.263Z
Answers.LastActivityDate: 2017-06-13T18:22:13.263Z
Answers.UserDisplayName: Titus Brown
Body: <p>When I use blastn and prokka (I will detail exactly how I did so below) on a 2.8 million bp fasta file I get output start/end numbers that do not seem to cover the entire genome. </p>
  
  <p>Starting with a .fna genome such as genome.fna I ...</p>
  
  <p>1 - blastn</p>
  
  <p>Searched the genome for 3135 different 28-mers using BLASTn.</p>
  
  <pre><code>makeblastdb -in genome.fna -dbtype nucl -parse_seqids -out ./output/genome
  </code></pre>
  
  <p><em>next command in python</em></p>
  
  <pre><code>blast_tsv = NcbiblastnCommandline(query=Q, db=DB, perc_identity=100, outfmt=6, out=(OUT))
      stdout, stderr = blast_tsv()
  </code></pre>
  
  <p><em>Q is this list of k-mers. DB is the database created</em></p>
  
  <p>This outputs a list of search results but the important thing is that none of the start/end (columns 7 and 8) integers were greater than 100,000 yet the entire genome is 2.8 million base pairs long. I can provide the link to this file in the comments. Does this have to do with blastn stopping searching after finding one match? And if so how can I tell blastn to search for every match for each k-mer? (I'm open to using other blastn programs other than biopython)</p>
  
  <p>2 - prokka</p>
  
  <p>I used prokka to create gff files for the genome. </p>
  
  <pre><code>prokka-1.12/bin/prokka —setupdb
  prokka-1.12/bin/prokka -kingdom Bacteria -rfam -outdir ./prokka_database/genome -force -prefix "${genome/.fna/}" genome.fna
  </code></pre>
  
  <p>Th gff file produced by this command only described genes with start and end positions less than 200k. I can provide the link to this file as well. </p>
  
Tags: <sequence-annotation><genome><k-mer><python><blast>
Title: Why do BLASTn and prokka not seem to be searching the whole fasta file?
ThreadText: # Why do BLASTn and prokka not seem to be searching the whole fasta file?
  
  (Asked by: username_0 on 2017-06-13)
  
  When I use blastn and prokka (I will detail exactly how I did so below) on a 2.8 million bp fasta file I get output start/end numbers that do not seem to cover the entire genome.
  
  Starting with a .fna genome such as genome.fna I ...
  
  1 - blastn
  
  Searched the genome for 3135 different 28-mers using BLASTn.
  
  ```
  makeblastdb -in genome.fna -dbtype nucl -parse_seqids -out ./output/genome
  
  ```
  
  *next command in python*
  
  ```
  blast_tsv = NcbiblastnCommandline(query=Q, db=DB, perc_identity=100, outfmt=6, out=(OUT))
      stdout, stderr = blast_tsv()
  
  ```
  
  *Q is this list of k-mers. DB is the database created*
  
  This outputs a list of search results but the important thing is that none of the start/end (columns 7 and 8) integers were greater than 100,000 yet the entire genome is 2.8 million base pairs long. I can provide the link to this file in the comments. Does this have to do with blastn stopping searching after finding one match? And if so how can I tell blastn to search for every match for each k-mer? (I'm open to using other blastn programs other than biopython)
  
  2 - prokka
  
  I used prokka to create gff files for the genome.
  
  ```
  prokka-1.12/bin/prokka —setupdb
  prokka-1.12/bin/prokka -kingdom Bacteria -rfam -outdir ./prokka_database/genome -force -prefix "${genome/.fna/}" genome.fna
  
  ```
  
  Th gff file produced by this command only described genes with start and end positions less than 200k. I can provide the link to this file as well.
  
  **Tags:** sequence-annotation, genome, k-mer, python, blast
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2017-06-13. Score: 6
  
  With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of
  
  ```
  
  grep ^'>' fasta_file | wc -l
  
  ```
  
  ? This would show how many contigs you have in the file.
  
  > **Comments:**
  > **username_1**: Over at http://github.com/dib-lab/khmer we'd be happy to help you work up a script that does what you need. Just post an issue.  
  > **username_0**: I searched 3135 different 28-mers and about half of them actually matched. What I'm concerned with is why they only happen in the first 100 kpb of the genome. Using grep ^'>' fasta_file | wc -l , I found 32 contigs in the original fasta file.  
  > **username_0**: Link to blast results (searched 3135 unique k-mers on one genome): http://ix.io/xqN  
  > **username_1**: Hi username_0, so it looks like you have multiple contigs. 32 * 100kb is about 3.2 Mbp so many of those contigs could be less than 100kb. It looks like maybe all the accession numbers are the same, from the BLAST results. Weird. Not sure what's going on, would have to take a look at the genome file to be sure.  
  > **username_0**: I see! Here is the genome file I am using. ftp://ftp.patricbrc.org/patric2/current_release/AMR_genome_sets/Staphylococcus/methicillin/Resistant/1413323.3.fna  
  > **username_0**: I have tested a bunch of these genomes from that ftp directory though with the same blast command and they all search 1 contig only but it varies which contig.  
  
  ---
  
  