nuevo modo paralelo en BLAST+ 2.12.0

Hola,

BLAST+ es una de las pocas herramientas esenciales para cualquier persona que haga biología computacional, todos la hemos usado alguna vez. Simplificando, creo que el secreto es que hace una sola cosa (alinear localmente secuencias) y la hace bien. Eso, y que BLAST evoluciona, y de vez en cuando introduce novedades interesantes, que podemos repasar en las notas de publicación periódicas disponibles en https://www.ncbi.nlm.nih.gov/books/NBK131777

La versión 2.12.0 (28 de junio de 2021) incluye un nuevo modo de paralelización por lotes de secuencias input que beneficia entre otros a los algoritmos clásicos BLASTN, BLASTP y BLASTX. Los detalles puedes verlos en  https://www.ncbi.nlm.nih.gov/books/NBK571452 , aquí comparto un ejemplo de invocación, que requiere el argumento -mt_mode 1:

$ blastp –db swissprot -query BIGFASTA.fsa –out test.out -num_threads 32 -mt_mode 1

La siguiente figura muestra la ganancia en tiempo de cálculo en función del número de cores y hebras empleado:

Espero que os sea útil, hasta pronto,

Bruno

minimap2 vs BLASTN

Hola,
recientemente Heng Li publicó un trabajo (https://doi.org/10.1093/bioinformatics/bty191) describiendo un nuevo alineador genérico de nucleótidos que se llama minimap2, que podéis descargar en https://github.com/lh3/minimap2.

Figura tomada de https://doi.org/10.1093/bioinformatics/bty191

En el artículo se compara minimap2 en diferentes escenarios contra otros softwares alternativos, incluyendo su antecesor BWA mem y se destaca su velocidad y su versatilidad, ya que es capaz de alinear lecturas cortas, secuencias largas e incluso también puede alinear saltando intrones.

Yo lo que he hecho ha sido una prueba rápida para compararlo con BLASTN en el escenario habitual de GET_HOMOLOGUES-EST, donde se comparan por ejemplo todos los genes de una planta (Brachypodium distachyon) contra todos los genes de otra especie cercana (Oryza sativa). Esto es lo que he hecho:

# how many sequences
$ grep -c "^>" *fna
Bdistachyon.fna:36647
Osativa.fna:42189

# index and BLASTN search
$ ncbi-blast-2.6.0+/bin/makeblastdb -in Osativa.fna -dbtype nucl
$ ncbi-blast-2.6.0+/bin/blastn -query Bdistachyon.fna -db Osativa.fna \
   -out Bdistachyon.Osativa.blastn.tsv -dbsize 100000000 -evalue 1e-5 -outfmt 6
real	0m40.937s
user	0m40.280s
sys	0m0.636s

# index [assuming up 80% sequence identity] and minimap2 search
$ minimap2/minimap2 -x asm20 -d Oryza.mmi Osativa.fna
$ time minimap2/minimap2 Oryza.mmi Bdistachyon.fna > Bdistachyon.Osativa.minimap.paf
real	0m2.084s
user	0m3.360s
sys	0m0.300s

Ahora echemos un ojo a los alineamientos resultantes. Selecciono un par de secuencias, primero de BLASTN:

BdiBd21-3.2G0760100.1   LOC_Os01g70090.1        87.839  847     95      5       31      876     37      876     0.0     987
BdiBd21-3.2G0521100.1   LOC_Os01g37510.1        85.652  683     92      3       91      773     103     779     0.0     713

y ahora de minimap2, en formato PAF:

BdiBd21-3.2G0760100.1	876	155	776	+	LOC_Os01g70090.1	876	161	776	181	621	60	tp:A:P	cm:i:16	s
1:i:179	s2:i:0	dv:f:0.0980
BdiBd21-3.2G0521100.1	777	110	653	+	LOC_Os01g37510.1	783	122	659	87	543	60	tp:A:P	cm:i:10	s
1:i:85	s2:i:0	dv:f:0.1196

Al maneo para estos dos ejemplos podemos observar que:
i) el mejor hit de BLASTN y minimap coinciden
ii) los alineamiento de BLASTN son más largos


Hasta pronto, buenas vacaciones,
Bruno

HS-BLASTN as replacement of MEGABLAST

Hi,
this is hopefully the last of a series of posts where I explored software choices that might replace BLAST+ programs in some scenarios. Today I'll write about HS-BLASTN, a parallel nucleotide local aligner reported to accelerate the default BLASTN algorithm (megablast), while producing the same results:

HS-BLASTN algorithm overview, taken from http://nar.oxfordjournals.org/content/early/2015/08/06/nar.gkv784.long

Note that megablast is a fast choice for intra-species comparisons and typically retrieves sequence matches with nucleotide identities greater than 70% (see Figure 3 in http://eead-csic-compbio.github.io/get_homologues/manual-est).

We will benchmark HS-BLASTN (version hs-blastn-0.0.5+) using the same Hordeum vulgare and Arabidopsis thaliana sequences used in a previous post.

$ hs-blastn index bur-0.fasta
#[IndexBuilder] Time elapsed: 89.002952 secs.
#-rw-rw-r-- 1 contrera contrera 442M Jan  5 09:59 bur-0.fasta.bwt
#-rw-rw-r-- 1 contrera contrera  13M Jan  5 09:58 bur-0.fasta.header
#-rw-rw-r-- 1 contrera contrera 116M Jan  5 09:59 bur-0.fasta.sa
#-rw-rw-r-- 1 contrera contrera  58M Jan  5 09:58 bur-0.fasta.sequence

$ ncbi-blast-2.2.30+/bin/makeblastdb -in bur-0.fasta -dbtype nucl
#Adding sequences from FASTA; added 67259 sequences in 5.51543 seconds.
#-rw-rw-r-- 1 contrera contrera  12M Jan  5 10:06 bur-0.fasta.nhr
#-rw-rw-r-- 1 contrera contrera 789K Jan  5 10:06 bur-0.fasta.nin
#-rw-rw-r-- 1 contrera contrera  15M Jan  5 10:06 bur-0.fasta.nsq

It can be seen that indexing a FASTA file with 67K sequences is about 16x slower with HS-BLASTN than with standard NCBI mableblastdb, and produces much larger index files (629Mb vs 28Mb). Now let's review search performance:
 

hs-blastn align -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.hsblastn
#[HS-BLASTN] done. Elpased time: 9.0055 secs.

hs-blastn align -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.hsblastn -num_threads 20
#[HS-BLASTN] done. Elpased time: 1.6599 secs.

time ncbi-blast-2.2.30+/bin/blastn -task megablast -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -soft_masking true -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.blastn
#real 0m33.943s

perl _split_blast.pl 20 2000 ncbi-blast-2.2.30+/bin/blastn \
  -task megablast -query SBCC073_fLF.fasta -db bur-0.fasta -evalue 0.00001 -soft_masking true \
  -outfmt 6 -max_target_seqs 5 -out SBCC073_fLF.bur-0.blastn 
# runtime:  7 wallclock secs ( 2.81 usr  0.50 sys + 60.47 cusr  7.28 csys = 71.06 CPU)

It can be seen that HS-BLASTN is ~4x faster using a single thread. About the same speedup is obtained when 20 threads are used and BLASTN is parallelized with help from an external script (_split_blast.pl). There are small nuisances, such as the fact that no soft-masking is available, or the ocasional non-stable output order of hits with same score, but it seems worth it nevertheless,
have a good weekend,
Bruno