3
$\begingroup$

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.

Improve this question
$\endgroup$
  • $\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 Nov 4 '14 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$ – user2030378 Nov 4 '14 at 18:34
2
$\begingroup$
  • 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
Improve this answer
$\endgroup$
  • $\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$ – user2030378 Nov 4 '14 at 18:32
  • $\begingroup$ @user2030378 Do you want to look for redundant sequences ? $\endgroup$ – WYSIWYG Nov 5 '14 at 3:59
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.