DIAMOND as alternative to BLASTP

Hi,
in a recent post I discussed the performance of DIAMOND when aligning nucleotide sequences against a peptide database. However, as this tool can also do protein similarity searches, I wanted to test that too. Here I summarize those results.

We used all protein sequences encoded in ecotype Bur-0 of Arabidopsis thaliana (n=53,027) and the proteome of Mycobacterium tuberculosis strain H37Rv (n=3,906). We scanned these against Swissprot, downloaded from ftp.uniprot.org and formatted as follows:

$ makeblastdb -in uniprot_sprot.fasta -dbtype prot                   
$ diamond makedb --in uniprot_sprot.fasta --db uniprot_sprot.fasta  

As earlier, we used v0.8.25, obtained from https://github.com/bbuchfink/diamond.

The actual sequence similarity searches were performed as follows:

$ diamond blastp -p 30 -d swissprot -q bur-0.faa -o bur-0.diamond.tsv \
  --max-target-seqs 1000 --evalue 0.00001 \
  --outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend \
  sstart send evalue bitscore qcovhsp sseq 

$ diamond blastp -p 30 -d swissprot -q bur-0.faa -o bur-0.diamond.sensitive.tsv \
  --sensitive ...

$ diamond blastp -p 30 -d swissprot -q bur-0.faa -o bur-0.diamond.more.tsv \
  --more-sensitive ...

$ time ncbi-blast-2.2.30+/bin/blastp -num_threads 30 -db swissprot \
  -query bur-0.faa -out bur-0.blastp.tsv 
  -max_target_seqs 1000 -evalue 0.00001 -outfmt '6 std qcovhsp sseq' 

# calling _split_blast.pl as used in get_homologues

$ get_homologues-est/_split_blast.pl 30 250 ncbi-blast-2.2.30+/bin/blastp \
  -db swissprot -query bur-0.faa -out bur-0.split.blastp.tsv \
  -max_target_seqs 1000 -evalue 0.00001 -outfmt \'6 std qcovhsp sseq\' 

Those commands were then run replacing bur-0.faa with Mtuberculosis_H37Rv.faa and using 8 CPU cores instead of 30. The obtained alignments were compared with a script (available upon request). Wall-clock time and matched hits are computed in relation to BLASTP (control). Alignment comparisons were performed with all BLASTP hits (red) and the subset of hits with identity >= 50% (blue).

The take home message is that DIAMOND "more sensitive" is 20x to 100x faster than BLASTP in these tests, with roughly 15% less sensitive overall, which is reduced to 5-9% when more remote homologues matter.

     search_strategy    time(s) matched =length longer shorter matched =length longer shorter

control (NCBI BLASTP)     16440  [         all hits          ]  [       IDENTITY:50%        ]
bur-0.diamond                55   0.402   0.804  0.109   0.088   0.848   0.916  0.063   0.021
bur-0.diamond.sensitive     192   0.774   0.794  0.116   0.090   0.906   0.916  0.062   0.022
bur-0.diamond.more          263   0.832   0.797  0.115   0.088   0.910   0.916  0.062   0.022
bur-0.split.blastp         6915   1.000    
                                          
control (NCBI BLASTP)      2657  [         all hits          ]  [       IDENTITY:50%        ]
Mtub.diamond                 34   0.483   0.853  0.075   0.072   0.930   0.932  0.051   0.017
Mtub.diamond.sensitive      124   0.832   0.823  0.091   0.087   0.957   0.931  0.052   0.017
Mtub.diamond.more           141   0.853   0.821  0.092   0.086   0.958   0.930  0.052   0.017
Mtub.split.blastp          2334   1.000 

Bruno

PS I post another test done by David A Wilkinson, from Massey University, New Zealand:

"I ran get homologues in conjunction with blastp, diamond in "standard mode" and diamond in "more-precise mode" for a small dataset.The dataset includes 20 genomes from the genus Leptospira - importantly, they are not all from the same "species", and for some pairwise comparisons will have relatively distant genomes and gene identities...

I used the ORTHOMCL algorithm in get-homologues for clustering. My numbers are in total agreement with published literature.Here are the comparative results:

	#clusters	#core	#softcore	#shell	#cloud
BLASTP	10398	1523	2197	2269	5932
DIAMOND_STANDARD	10829	1421	2146	2328	6355
DIAMOND_MORE_PRECISE	10410	1535	2215	2269	5926
	%clusters	%core	%softcore	%shell	%cloud
BLASTP	100%	100%	100%	100%	100%
DIAMOND_STANDARD	104%	93%	98%	103%	107%
DIAMOND_MORE_PRECISE	100%	101%	101%	100%	100%

"

Chuleta de Python 3 para principiantes

