blast eficiente con alias de nr

Hola,
en esta entrada quisiera mostraros una manera de buscar secuencias de proteínas dentro de la base de datos estrictamente no redundante nr del NCBI, restringiendo la búsqueda a grupos taxonómicos arbitrarios, sin necesidad de duplicar secuencias en otros archivos. De hecho, por su gran tamaño, 45 millones de secuencias a dia de hoy, desde hace tiempo nosotros ya no trabajamos con el archivo FASTA de nr, si no con la base de datos ya formateada para blast que descargamos regularmente con update_blastdb.pl desde ftp://ftp.ncbi.nlm.nih.gov/blast/db .

La receta para realizar búsquedas sobre subconjuntos de nr la hemos sacado de https://www.biostars.org/p/6528, y aquí os mostraré un ejemplo con secuencias de plantas:

1) Si no la tienes ya descárgate con el script update_blastdb.pl una copia de nr. Ojo, son una serie de archivos que en total ocupan 42G:
$ perl update_blastdb.pl nr


2) Encontrar el identificador taxonómico que mejor represente al subconjunto de nr que quieres definir, en nuestro caso Viridiplantae:
http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=33090

3) Buscar todas las proteínas de plantas anotadas en el NCBI y exportar (Send to File => GI List) sus identificadores de GenBank (GIs) a un archivo, por ejemplo 'Viridiplantae_protein.gi.txt':
http://www.ncbi.nlm.nih.gov/protein/?term=txid33090[ORGN]

4) Hacer un alias de nr que contenga sólo proteínas de plantas:
$ ncbi-blast-2.2.29+/bin/blastdb_aliastool -gilist Viridiplantae_protein.gi.txt -db nr -out nr_plants -title nr_plants

En mi ejemplo produce el siguiente output:

Converted 3741334 GIs from Viridiplantae_protein.gi.txt to binary format in nr_plants.p.gil
Created protein BLAST (alias) database nr_plants with 2041649 sequences


Ahora veamos qué diferencias y ventajas tiene usar este subconjunto de nr en comparación con la base de datos completa:

[ NR ]
$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr

 Database: All non-redundant GenBank CDS translations+PDB+SwissProt+PIR+PRF
excluding environmental samples from WGS projects
           45,360,603 sequences; 16,206,973,370 total letters

Query= test

Length=312
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

