Differences in BLAST-ings

[kubu4]

I'll start out by qualifying this post as possibly being moot from the start, as the following is definitely not an "apples to apples" comparison. Okay, with that out of the way...

Recently, Steven and I took two different approaches to annotating (via BLASTing) the C.virginica genome and we got very different results in the number of BLAST matches.

@sr320's approach (Notebook)

1. Use translated CDS FastA as queries (GCF_002022765.2_C_virginica-3.0_translated_cds.faa.gz)
   • 60,213 sequences (grep -c "^>" GCF_002022765.2_C_virginica-3.0_translated_cds.faa)
   • 34,608 unique genes (grep -o '\[gene=[^]]*\]' | grep -oP '(?<=\[gene=)[^]]*' GCF_002022765.2_C_virginica-3.0_translated_cds.faa | sort -u | wc -l)
2. Use reviewed UniProt proteins FastA as BLAST database (uniprot_sprot.fasta.gz).
3. Use NCBI BLAST to make BLAST database.
4. Use NCBI BLASTp to BLAST.
5. Do some gene ontology stuff.

Results:

BLAST

• 46,130 matches (wc -l Cvir_transcds-uniprot_blastp.tab)
• 10,370 unique SwissProt IDs (awk -F "|" '{print $3}' Cvir_transcds-uniprot_blastp.tab | sort -u | wc -l)

Gene Ontology

• 17,971 genes annotated

---

Sam's approach

1. Extract gene sequences as FastA (Notebook).

   • 38,838 genes (wc -l C_virginica-3.0_Gnomon_genes.bed)

2. Use reviewed UniProt proteins FastA as BLAST database (uniprot_sprot.fasta.gz).

3. Use DIAMOND to make DIAMOND BLAST database.

4. Use DIAMOND BLASTx to BLAST. (Notebook)

Results:

BLAST

• 5,869 matches (wc -l GCF_002022765.2_C_virginica-3.0-genes.blastx.outfmt6)
• 3,113 unique SwissProt IDs (awk -F "|" '{print $2}' GCF_002022765.2_C_virginica-3.0-genes.blastx.outfmt6 | sort -u | wc -l)

---

Sooo, what do people think about this? Anyone have any thoughts/opinions? E-value cutoffs for BLASTing for both approaches were 1e^-20 and 1e^-25 for each approach, respectively.

[sr320]

One thing I would be interested in is head to head comparison of blast and diamond blast (same query/db/sertings)

[kubu4]

I ran an NCBI BLASTx, which is comparable to DIAMOND BLASTx (output directory with files).

I say "comparable" because I used the same e-value cutoff and the same UniProt protein database (i.e. constructed from the same starting FastA). However, DIAMOND BLASTx has some "sensitivity" settings that can be set and I'm not sure how those compare to NCBI BLASTx...

Anyway, here is a table of the two results:

	DIAMOND	NCBI
Matches	5,869	6,025
Unique SwissProt	3,113	3,029
Run time	00:01:30	17:07:02

The difference in run times is wild! Haha!

DIAMOND: 1 minute and 30secs
NCBI: 17hrs

NOTE: The NCBI BLASTx had multiple hits per query, despite using -num_target_seqs 1 set... The number of matches reported above are unique matches.

[sr320]

Very interesting! Note I think '-max_hsps 1' added to max target will get you uniques.

[kubu4]

Note I think '-max_hsps 1' added to max target will get you uniques

Thanks. I haven't used NCBI BLAST in so long that I couldn't remember the magic combination to get unique matches.

[kubu4]

Personally, I'm happy to see similar results between the two when using BLASTx, but it doesn't clear up why there's such a large difference between the BLASTp and BLASTx results - which is what prompted this discussion. I still find that difference surprising.