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