Hoy quiero compartir con vosotros una chuleta de Python 3 que he creado para mis estudiantes, he intentado recopilar en dos caras de DIN-A4 los tipos de datos, operadores, métodos, funciones y otros contenidos útiles que uno necesita tener a mano cuando comienza a programar en Python 3 (incluso meses más tarde).

Puedes descargarte la chuleta en formato PDF de mi blog en inglés o usarla online:

DIAMOND as replacement of BLASTX

Hi,
about a year ago my colleague Javier Tamames told me about a piece of software called DIAMOND which could be thought of as a replacement of some BLAST tools, particularly BLASTP and BLASTX.

Author's benchmark of DIAMOND, taken from http://www.nature.com/nmeth/journal/v12/n1/full/nmeth.3176.html .

Recently we tested in-house the BLASTX alternative and I summarize here the results of finding protein-coding segments in transcripts from Arabidopsis thaliana (n=67,259) and barley (Hordeum vulgare, n=76,362) by comparing them to Swissprot, downloaded from ftp.uniprot.org.

You can get or build DIAMOND from https://github.com/bbuchfink/diamond.

First, we had to format Swissprot to support searches:

$ makeblastdb -in uniprot_sprot.fasta -dbtype prot                   # produces 3 files
$ diamond makedb --in uniprot_sprot.fasta --db uniprot_sprot.fasta   # produces 1 file

Then we could run the actual nucleotide-to-peptide sequence searches allocating 30 CPU cores. Note that files bur-0.fasta and SBCC073_fLF.fasta contain the A.thaliana and barley transcripts:

$ diamond blastx -p 30 -d swissprot -q bur-0.fasta -o bur-0.diamond.tsv \
  --max-target-seqs 1 --evalue 0.00001 \
  --outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovhsp sseq 
#Total time = 36.2444s

$ diamond blastx -p 30 -d swissprot -q bur-0.fasta -o bur-0.diamond.sensitive.tsv \
  --sensitive ...

#Total time = 144.636s

$ diamond blastx -p 30 -d swissprot -q bur-0.fasta -o bur-0.diamond.more.tsv \
  --more-sensitive ...

#Total time = 194.35s

$time ncbi-blast-2.2.30+/bin/blastx -num_threads 30 -db ~/db/swissprot \

  -query bur-0.fasta -out bur-0.blastx.tsv \
  -max_target_seqs 1 -evalue 0.00001 -outfmt '6 std qcovhsp sseq' 
#real    351m30.313s
#user    8294m28.727s
#sys    16m26.575s

Similar searches were performed with barley sequences and the results then compared with help from a Perl script which required aligments to be 5% longer/shorter to be considered significantly loger or shorter, respectively. The total BLASTX alignments were 44,970 for A.thaliana and 32,887 for barley:

diamond_strategy                 matched same_hit  same_length longer  shorter
bur-0.diamond.tsv                  0.935    0.895        0.921  0.031    0.047
bur-0.diamond.sensitive.tsv        0.973    0.902        0.919  0.037    0.044
bur-0.diamond.more.tsv             0.973    0.902        0.918  0.037    0.045

SBCC073_fLF.diamond.tsv            0.889    0.807        0.852  0.071    0.077
SBCC073_fLF.diamond.sensitive.tsv  0.960    0.831        0.856  0.076    0.067
SBCC073_fLF.diamond.more.tsv       0.961    0.831        0.856  0.076    0.067 

All in all, our tests suggest a reduction in computing time of 2 orders of magnitude with a 4% loss of sensitivity and alignments to protein sequences of the same length in most cases.

See you,
Bruno

Correlaciones y regresiones en R

Hola,
al analizar datos es frecuente que recurrir al cálculo de correlaciones, que permiten establecer si dos variables covarían, y de regresiones, que permiten modelar de qué manera conociendo una variable podemos estimar la otra. De hecho ya habíamos tocado este tema en este blog.

Hoy quisiera compartir el estupendo material sobre estas materias que ha preparado nuestro colega Pablo Vinuesa, que podéis explorar en línea y con ayuda de Rstudio, por ejemplo, para probar todo el código de los ejemplos:

