Sequence extract from multifasta using blastall coordinates

Question

[chland]

Hi everyone.

This question has been asked in various forms before, however am new to the field and trying to learn and the answers posted don't seem to be working. Sorry, I need to specify that I'm trying to extract the sequence corresponding to the Subject start and subject end coordinates for each hit from a multifasta file of concatenated sequences. Thanks for the useful command though, I will make note of that for other purposes.

I ran blastall search that returned a file of hits that look like this with up to 500 hits ( -b 500) to the database of concatenated fasta files. I'm looking for a method to extract the sequence using the s start and s end for each hit along with the corresponding subject ID from the database of concatenated sequences.
Thanking you in advance. Chandly

# BLASTN 2.2.26 [Sep-21-2011]
# Query: querysequence.fasta
# Database: blast_db.fna
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
NCTC11168   NZ_FPEE01000064.1   100.00  960 0   0   1   960 16113   15154   0.0 1861
NCTC11168   NZ_PSND01000007.1   100.00  960 0   0   1   960 49840   50799   0.0 1861
NCTC11168   NZ_PIBG01000017.1   100.00  960 0   0   1   960 12161   13120   0.0 1861
NCTC11168   NZ_PHZN01000005.1   100.00  960 0   0   1   960 49939   50898   0.0 1861
NCTC11168   NZ_NFPZ01000033.1   100.00  960 0   0   1   960 16585   15626   0.0 1861

[GenoMax]

Something like awk 'BEGIN{OFS="\t"}{print $2,$9,$10}' your_blast_file will get you the start and stops. Are the subject sequences in a local multi-fasta file? You can find many threads here to extract sequences from that point on.

[chland]

Thank you genomax for trying to help. I've been looking though answers posted and am not finding an answer that works and that I can understand. Yes, the sequences are in a local multi-fasta file created by running the cat command on sequence data ( complete, contigs/scaffold) from NCBI assembly.

[Joe]

BioPythons slice notation will be want you need here.

I have a gist which employs a blastfile to slice up genbanks, but you could do the same with a Rasta read in to Biopython as a SeqRecord:

[chland]

Thank you, I've downloaded the script from github, but by any chance would you be able to provide a usage statement example? Also, the assemblies were grabbed from refseq in case that matters. I have a database consisting of multifasta sequences called something like mydb.fasta and blastall hits as in the original post.

[Joe]

I’m not in front of a computer right now to put together the full solution.

The script will explain how it works if you run python script.py -h. But it wont work for fast as as it is. You’ll need to write something a little bespoke yourself.

The line you’ll be most interested in is 89, the slice function

[chland]

Thanks all, but I'm starting from the very beginning here and have not made it to any kind of scripting so it's all greek.

Answer 1

[h.mon]

My experience (but I never seriously benchmarked) is that, for big fasta files, method 2 bellow is faster than method 1.

Method 1:

a) index the fasta file containing the original sequences:

samtools faidx file.fasta

b) extract sequences interval of interest:

samtools faidx file.fasta $(cat blast.txt | grep -v "^#" | awk 'BEGIN{OFS="\t"}{print $2":"$9"-"$10}')

Method 2:

a) extract the subject hit accessions:

cat blast.txt | grep -v "^#" | cut -f2 | sort | uniq > accessions.txt

b) extract the subject hit sequences with formatdb (I will illustrate here with blastdbcmd, which is the equivalent command for the Blast+ suite):

blastdbcmd -db nt -entry_batch accessions.txt -out subjects.fasta

c) index the subjects fasta:

samtools faidx subjects.fasta

d) extract sequences interval of interest:

samtools faidx subjects.fasta $(cat blast.txt | grep -v "^#" | awk 'BEGIN{OFS="\t"}{print $2":"$9"-"$10}')

[chland]

Spoke too soon. Thanks @h.mon. I've used the blast suite of tools and will try this method. The python script is too advanced for a complete beginner

[baurumon]

Hi, Thank you for this useful method @h.mon

what if the $9>$10. As you can see above blast out. when i tried with my data i am getting

665:1927-1645 [faidx] Zero length sequence: 665:1927-1645

when there is $9>$10.