Question: How do I retrieve a specific region of a chromosome from a local blast database
0
gravatar for b10hazard
23 months ago by
b10hazard0
b10hazard0 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 • 908 views
ADD COMMENTlink modified 23 months ago • written 23 months ago by b10hazard0
0
gravatar for st.ph.n
23 months ago by
st.ph.n2.5k
Philadelphia, PA
st.ph.n2.5k wrote:

You can do this using samtools.

samtools faidx /opt/human_genome_database/genome_GRCh38/GRCh38.primary_assembly.genome.fa 1:24957500-24958000
ADD COMMENTlink modified 23 months ago • written 23 months ago by st.ph.n2.5k

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 REPLYlink modified 23 months ago • written 23 months ago by b10hazard0

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 REPLYlink modified 23 months ago • written 23 months ago by st.ph.n2.5k

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

ADD REPLYlink written 23 months ago by b10hazard0

Please accept the answer if it satisfies the OP.

ADD REPLYlink written 23 months ago by st.ph.n2.5k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1043 users visited in the last hour