• correlación: [http://www.ccg.unam.mx/~vinuesa/R4biosciences/docs/Tema8_correlacion_presentacionR.html#/]


• regresión: [http://www.ccg.unam.mx/~vinuesa/R4biosciences/docs/Tema9_regresion_presentacionR.html#/]

Un saludo,
Bruno

extraer TSV de todas las hojas de un libro excel

Hola,
cuando intercambiamos datos en al laboratorio a menudo usamos libros Excel, en formato XLSX, con varias hojas. Ahora bien, si luego queremos convertirlos a ficheros con valores separados por tabuladores (TSV) o comas (CSV), MS Excel sólo te permite hacerlo de hoja en hoja, lo cual es un poco latoso para libros gordos. Aquí os pongo un script en Perl para hacer esta tarea desde el terminal, que requiere instalar el módulo Spreadsheet::ParseXLSX, por ejemplo con:

$ sudo cpan -i Spreadsheet::ParseXLSX

El código es:

#!/usr/bin/perl -w
# get separate tab-separated (TSV) files from sheets in Excel .xlsx file 
use strict;
use Spreadsheet::ParseXLSX;

die "# usage: $0 <infile.xlsx> \n" if(!$ARGV[0]);

my $parser = Spreadsheet::ParseXLSX->new();
my $book = $parser->parse($ARGV[0]);
foreach my $sheet (@{$book->{Worksheet}})
{
  open(TSV,'>',"$sheet->{'Name'}.tsv") || 
    die "# cannot create $sheet->{'Name'}.tsv: $!\n";

  foreach my $row ($sheet->{'MinRow'} .. $sheet->{'MaxRow'})
  {
    foreach my $col ($sheet->{'MinCol'} ..  $sheet->{'MaxCol'})
    {
      print TSV "$sheet->{'Cells'}->[$row]->[$col]->{'Val'}\t";
    } print TSV "\n";
  }

  close(TSV);
}

Si guarmos el código en un fichero de nombre xlsx2multitab.pl, podemos invocarlo de la siguiente manera:

$ perl xlsx2multitab.pl libro.xlsx

En el directorio actual se guardará un fichero .tsv por cada hoja del libro,
hasta luego,
Bruno

Seminario de introducción a la metagenómica

Os quiero anunciar que el 19 de Octubre impartiré el workshop titulado "Introduction to Bioinformatics applied to Metagenomics and Community Ecology" como parte de la conferencia ​Community Ecology for the 21st Century (Évora, Portugal).

Si estáis interesados, podéis contactar a los organizadores de la conferencia en el siguiente enlace, todavía hay plazas disponibles en el workshop.

Durante el curso presentaré la nueva herramienta AmpliTAXO para el análisis sencillo y online de datos de NGS de RNA ribosomal y otros marcadores.

El curso consistirá en 2 partes, la primera teórica donde se expondrán los retos de la metagenómica, las posibilidades de las nuevas técnicas de secuenciación y el funcionamiento de las herramientas de análisis más habituales (UPARSE, QIIME, MOTHUR). La segunda parte será práctica y consistirá en el análisis de datos metagenómicos reales obtenidos por NGS.

Podéis encontrar más información en inglés en mi nuevo blog y próximamamente en la página de la conferencia (pendiente de actualizar).

Os dejo un pequeño adelanto de los contenidos...

PHMMER como alternativa a BLASTP

Hola,
hoy quería hablar de herramientas de búsqueda de secuencias de proteína, uno de los caballos de batalla en nuestro trabajo. Es un tema ya tocado en este blog, por ejemplo cuando hablamos de deltablast y cs-blast, pero que sigue siendo de actualidad.

En esta ocasión ha sido un artículo de Saripella et al el que me lo ha vuelto a recordar, comparando algoritmos que usan perfiles de secuencia (CS-BLAST, HHSEARCH and PHMMER) con algoritmos convencionales (BLASTP, USEARCH, UBLAST and FASTA). Es uno de esos artículos aburridos pero necesarios, donde se compara por vía independiente la validez de los mejores algoritmos para esta tarea, y se guardan los tiempos de cálculo empleados por cada uno de ellos (usando un core de CPU), para anotar dominios conocidos de secuencias de SwissProt extraídos de 3 repositorios complementarios (Pfam, Superfamily y CATH):

Figura original de http://bioinformatics.oxfordjournals.org/content/32/17/2636.full.

En el artículo se calculan áreas bajo curvas ROC para caracterizar a los diferentes algoritmos. De todos ellos destacaría dos:

1) BLASTP por ser muy rápido y muy preciso, con áreas de 0.908, 0.857 y 0.878 para las 3 colecciones de dominios y el tiempo que se muestra en la gráfica para buscar 100 secuencias al azar.

2) PHMMER por ser marginalmente más lento que BLASTP pero con una ganancia en área de aproximadamente el 3% (0.922,0.903,0.903).

Además, PHMMER es muy fácil de usar. Lo descargas de http://eddylab.org/software/hmmer3/CURRENT/ y solamente necesitas un archivo FASTA con la(s) secuencia(s) problema y otro con un repositorio de secuencias como SwissProt:

$ hmmer-3.1b1-linux-intel-x86_64/src/phmmer problema.faa uniprot_sprot.fasta

