I generated a series of local BLAST databases using makeblastdb of metagenomic data and am searching for the presence of a particular gene. While I can do the normal BLAST analysis looking at e-values, % identity, etc. I haven't been able to figure out how to determine how many times a given sequence exists in the database.

So: given a BLAST database of metagenomic data where a particular gene may exist multiple times, how do I determine the copy number of a given BLAST hit?

Thanks in advance.

  • $\begingroup$ That depends on how you define a hit. Are we talking about homologous genes from different species? Multiple, slightly different copies of the same gene in a single genome? Multiple, identical copies of a single gene which is repeated in the database? $\endgroup$
    – terdon
    Commented Nov 4, 2014 at 15:01
  • $\begingroup$ I'm defining a hit based on the BLAST results -- I blasted seq1 against the database and it returned seq2...seqn. Depending on what cutoffs I use for bit scores / identity scores, etc. I might have three "hits" - seq2, seq3, and seq4. I want to know how many times seq2 showed up in the metagenomic data. $\endgroup$ Commented Nov 4, 2014 at 18:34

2 Answers 2

  • Define a "Hit" (based on some cutoff- evalue, score etc)
  • Get output in the tabular format
  • Count number of hits per query — it is usually given in the header; if you want to look for some selected hits (based on some cutoff, then you can parse the file and find out)

Example file (header):

# BLASTN 2.2.27+
# Query: TCONS_00036712 gene=XLOC_017996
# Database: ../nt_db/nt
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 1014 hits found

If you want to parse the header to look for number of hits then you can run this command in the terminal (if you have awk — there by default in linux and other unix based systems. For windows install gnuwin32).

awk -v OFS="\t" '/^# Query/{q=$3 FS $4} /^# .* hits found/{h[q]=$2} END{for(i in h){print i,h[i]}}' blastoutput.txt

To get the number of hits for a given criterion (for e.g. I just defined bitscore [12th column] to be >400)

awk -F "\t" '!/^#/ && $12>400{a[$1]++} END{for(i in a){print i,a[i]}}' blastoutput.txt
  • $\begingroup$ Great, but I have that already -- I know how many "hits" i.e. unique sequences there are. What I want to know is how many times that given hit's sequence is present in my assembled metagenomic dataset. $\endgroup$ Commented Nov 4, 2014 at 18:32
  • $\begingroup$ @user2030378 Do you want to look for redundant sequences ? $\endgroup$
    Commented Nov 5, 2014 at 3:59

BLAST is not the correct program to use in metagenomics analyses for determining the number of reads that map to a certain region of the genome, as it is optimised for single hits instead of a large number of different hits.

The sheer number of hits within a metagenomics database (which contains a large amount of redundant data) is likely to cause inefficient scaling and high memory usage when the detection of all possible hits is required.

If the number of hits required is low, it is possible to have it return those hits, but if you are looking at >1000 hits, it is probably recommended to use a mapper program such as BWA or Bowtie. They are designed for the express purpose of mapping reads (such as from an RNASeq experiment) onto genomes.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .