How do I retrieve a specific region of a chromosome from a local blast database
1
0
Entering edit mode
3.7 years ago
b10hazard • 0

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
0
Entering edit mode
3.7 years ago
st.ph.n ★ 2.6k

You can do this using samtools.

samtools faidx /opt/human_genome_database/genome_GRCh38/GRCh38.primary_assembly.genome.fa 1:24957500-24958000
ADD COMMENT
0
Entering edit mode

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
0
Entering edit mode

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
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Please accept the answer if it satisfies the OP.

ADD REPLY

Login before adding your answer.

Traffic: 1667 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6