Hasta pronto,
Bruno

One-liner para calcular el N50 de un ensamblaje

Hola,
el estadístico N50 se usa a menudo para resumir a grosso modo un ensamblaje de secuencias, que generalmente es un fichero FASTA con una serie de contigs.
Modificando la definición de la wikipedia, N50 es la longitud de contigs tal que usando contigs de igual o mayor tamaño recuperamos la mitad de las bases de ese ensamblaje.

Figura tomada de http://bmpvieira.github.io/assembly14.

Hoy solamente quería compartir un comando one-liner de Perl que nos permite calcularlo:

$ perl -lne 'if(/^>(\S+)/){$h=$1} else {$l=length($_);$TL+=$l;$L{$h}+=$l} \

END{ foreach $s (sort {$L{$b}<=>$L{$a}} keys(%L)){ $t+=$L{$s}; \

if($t>$TL/2){ print "N50=$L{$s}";exit }}}' assembly.fasta

Cuidado al copiarlo al terminal desde el blog, mejor hacerlo línea a línea,
hasta pronto, Bruno

Las proteínas dúctiles

Hola,
ayer me regaló mi colega Inma Yruela una copia de su reciente libro "Las proteínas dúctiles". Es éste un libro de divulgación científica donde se explica de manera amena cómo las proteínas desordenadas, que Inma y yo estudiamos en las plantas (ver artículos 1 y 2), se salen del paradigma dominante "secuencia -> estructura 3D -> función", dado que desempeñan funciones importantes en los seres vivos, al parecer más en los eucariotas, sin tener una estructura definida en al menos una parte significativa de su secuencia.

Enlaces en CSIC y Amazon.

Desorden intrínseco de los extremos N y C terminal (en rojo) de la proteína supresora tumoral p53. Figura tomada de doi: 10.1146/annurev-biochem-072711-164947.

Si bien en la literatura está ya aceptado el término de desorden intrínseco (IDP en inglés) para definir a partes de estas proteínas, Inma propone reemplazarlo por proteínas dúctiles, que muestra de manera explícita que son regiones altamente moldeables en vez de destacar que carecen de órden. Os invito a que lo leáis porque es entretenido, por su exposición sencilla de las metodologías bioquímicas, biofísicas y bioinformáticas que se usan para estudiarlas, y por la descripción de los procesos biológicos en los que participan, por ejemplo su relación con el splicing. Lo más fascinante de estas proteínas es que seguramente nos falta mucho por saber de ellas,
un saludo,
Bruno

Cómo elegir software para simular reads

Hola,
si trabajas con cierta frecuencia con datos de ultrasecuenciación, es decir, con ficheros en formato FASTQ, a menudo te habrás encontrado con la situación de que te gustaría tener más datos para validar un programa publicado, o para diseñar tu propia tubería de análisis. Pues bien, una posibilidad real desde hace varios años es simular tus propios reads o lecturas a partir de secuencias genómicas del organismo en cuestión y parámetros como la tasa de error o la plataforma de secuenciación empleada. Como ocurre frecuentemente en Biología Computacional, hay una gran variedad de software publicado para esta tarea y elegir no es trivial. El diagrama a continuación es un árbol de decisión que os guiará para esta tarea:

Figura prestada de Escalona, Rocha & Posada (2016) doi:10.1038/nrg.2016.57 http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html

Buen finde,
Bruno

email desde servidor TLS con perl

