<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>

CC BY-SA 3.0
2017-06-13T18:35:10.313Z
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.
Titus Brown
CC BY-SA 3.0
2017-06-13T19:32:25.940Z
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.
Daniel Harris
CC BY-SA 3.0
2017-06-13T19:35:45.283Z
Link to blast results (searched 3135 unique k-mers on one genome): http://ix.io/xqN
Daniel Harris
CC BY-SA 3.0
2017-06-13T20:52:39.003Z
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.
Titus Brown
CC BY-SA 3.0
2017-06-14T14:09:26.443Z
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
Daniel Harris
CC BY-SA 3.0
2017-06-14T14:10:20.777Z
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.
Daniel Harris
CC BY-SA 3.0
2017-06-13T18:22:13.263Z
2017-06-13T18:22:13.263Z
Titus Brown
<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>

<sequence-annotation><genome><k-mer><python><blast>
Why do BLASTn and prokka not seem to be searching the whole fasta file?
# 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 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.  
> **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.  

---

