How do I retrieve a specific region of a chromosome from a local blast database

Question: How do I retrieve a specific region of a chromosome from a local blast database

3.4 years ago by

b10hazard • 0 wrote:

I created a local blast database of the human genome (version GRCh38) using a fasta file I downloaded from NCBI. The fasta file records are named after the chromosomes so chromosome 1 has the record "chr1 1". I would like to query this database to get base pairs 24957500-24958000 in a fasta file format. How do I do this? I did some reading and found that there is a command line tool called blastdbcmd. I tried this....

blastdbcmd -db /opt/human_genome_database/genome_GRCh38/GRCh38.primary_assembly.genome.fa -entry "chr1 1" -range 24957500,24958000 -out output.fasta

However, this gives me the error...

Error: [blastdbcmd] Entry not found: chr1 1
Error: [blastdbcmd] Entry or entries not found in BLAST database

What is wrong with my command? Or am I using the wrong application for this? Thanks!

blast • 1.4k views

ADD COMMENT • link •

modified 3.4 years ago • written 3.4 years ago by b10hazard • 0

3.4 years ago by

st.ph.n • 2.6k

Philadelphia, PA

st.ph.n • 2.6k wrote:

You can do this using samtools.

samtools faidx /opt/human_genome_database/genome_GRCh38/GRCh38.primary_assembly.genome.fa 1:24957500-24958000

ADD COMMENT • link modified 3.4 years ago • written 3.4 years ago by st.ph.n • 2.6k

I install samtools using apt-get (I am on ubuntu). I ran the command you posted but it did not produce a fasta file. All I got was this output on the command line....

>1:24957500-24958000

I grabbed the whole assembly because I need to use it for other projects that involve blast searches.

ADD REPLY • link modified 3.4 years ago • written 3.4 years ago by b10hazard • 0

Double check what is right after the '>' in the fasta, which goes before the ':' in the command. I tested with, chr 1, and this is what the header looks like, and the output.

>1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF

-bash-4.1$ samtools faidx Homo_sapiens.GRCh38.dna.chromosome.1.fa 1:1-10

>1:1-10
NNNNNNNNNN

ADD REPLY • link modified 3.4 years ago • written 3.4 years ago by st.ph.n • 2.6k

Ah, you're correct, that was it. The fasta has chr1 as the record ID. It works now. Thanks!

ADD REPLY • link written 3.4 years ago by b10hazard • 0

Please accept the answer if it satisfies the OP.

ADD REPLY • link written 3.4 years ago by st.ph.n • 2.6k

Similar posts • Search »