Hola,
en la entrada de hoy muestro una manera de enviar correo electrónico desde un servidor seguro (STARTTLS) por medio de un script Perl. Hasta hace poco lo hacía con el vetusto módulo Net::SMTP::TLS, pero esta solución era muy frágil y se rompía cada vez que otros módulos (fundamentalmente IO::Socket::SSL) no iban a la par en sus versiones. De hecho, en un sistema Ubuntu 14 no he logrado hacerla funcionar actualizando esos módulos, pero cambio he descubierto que con el módulo Net::SMTPS es más sencillo. Aquí queda el código:

 use strict;  
 use MIME::Lite;  
 use Net::SMTPS;  
   
 my $HELLOIP   = 'xxx.xxx.xxx.xxx'; # your IP, must be allowed by $SECURESERVER  
 my $SECURESERVER = 'smtp.domain.com'; # can also be an IP address  
 my $SECUREPORT  = 587;  
 my $SECUREUSER  = 'secure_user';   # user login in secure server  
 my $SECUREPASSWD = 'zzzzzzzzzzz';   # user password, better read from safe file  
 my $SECURESENDER = 'user@domain.com'; # must be owned by $SECUREUSER  
     
 send_email(   
  $HELLOIP,$SECURESERVER,$SECUREPORT,  
  $SECUREUSER,$SECUREPASSWD,  
  $SECURESENDER,'recipient@test.domain.org',   
  'probando','texto de prueba'  
  );  
   
 sub send_email  
 {  
  my ($hello,$host,$port,$user,$pass,$sender,$recipient,$subject,$text) = @_;  
   
  # connect to server  
  my $smtp = new Net::SMTPS(  
   $host,  
   Hello  =>$hello,  
   Port  =>$port,  
   doSSL  =>'starttls',  
   Debug  =>0 # change to 1 to see behind the curtains  
  );  
   
  $smtp->auth($user,$pass); # login  
   
  # Create the multipart container  
  my $msg = MIME::Lite->new (  
   From  => $sender,  
   To   => $recipient,  
   Subject => $subject,  
   Type  =>'multipart/mixed'  
  ) or die "# send_email : error creating multipart container: $!\n";  
   
  # add main-text if required  
  if($text)  
  {  
   $msg->attach (  
    Type => 'text/plain',  
    Data => $text  
   ) or die "# send_email : error adding the text message part: $!\n";  
  }   
   
  # actually send email  
  $smtp->mail($sender);  
     
  if($smtp->to($recipient))  
  {   
   $smtp -> data();  
   $smtp -> datasend( $msg->as_string() );   
   $smtp -> dataend();  
  }  
  else{ print "# ERROR: ". $smtp->message() }   
   
    $smtp->quit();  
 }  
   
   

Hasta la próxima,
Bruno

pseudoalineamientos y conteos de tránscritos

Hola,
durante el análisis de unos experimentos de RNAseq mi colega Carlos Cantalapiedra y yo nos desesperábamos al calcular los valores de expresión de cada transcrito en una serie de condiciones, dado que para cada una de ellas era necesario alinear las lecturas/reads originales (en formato FASTQ) contra los transcritos ensamblados. Al menos éste era el protocolo habitual antes de aparecer el algoritmo genérico de pseudoalineamiento, que se describe en la siguiente figura:

Pseudoalineamiento de reads a partir de su composición en k-meros (k=3 en el ejemplo) y un  vector de sufijos de un transcriptoma. Figura tomada de http://bioinformatics.oxfordjournals.org/content/32/12/i192.full.  

Este tipo de algoritmos, implementados en software como kalllisto o RapMap (éste último con licencia GPL), permiten estimar de manera muy eficiente y precisa con qué transcritos son compatibles las lecturas de un archivo FASTQ, es decir, a qué secuencias de un transcriptoma mapean, sin necesidad de alinear base a base los reads. Los alineamientos se pueden calcular, si fuera preciso, después del mapeo. 

Hasta la fecha se han propuesto al menos dos implementaciones del algoritmo genérico de pseudoalineamiento, que se diferencian fundamentalmente en las estructuras de datos utilizadas, como se explica en estos blogs (1 , 2) [en el segundo se hace un repaso a la evolución de este tipo de algoritmos y a su nomenclatura].

En los próximos párrafos muestro dos aplicaciones usando ejecutables de RapMap y como datos un fichero (tr.fna) con 67K transcritos de Arabidopsis thaliana, ensamblados de novo con trinity, y otro (1.f1) con 89M de lecturas Illumina SE en formato FASTQ.

1. Mapeo de reads (20 cores, fichero de salida SAM)

          operación   tiempo   comando      

BWAmem       índice    1m28s   bwa-0.7.12/bwa index tr.fna    

BWAmem        mapeo    9m57s   bwa-0.7.12/bwa mem -t 20 tr.fna 1.fq > 1.sam 

RapMap       índice      58s   RapMap-0.2.2/bin/rapmap quasiindex -t tr.fna -i ./ 

RapMap        mapeo    5m56s   RapMap-0.2.2/bin/rapmap quasimap -t 20 -i ./ -r 1.fq \

                                 -o 1.sam 

Por comparación con BWA queda claro que incluso generando alineamientos SAM estos algoritmos son mucho más rápidos. 


2. Conteo de transcritos (20 cores, implementación Sailfish)

          operación   tiempo   comando      

Sailfish     índice    1m06s   SailfishBeta-0.10.0/bin/sailfish index -t tr.fna -o qindex \

                                  -p 20   

Sailfish      mapeo    1m09s   SailfishBeta-0.10.0/bin/sailfish quant -i qindex/ -r 1.fq \

                                  -p 20 -o 1.bur-0.sf --auxDir tmpsf --dumpEq --libType SF
  

Como se muestra en la tabla, la estimación de valores de expresión de transcritos por pseudoalineamiento y conteo de reads es muy eficiente. En este ejemplo se genera un archivo de salida (quant.sf) con este contenido:

