3
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Over at 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. $\endgroup$
    – Titus Brown
    Jun 13, 2017 at 18:35
  • $\begingroup$ 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. $\endgroup$
    – Daniel Harris
    Jun 13, 2017 at 19:32
  • $\begingroup$ Link to blast results (searched 3135 unique k-mers on one genome): ix.io/xqN $\endgroup$
    – Daniel Harris
    Jun 13, 2017 at 19:35
  • $\begingroup$ 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. $\endgroup$
    – Titus Brown
    Jun 13, 2017 at 20:52
  • $\begingroup$ I see! Here is the genome file I am using. ftp.patricbrc.org/patric2/current_release/AMR_genome_sets/… $\endgroup$
    – Daniel Harris
    Jun 14, 2017 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.