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

conjunto diferencia entre listas con Perl

Hola de nuevo,
sirva esta entrada para compartir una receta eficiente para calcular el conjunto diferencia entre dos listas o arrays en lenguaje Perl5.

Conjunto diferencia, tomado de https://es.wikipedia.org/wiki/Diferencia_de_conjuntos.

Para ello podemos definar la siguiente subrutina, tomada del módulo Array::Utils:

my @a = 0..10000; 
my @b = 5000..10000; 

array_minus(@a, @b);

sub array_minus(\@\@) {
    my %e = map{ $_ => undef } @{$_[1]};
    return grep( ! exists( $e{$_} ), @{$_[0]} ); 
}

Podéis ver otras alternativas en reddit, incluyendo soluciones en Perl6 y python,
un saludo,
Bruno

estructuras del PDB parecidas en secuencia (REST)

Hola,
hace unos días un usuario preguntaba en la lista de usuarios del Protein Data Bank (PDB) cómo usar la interfaz de servicios REST, documentada en https://www.rcsb.org/pdb/software/rest.do, para hacer búsquedas BLAST.

Mientras otro usuario compartía un script escrito en Python3, llamado sequenceSimilarity.py, que requiere mmtf-pyspark para hacer consultas PSI-BLAST en tiempo real contra el PDB, a mi me llamó la atención la simplicidad del servicio sequenceCluster, que para cualquier cadena de una estructura depositada en el PDB permite obtener el clúster de secuencias con un porcentaje de identidad controlado por el usuario:

https://www.rcsb.org/pdb/rest/sequenceCluster?cluster=70&structureId=9ant.A

Esto devuelve una lista de cadenas de estructuras PDB en formato XML que podemos procesar por ejemplo como se explica en

https://developer.atlassian.com/server/fisheye-crucible/writing-a-rest-client-in-perl

Hasta luego, buen verano,
Bruno

updated footprintDB motifs

Hi,
I have just updated the format of footprintDB motifs used in RSAT (www.rsat.eu), which can de downloaded from http://floresta.eead.csic.es/footprintdb/download. This was done to ensure compatibility with compare-matrices in RSAT. Now, motifs look this:

AC  AY750993/VRN1/EEADannot
XX
ID  AY750993:VRN1:EEADannot
XX
NA  AY750993
XX
DE  AY750993
XX
OS  Hordeum vulgare cv Strider
XX
BF  10127;
XX
P0      A      C      G      T
01     19     72      7      2      C
02     10     71      7     12      C
03     57     17     14     12      a
04     56      6     32      6      r
05     94      2      4      0      A
06     84      1      3     12      A
07     79      7      7      7      A
08     57      4      9     30      w
09     24      0     76      0      G
10     15      0     85      0      G
XX
LN  http://floresta.eead.csic.es/footprintdb/index.php?motif=AY750993:VRN1:EEADannot
XX
RN  [1];
RL  Deng W, Casao Mc, Wang P, Sato K, Hayes PM, Jean Finnegan E, Trvaskis B (2015) Direct links between the vernalization response and other key 
traits of cereal crops. Nat Comm 6:5882
RN  [2];
RX  PUBMED: 24234003
RL  Sebastian, A., Contreras-Moreira, B. footprintDB: a database of transcription factors with annotated cis elements and binding interfaces. Bio
informatics 30, 258-65 (2014).
XX

Note that the ID field was added and that a URL pointing to the relevant database entry identified with it is attached. The accession code replaces ':' with '/' to ensure that downstream RSAT analyses can be carried out just fine,
cheers,
Bruno

Identificar tránscritos no codificantes

Hola,
recientemente he leído diferentes trabajos sobre cómo solamente una fracción de los tránscritos totales humanos realmente son codificantes. Estas observaciones tienen consecuencias prácticas que supongo se pueden extrapolar a otras especies.

1. Cómo identificar la isoforma principal de un gen
La base de datos APPRIS (http://appris-tools.org) aplica una serie de filtros para identificar las isoformas principales de cada gen humano en base a la combinación de varios criterios (leer artículo):

• conservación de la estructura de exones
• evolución no neutral
• alineamiento sin inserciones con estructuras homólogas
• conservación de residuos funcionales
• alineamiento completo con secuencias de otros vertebrados

Anotación de 3 isoformas según APPRIS, tomada de https://academic.oup.com/nar/article/46/D1/D213/4561658

2. Como identificar un transcrito no codificante
Una vez hemos ensamblado tránscritos de uno o más tejidos o condiciones puede ser útil clasificarlos como codificantes o no. En un trabajo nuestro (leer aquí) lo hacíamos con CPC y el script transcripts2cdsCPP.pl de GET_HOMOLOGUES-EST. Ahora, un trabajo reciente (leer aquí) propone los siguientes criterios, que comento en algunos casos:

• El tránscrito debe abarcar al menos un intrón y tener un nivel de expresión > 1 tránscrito por millón (TPM).
• Si sólo comprende un exón debe expresarse al menos igual que los  tránscritos mejor descritos (TPM > 13.87 en humanos).
• No debe estar contenido en otro tránscrito.
• Debe codificar un marco de lectura (ORF) de al menos 60 aminoácidos. OJO: esto podría dejar fuera proteínas pequeñas importantes.
• El ORF no debe solapar con elementos repetidos/transposones (LINe, LTR, etc) ni loci rRNA.
• El E-valor producido por BLASTX al comparar el transcrito con proteínas de mamíferos en GenBank y UniProt debe ser < 10E-15. OJO: el E-valor cambiará a medida que crezcan las bases de datos, este criterio debería expresarse mejor como cobertura o bit score. Ver siguiente criterio.
• El mejor ORF del tránscrito debe alinear al menos el 75% de otras proteínas conocidas, para eliminar pseudogenes, que suelen estar truncados.
• Si dos tránscritos de hebras contrarias solapan, nos quedaremos con el que se parezca a proteínas con función conocida.

Espero que esto os sea útil, un saludo,
Bruno