Name      Length EffectiveLength TPM NumReads
TR1|c0_g1_i1 517 316.623      5.3984 214.782
TR2|c0_g1_i1 390 191.501     3.36608     81
TR3|c0_g1_i1 699 498.611      10.502 658
...

Entre las opciones del program está la posibilidad de calcular los conteos con bootstraping, como hace kallisto, aunque en la versión que yo he probado es experimental.

En la siguiente entrada veremos más aplicaciones,
un saludo,
Bruno

Integration of footprintDB and RSAT

Hi,
it's been a while since the first public release of footprintDB, our database of DNA regulatory motifs and their cognate transcription factors with annotated interface residues. This resource can be reached at:

http://floresta.eead.csic.es/footprintdb

With help from Álvaro Sebastián, who left the lab for a postdoc abroad, we have been updating it so that it currently contains 5,084 unique transcription factors, 7,685 position specific scoring matrices (PSSMs) and 19,792 DNA binding sites, resulting from the integration of 17 motif collections, including our own 3D-footprint which is still used to annotate interface protein residues, those in direct contact with nitrogen bases of DNA.

The most recent news about footprintDB is its integration with RSAT, the well known suite of tools for the analysis of regulatory sequences designed by the team of Jacques van Helden. This means that all footprintDB motifs can now be used to annotate regulatory motifs using RSAT tools such as compare-matrices.

Since we work mostly with plants, a subset of plant-related motifs is also available, which should fit the needs of researchers concerned with plant genomics. Indeed we are in charge of the plant RSAT mirror (RSAT::Plants), which currently supports 41 plant genomes, with sequences and gene models obtained from Ensembl Plants and Phytozome.

We hope these new tools will be on interest to users,
Bruno

clustal one-liner con parallel

Hola,
hoy comparto un comando que a veces utilizo cuando necesito calcular muchos alineamientos múltiples a partir de una colección de ficheros de secuencias en formato FASTA. Como mi máquina, igual que la de casi todos, tiene amplia RAM y muchos cores, es un trabajo ideal para parallel. Supongamos que los archivos de salida están en la carpeta 'entrada' y queremos guardar los ficheros de salida en la carpeta 'path/to/salida', y que tenemos 20 cores disponibles:

$ mkdir /path/to/salida/

$ cd entrada

$ ls -1 *fasta | parallel --gnu -j 20 ~/soft/clustal-omega-1.2.1/src/clustalo \

--threads=1 -i {} -o /path/to/salida/{} :::

Este comando pondrá a trabajar 20 cores del sistema hasta que todos los archivos FASTA de la carpeta entrada estén alineados, con ganancias de tiempo de ejecución importantes en un experimento con 100 ficheros:

| cores (-j) | time(real) | time(user) | time(sys) |
|   1        | 4m34.440s  | 4m5.180s   | 0m2.168s  |
|  10        | 0m29.358s  | 3m57.768s  | 0m2.400s  |
|  20        | 0m23.248s  | 5m6.204s   | 0m3.364s  |

Un saludo,
Bruno

mapeo fino de genes por NGS

Buenas,
esta semana copio aquí una reseña de un trabajo recientemente publicado de Carlos P Cantalapiedra, autor habitual de este blog y próximo doctor del grupo, donde se explica el proceso para localizar un loci responsable de una resistencia a infección por parte de hongos, combinando genética clásica y secuenciación de nueva generación: http://www.eead.csic.es/spreading/showspreading?Id=416

Pongo aquí una de las figuras del artículo:

Genotipo de varias líneas de cebada en torno al locus que confiere resistencia. En naranja, genotipos como los del parental resistente. En verde, genotipos como los del parental susceptible. La captura de exoma permite reducir la zona de búsqueda al punto donde se unen ambos genotipos (punto 211721 dentro del recuadro). Adaptada de https://dl.sciencesocieties.org/publications/tpg/first-look/pdf/plantgenome2015.10.0101.pdf.

La referencia del artículo completo, en inglés, es:

Cantalapiedra CP, Contreras-Moreira B, Silvar C, Perovic D, Ordon F, Gracia MP, Igartua E, Casas A. (2016) A cluster of NBS-LRR genes resides in a barley powdery mildew resistance QTL on 7HL. The Plant Genome. Early access. DOI: 10.3835/plantgenome2015.10.0101. URL.

Hasta luego,

Bruno

Calculando experimentos de secuenciación

Buenas,
hoy necesitábamos calcular cuántos individuos (de una especie monocotiledónea) podríamos secuenciar con cierta profundidad en un secuenciador Illumina, pensando en el HiSeq2500 en concreto. Al final decidimos apostar por una profundidad promedio de 80x, para is sobre seguro, usando parejas de lecturas de 2x125b. Buscando en Internet encontré rápidamente una calculadora del propio fabricante que igual algunos no conocéis y puede ayudar a hacer esto rápidamente.