dbj|BAK07057.1|  predicted protein [Hordeum vulgare subsp. vulgare]    156    8e-47
dbj|BAJ89957.1|  predicted protein [Hordeum vulgare subsp. vulgar...  62.0    5e-10
ref|WP_014745933.1|  asparagine synthase [Tistrella mobilis] >ref  42.7    0.020
ref|WP_006720063.1|  hypothetical protein [Collinsella stercoris]...  35.0    5.0 
ref|XP_001817163.2|  hypothetical protein AOR_1_2924174 [Aspergil...  35.0    5.3 
dbj|BAE55161.1|  unnamed protein product [Aspergillus oryzae RIB40]   35.0    6.9 
ref|WP_021178774.1|  Phosphate transport system permease protein ...  34.7    7.6 
ref|WP_021807128.1|  Phosphate transport system permease protein ...  34.7    7.9 
ref|WP_026166948.1|  glucan biosynthesis protein D [Novosphingobi...  34.3    8.8

[...]

real    5m28.923s
user    5m24.564s
sys    0m3.596s

[ NR_plants ]

$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr_plants

Database: nr_plants

           2,041,649 sequences; 772,729,335 total letters

Query= test

Length=312
                                                                      Score     E

Sequences producing significant alignments:                          (Bits)  Value

dbj|BAK07057.1|  predicted protein [Hordeum vulgare subsp. vulgare]    156    4e-48
dbj|BAJ89957.1|  predicted protein [Hordeum vulgare subsp. vulgar...  62.0    2e-11
gb|EMS57556.1|  hypothetical protein TRIUR3_28581 [Triticum urartu]   33.1    0.95 
ref|XP_002311179.2|  hypothetical protein POPTR_0008s05850g [Popu...  33.5    1.1  
ref|XP_004509944.1|  PREDICTED: uncharacterized protein LOC101513...  32.0    3.0  
ref|XP_004509945.1|  PREDICTED: uncharacterized protein LOC101513...  32.0    3.0  
gb|EMT15944.1|  hypothetical protein F775_17477 [Aegilops tauschii]   32.0    3.3  
dbj|BAK00277.1|  predicted protein [Hordeum vulgare subsp. vulgare]   32.0    3.7  
ref|XP_004234300.1|  PREDICTED: uncharacterized protein LOC101244...  31.2    4.7  
ref|XP_006358717.1|  PREDICTED: LRR receptor-like serine/threonin...  31.2    5.2  
ref|XP_002316273.2|  hypothetical protein POPTR_0010s20870g [Popu...  31.2    6.6  
ref|XP_003571199.1|  PREDICTED: probable tocopherol cyclase, chlo...  30.8    8.4  
ref|XP_008451998.1|  PREDICTED: respiratory burst oxidase homolog...  30.4    9.9

[...]

real    0m16.805s
user    0m16.332s
sys    0m0.372s

Como se puede ver en ambos casos las dos primeras secuencias encontradas son las mismas, con la misma puntuación en bits y distinto E-valor, al ser universos de búsqueda de tamaños  distintos.

Además,  el tiempo de búsqueda con un sólo procesador es del orden de 20x más rápido con el alias nr_plants, y a la vez el consumo de RAM en tiempo real mucho menor.

Hasta luego,
Bruno

regenerando DNA degenerado

Hola,
durante la reciente visita de mi colega Pablo Vinuesa al laboratorio pasamos ratos escribiendo código en Perl con el fin de diseñar de manera automática, para grandes conjuntos de secuencias de DNA previamente alineadas, parejas de primers (cebadores de PCR) adecuadas para estudios de microbiología ambiental, siguiendo principios bien conocidos ya. En cuanto probamos el código con unos cuantos ejemplos observamos que a menudo los primers para algunas regiones de los alineamientos múltiples eran degenerados. Como se muestra en la figura, eso significa que algunas posiciones de los cebadores no podían ser nucleótidos únicos si no que tenían que ser combinaciones de 2 o más para poder hibridarse con las secuencias utlizadas para su diseño.

Pareja de primers degenerados que definen un amplicón a partir de secuencias de DNA alineadas. De acuerdo con la nomenclatura IUPAC, el de la izquierda (fwd) puede escribirse como GAYTST y el de la derecha (rev) como GHGKAG. Figura tomada de http://acgt.cs.tau.ac.il/hyden.

A la hora de evaluar parejas de primers nos encontramos con que los degenerados en realidad tendrían que ser examinados por medio de cada una de las secuencias implícitas en el código degenerado. Y por tanto, tuvimos que buscar una manera de regenerar las secuencias de nucleótidos correspondientes a un primer degenerado, que es de lo que va esta entrada. El siguiente código en Perl hace precisamente esto:

 #!/usr/bin/perl  
 use strict;  
 my $degenerate_sequence = $ARGV[0] || die "# usage: $0 <degenerate DNA sequence>\n";  
 my $regenerated_seqs = regenerate($degenerate_sequence);  
 foreach my $seq (@$regenerated_seqs)  
 {  
    print "$seq\n";  
 }  
 # returns a ref to list of all DNA sequences implied in a degenerate oligo  
 # adapted from Amplicon (Simon Neil Jarman, http://sourceforge.net/projects/amplicon)  
 sub regenerate  
 {  
    my ($primerseq) = @_;  
    my %IUPACdegen = (   
    'A'=>['A'],'C'=>['C'], 'G'=>['G'], 'T'=>['T'],   
    'R'=>['A','G'], 'Y'=>['C','T'], 'S'=>['C','G'], 'W'=>['A','T'], 'K'=>['G','T'], 'M'=>['A','C'],  
    'B'=>['C','G','T'], 'D'=>['A','G','T'], 'V'=>['A','C','G'], 'H'=>['A','C','T'],   
    'N'=>['A','C','G','T']   
    );  
    my @oligo = split(//,uc($primerseq));   
    my @grow = ('');  
    my @newgrow = ('');  
    my ($res,$degen,$x,$n,$seq);  
    foreach $res (0 .. $#oligo){  
       $degen = $IUPACdegen{$oligo[$res]};   
       if($#{$degen}>0){  
          $x = 0;  
          @newgrow = @grow;  
          while($x<$#{$degen}){  
             push(@newgrow,@grow);  
             $x++;     
          }     
          @grow = @newgrow;  
       }  
       $n=$x=0;   
       foreach $seq (0 .. $#grow){  
          $grow[$seq] .= $degen->[$x];   
          $n++;  
          if($n == (scalar(@grow)/scalar(@$degen))){  
             $n=0;  
             $x++;  
          }     
       }  
    }  
    return \@grow;     
 }            

Podemos probarlo con las secuencias de la figura:

$ perl degenprimer.pl GAYTST
GACTCT
GATTCT
GACTGT
GATTGT

 y

$ perl degenprimer.pl GHGKAG
GAGGAG
GCGGAG
GTGGAG
GAGTAG
GCGTAG
GTGTAG

Hasta otra,
Bruno

postdoc computational systems biology Luxemburg

Postdoctoral Fellow in Computational Systems Biology at the Luxembourg Centre for Systems Biomedicine in collaboration with the Black Family Stem Cell Institute of the Mount Sinai School of Medicine

http://www.nature.com/naturejobs/science/jobs/425291-postdoctoral-fellow-in-systems-biology

The Computational Biology Group seeks a highly skilled and motivated Postdotoral Fellow to work on an exciting project on the application of network biology approaches to study the role of gene expression stochasticity in cell fate commitment. In particular, the selected candidate shall develop a computational model that integrates transcriptomics and epigenomics data to describe heterogeneity in embryonic stem cells and its implication in differentiation. Single cell and cell perturbation experiments will be performed in order to validate the predictions generated by the model. This project will be carried out in collaboration with Profs. Ihor Lemischka and Kateri Moore at the Black Family Stem Cell Institute of the Mount Sinai School of Medicine. The selected applicant will have the opportunity to visit the experimental labs of our collaborators at the Mount Sinai School of Medicine.

Requirements of the ideal candidate:

• Ph.D. in Bioinformatics, Computer Science, Biology or a related discipline
• Strong computational skills
• Prior experience in mathematical modelling of biological networks, especially in network inference and analysis
• A strong first-author publication record in the fields of Bioinformatics and Computational Biology
• Excellent working knowledge in English.

We offer:

• Opportunity to do applied research to medical problems within a highly dynamic research institution (LCSB) and in collaboration with internationally recognized partners
• An exciting international environment
• A very competitive salary

For further information, please contact:

Prof. Dr. Antonio del Sol, Luxembourg Centre for Systems Biomedicine
E-mail: antonio.delsol@uni.lu

posdoc bioinformática/proteómica (Oxford)

Location:	Oxford, UK
Salary:	£29,837 to £36,661 discretionary range to £40,046 p.a.
Hours:	Full Time
Contract:	Contract / Temporary

Placed on:	20th May 2014
Closes:	13th June 2014
Job Ref:	113376

Sir William Dunn School of Pathology, South Parks Road, Oxford

The Central Proteomics Facility at the Sir William Dunn School of Pathology, University of Oxford is looking for an enthusiastic team member experienced in computational and statistical methods to provide bioinformatics support to researchers using proteomics (http://www.proteomics.ox.ac.uk).

The successful candidate will work collaboratively with existing bioinformaticians, facility staff and academic researchers. The position will cover two main roles. The first is to advise research scientists on statistically rigorous experimental design and data analysis procedures. The second is to maintain and expand the range of computational and statistical tools the facility uses for the analysis of proteomics data. 

You will have a degree in statistics, computing, bioinformatics or equivalent with considerable demonstrable practical experience. You should have strong communication and organisational skills, experience with analysing high dimensional datasets, statistical programming and experience with at least one scripting language. While it is desirable to have previous experience with proteomics and mass spectrometry data this is not essential as training will be provided.

The post is available as a fixed-term contract for 2 years in the first instance. If you are interested in this role, please apply online. You will be required to upload a CV and supporting statement as part of your online application. 

Interviews will be held during the week commencing 23 June 2014.

http://www.jobs.ac.uk/job/AIU252/postdoctoral-bioinformatician-biostatistician

Cuando Blastn no es Blastn

BLAST cambió hace ya algún tiempo a mejor con su versión BLAST+ pero por el camino se olvidó de algún detalle que puede confundir a más de uno.

El antiguo BLAST se ejecutaba con el comando 'blastall':

Blastall
--------

Blastall may be used to perform all five flavors of blast comparison. One
may obtain the blastall options by executing 'blastall -' (note the dash). A
typical use of blastall would be to perform a blastn search (nucl. vs. nucl.) 
of a file called QUERY would be:

blastall -p blastn -d nr -i QUERY -o out.QUERY

The output is placed into the output file out.QUERY and the search is performed
against the 'nr' database.  If a protein vs. protein search is desired,
then 'blastn' should be replaced with 'blastp' etc.

De esta forma un alineamiento de proteínas comenzaría como 'blastall -p blastp' y uno de ácidos nucleicos como 'blastall -p blastn' y si queremos usar MEGABLAST tenemos el comando diferente 'megablast'.

Sin embargo en la 'nueva' versión BLAST+, se separaron el alineamiento de proteínas y el de ácidos nucleicos en dos comandos: 'blastp' y 'blastn' (ver manual). Hasta aquí todo parece lógico y normal, lo que no todo el mundo sabe es que LA OPCIÓN POR DEFECTO DE BLASTN ES MEGABLAST, si ejecutamos 'blastn -help' encontraremos lo siguiente:

 *** General search options
 -task                 'megablast' 'rmblastn' >
   Task to execute
   Default = `megablast'

Y es que no todo el mundo está interesado en la velocidad de búsqueda de alineamientos, muchos de los que todavía usamos Blastn es porque apreciamos su gran sensibilidad para detectar alineamientos. En la actualidad usamos Blast para alinear miles de secuencias en tiempos muy razonables de minutos e incluso segundos. Para ganar velocidad en alineamientos de millones de secuencias existen otras mejores alternativas como Bowtie2.

La CONCLUSIÓN de todo esto, si usamos Blastn y nos interesa la sensibilidad deberemos ejecutarlo como:

 blastn -task blastn

Si lo hacemos sin añadir esta opción estaremos ejecutando MEGABLAST y correremos el peligro de perder una gran sensibilidad y no encontrar los alineamientos que deseamos. Por ejemplo, busquemos homología entre la 2'beta microglobulina humana (NM_004048.2) y la de ratón (NM_009735.3) usando la herramienta online de Blastn con las opciones por defecto ('Highly similar sequences (megablast)'):

Sin embargo si cambiamos la opción de búsqueda a 'Somewhat similar sequences (blastn)':

La diferencia es considerable, ¿no creéis? pasamos de no encontrar similaritud a un alineamiento con E-valor de 1.5E-56!!!!