Figura tomada de http://www.danielecook.com/calculate-depth-coverage-bam-file.

Vayamos con un ejemplo con la calculadora
[ http://support.illumina.com/downloads/sequencing_coverage_calculator.html ]:

0. tipo de secuenciación: DNA             [se puede elegir RNA también]
1. protocolo: whole-genome sequencing  [otras: Nextera, Truseq, custom]
2. tamaño del genoma: 320Mbp
3. profundidad deseada: 80x
4. % de duplicados: 2%                [valor por defecto]
5. instrumento: HiSeq 1500/2500

Volumen total de secuenciación requerido: 26,1Gb   [26.122.448.980b]

En mi ejemplo, usando el protocolo v4, esto equivale a 0.42 líneas o lanes, lo que significa que podría poner hasta 2 muestras por línea.

Hasta luego,
Bruno

PD Me dicen mis colegas Dave Des Marais y Pat Edger que la longitud de un genoma (de plantas en este caso)  puede estimarse aproximadamente a partir del contenido en DNA del núclo usando la fórmula long = 1C * 980.

R one-liners

Hola,
en esta entrada quería compartir unos ejemplos para aprovechar las capacidades del lenguaje R para analizar de manera rápida datos desde el terminal, espoleado por el primero de ellos, que me pasó Carlos Cantalapiedra la semana pasada. Pero vayamos al grano usando como ejemplo un fichero de salida tabular de BLAST, al que llamarmos 'all.blast'.

Receta 1. Tienes un fichero con columnas de números y quieres saber la media de una de ellas, por ejemplo la tercera. Para la la desviación estándar cambia 'mean' por 'sd':

$ cut -f 3 all.blast | Rscript -e 'data=scan(file="stdin"); mean(data)'

Receta 2. Tienes un fichero con columnas de números y quieres calcular un histograma de una de ellas, por ejemplo la tercera:

$ cut -f 3 all.blast | \

Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)' 

Se puede modificar mostrando el PDF creado inmediatamente:

$ cut -f 3 all.blast | \

Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)'; evince hist.pdf 

Receta 3. Tienes un fichero con columnas de números y quieres calcular boxplots con dos de ellas, por ejemplo la primera y la segunda:

cut -f 1,2 all.blast | \

Rscript -e 'd=matrix(scan(file="stdin"),ncol=2); pdf("box.pdf");boxplot(d[,1],d[,2])'

Espero que os sirvan de inspiración,
un saludo,
Bruno

Expresión regular de la familia de las O-fucosiltransferasas

Hola, 

el pasado 8 de febrero se publicó en la revista Nature Chemical Biology  (http://dx.doi.org/10.1038/nchembio.2019) un artículo donde se  describen las bases moleculares de la reacción de O-fucosilación. Ésta es una modificación postraduccional poco frecuente, que realizan las enzimas O-fucosiltransferasas, como nos explica la investigadora de nuestro grupo Inmaculada Yruela, una de las autoras del trabajo [reseña completa en www.eead.csic.es]:

"Esta reacción resulta esencial en algunas rutas metabólicas de los organismos eucariotas, incluidas las plantas, para mantener las funciones básicas de las células. El artículo describe el mecanismo por el cual la enzima O-fucosiltransferasa 2 (POFUT2) reconoce sin errores una secuencia de aminoácidos (TSR) de la proteína receptora y le transfiere una molécula de azúcar tipo fucosa –así resulta fucosilada–. Hasta la fecha no se conocía cómo estas enzimas reconocen y se unen a sus sustratos proteicos. La O-fucosilación es esencial para el correcto plegamiento y estabilidad del dominio TSR y el reconocimiento molecular de POFUT2."

Como se ve en el alineamiento, el dominio TSR contiene tres puentes disulfuro y una secuencia consenso en torno a los dos primeras cisteínas CX{2,3}[S|T]CX{2}G , lo que se llama un motivo, como los del repositorio Prosite:

Fragmento de un alineamiento múltiple de dominios TSR, tomado de http://dx.doi.org/10.1038/nchembio.2019

Tras alinear algunas secuencias de ejemplo y perfeccionar la expresión regular, ésta se empleó para identificar todas las secuencias con dominios TSR en los proteomas de Homo sapiens y Caenorhabditis elegans. Más generalmente, el siguiente trozo de código permite localizar dentro de un fichero FASTA todas las secuencias reconocidas por una expresión regular:

 #!/usr/bin/perl  
 use strict;  
   
 # script that takes a FASTA file(s) and scans all protein sequences   
 # looking for matches of a chosen motif expressed as a regular expression  
 # 11042015: edited following comments by JF
   
 # constant as indices to access read_FASTA_file arrays  
 use constant NAME => 0;   
 use constant SEQ => 1;  
   
 #my $motifRE = '.*C[^C]{0,21}.C[^C]{2,6}[S|T].C[^C]{4,75}C[^C]C[^C]{4,15}C.*'; 
 my $motifRE = 'C[^C]{0,21}.C[^C]{2,6}[ST].C[^C]{4,75}C[^C]C[^C]{4,15}C';
   
 if(!@ARGV){ die "# usage: $0  ... \n"; }  
 else{ print "# motif: $motifRE\n"; }  
   
 my ($n_of_matches,@matches) = (0);  
 foreach my $infile (@ARGV)  
 {  
   next if($infile !~ /.fa/);  
   my (@FASTA,$start,$end,$l);  
   my $n_of_sequences = -1;  
   open(FASTA,"<$infile") || die "# cannot read $infile $!:\n";  
   while()  
   {  
    next if(/^$/);  
    if(/^\>(.*?)[\n\r]/)  
    {  
      $n_of_sequences++;   
      $FASTA[$n_of_sequences][NAME] = $1;  
    }                
    else  
    {  
      s/[-\s\n.]//g; #s/[\s|\n|\-|\.]//g;  
      $FASTA[$n_of_sequences][SEQ] .= $_;  
    }  
   }  
   close(FASTA);  
   
   foreach my $seq ( 0 .. $n_of_sequences )  
   {  
    while($FASTA[$seq][SEQ] =~ m/($motifRE)/gi)   
    {  
      #$start=length($`);  
      #$l=length($1);  
      #$end=$start+$l;  

      $start = $-[1];
      $end   = $+[1];
      $l     = $end - $start + 1;

      print ">$FASTA[$seq][NAME] match=($start,$end,$l)\n$1\n";           
    }  
   }  
 }  

Un saludo,
Inma y Bruno

Word + bioinformática = CRISPR

Hola,
este será un artículo muy breve, porque tengo poco que escribir yo mismo, a parte de volver a comprobar una vez más que la ciencia básica, tan denostada a menudo, es a menudo fundamental para los avances más espectaculares de la ciencia aplicada.  En vez de contar yo aquí una historia hoy os invito a leer la increíble historia del descubrimiento de los sistemas CRISPR-Cas9, que ya están revolucionando la ciencia y pronto la medicina.

La fuente original, Cell, bajo subscripción:
http://www.cell.com/cell/fulltext/S0092-8674%2815%2901705-5

Un resumen en español en El Confidencial:
http://www.elconfidencial.com/tecnologia/2016-01-16/crispr-francis-mojica-charpentier-doudna-edicion-genomica_1136337/

Reconocimiento de una hebra sencilla de DNA por hibridación de un RNA guía de 32 bases (crRNA) por parte del complejo de proteínas Cascade. Tomado de [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4427192/]. Como se ve en nuestro repositorio 3D-footprint, la especificidad de reconocimiento es de las más altas observadas entre proteínas naturales: http://floresta.eead.csic.es/3dfootprint/complexes/4qyz_ABCDEFGHIJ.html

 Ah, y por cierto, esta historia una cura de humildad para los que nos dedicamos a inventar y probar complicados algoritmos. A veces, la herramienta "Buscar y reemplazar" de Word es todo lo que hace falta para hacer descubrimientos en Biología, como a veces me recuerda mi colega Ana Casas, y eso también es Biología computacional,
un saludo,
Bruno

Bioinformática en C (y perl6 async)

Hola, y feliz año!
tras la última entrada, y antes de que escriba algo más completo sobre perl6,
parece que la gestión de tareas asíncronas tiene mucho potencial...

Pero vamos a lo que vamos. Hoy quería compartir un documento sobre cómo programar en C con los compiladores actuales, porque algunos todavía seguimos usando sintaxis anticuada y esto es relevante en nuestro campo, dado que hay muchas aplicaciones bioinformáticas escritas en C, normalmente por cuestiones de eficiencia.

El documento en cuestión es: https://matt.sh/howto-c

y ofrece un montón de trucos y recomendaciones para modernizar y simplificar la escritura de código en C, con ejemplos sencillos como éste para vectores dinámicos:

 uintmax_t arrayLength = strtoumax(argv[1], NULL, 10);
    void *array[arrayLength];

 /* Hace innecesario liberar el puntero al final.

    Ojo: debe ser un tamaño razonable para el sistema donde se ejecute. */

Ha sido traducido al español aquí y comentado/criticado aquí y en menéame,
un saludo,
Bruno