mapeando genes con liftoff

Hola,

esta semana en Ensembl tuvimos que probar a mapear genes de un ensamblaje genómico de arroz sobre otro ensamblaje más reciente. Esto se llama lift-over en la literatura. Para ello probamos un software que se llama Liftoff, publicado en Bioinformatics, con código fuente en https://github.com/agshumate/Liftoff

Como resume la figura para el tránscrito humano ENST00000598723.5, Liftoff calcula por medio de minimap2 alineamientos parciales al genoma de referencia y luego calcula el grafo más corto que conecte los alineamientos:

Resumo ahora cómo instalé este software escrito en python y cómo lo probé.

## pyenv 

# set up a pyenv folder
export PYENV_ROOT="$HOME/.pyenv"
 
curl https://pyenv.run | bash
 
pyenv install 3.7.9
 
# create dedicated env
pyenv virtualenv 3.7.9 liftoff

# install Liftoff inside that environment
pyenv shell 3.7.9 liftoff
pip install --upgrade pip
pip install Liftoff

# add the following to your .bashrc:
export PATH="$PYENV_ROOT/bin:$PATH"
eval "$(pyenv init -)"
eval "$(pyenv virtualenv-init -)"

## Test run:

pyenv shell 3.7.9 liftoff
liftoff -g old.gff3 -o old.new.liftoff.gff -p 4 new.genome.fna old.genome.fasta

 

El fichero de salida es un GFF con contenido como éste:

1       Liftoff gene    2903    10817   .       +       .       ID=LOC_Os01g01010;Name=LOC_Os01g01010;Note=TBC domain containing protein, expressed;coverage=1.0;sequence_ID=1.0;extra_copy_number=0;copy_num_ID=LOC_Os01g01010_0
1       Liftoff mRNA    2903    10817   .       +       .       ID=LOC_Os01g01010.1;Name=LOC_Os01g01010.1;Parent=LOC_Os01g01010;extra_copy_number=0
1       Liftoff exon    2903    3268    .       +       .       ID=LOC_Os01g01010.1:exon_1;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    3354    3616    .       +       .       ID=LOC_Os01g01010.1:exon_2;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    4357    4455    .       +       .       ID=LOC_Os01g01010.1:exon_3;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    5457    5560    .       +       .       ID=LOC_Os01g01010.1:exon_4;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    7136    7944    .       +       .       ID=LOC_Os01g01010.1:exon_5;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    8028    8150    .       +       .       ID=LOC_Os01g01010.1:exon_6;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    8232    8320    .       +       .       ID=LOC_Os01g01010.1:exon_7;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    8408    8608    .       +       .       ID=LOC_Os01g01010.1:exon_8;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    9210    9617    .       +       .       ID=LOC_Os01g01010.1:exon_9;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    10104   10187   .       +       .       ID=LOC_Os01g01010.1:exon_10;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    10274   10430   .       +       .       ID=LOC_Os01g01010.1:exon_11;Parent=LOC_Os01g01010.1;extra_copy_number=0
1       Liftoff exon    10504   10817   .       +       .       ID=LOC_Os01g01010.1:exon_12;Parent=LOC_Os01g01010.1;extra_copy_number=0

Hasta pronto, 

Bruno

eliminar redundancia en grandes ficheros de nucleótidos (linclust)

Hola,

recientemente me vi en la necesidad de eliminar redundancia de un fichero con millones de secuencias de nucleótidos. En mi caso se trataba de secuencias repetidas del genoma de cebada (n=4.638.834), y lo intenté con CD-HIT-EST, una herramienta que he probado muchas veces, y que tenía por muy eficiente.

Sin embargo, para esta tarea, a pesar de reservar más de 10GB de RAM, no terminá en varias horas usando 20 cores.

Por tanto, me puse a buscar opciones:

• dnaclust, usa demasiada RAM
• SigClust, muy rápido, pero usa el algoritmos K-medias y por tanto le tienes que decir cuántos clusters quieres, que en mi caso es lo que quiero averiguar
• swarm, es demasiado fino separando, solamente es eficiente para eliminar copias idénticas

Mi colega Pablo Vinuesa me habló de otras opciones que se ocupan de calcular distancias entre genomas o metagenomas, un problema relacionado:

• mash
• dRep
  

Después de mucho buscar encontré linclust, cuyo código fuente está en https://github.com/soedinglab/MMseqs2. Este programa se describió originalmente para secuencias de péptidos (artículo) pero los autores, entre ellos J Soding, añadieron después la posibilidad de agregar nucleótidos. En mis manos este es la mejor opción ahora mismo.

A continuación os muestra un banco de pruebas con secuencias repetidas de la planta Arabidopsis thaliana (n=66.752), analizadas con el binario más rápido distribuido por los autores (AVX2).

El programa lo ejecuté de esta manera:

$ command time -v mmseqs/bin/mmseqs easy-linclust \

  arabidopsis_thaliana.repeats.nondeg.fasta --threads 4 \

  arabidopsis_thaliana.48.repeats.nr0.95 ./ --min-seq-id 0.95

Verás que guarda datos temporales en el directorio actual (./) que luego debes eliminar.

Lo que observé en mis pruebas con cebada y A. thaliana es que a partir de un umbral de identidad del 70% el programa converge:

umbral	RAM (kbytes)	segundos	secuencias
0.99	311164	4.31	58285
0.95	303764	4.22	54353
0.90	307804	4.10	51359
0.80	306648	4.05	48722
0.70	304736	3.85	47770
0.50	307140	3.34	46433

Hasta pronto,

Bruno

Consultas GraphQL al Protein Data Bank

Hola, 

durante el mantenimiento anual de footprintDB he descubierto que el Protein Data Bank ha cerrado su interfaz de web servicios (WS) https://www.rcsb.org/pages/webservices/rest-fetch .  Leyendo en https://data.rcsb.org/index.html#data-api veo que ahora se pueden hacer consultas REST o GraphQL, y opté por la segunda por aprender un poco.

Hay documentación en https://data.rcsb.org/migration-guide.html#legacy-fetch-api para migrar las consultas WS antiguas a GraphQL. Las nuevas consultas tienen este aspecto:

query={
 entry(entry_id:"9ANT"{
  polymer_entities{
   rcsb_polymer_entity{
    pdbx_description
   }
   rcsb_entity_source_organism{
    scientific_name
   }
   rcsb_polymer_entity_container_identifiers{
    entry_id
    auth_asym_ids
    reference_sequence_identifiers{
     database_accession,
     database_name
    }
   }
  }
  struct{
   title
  }
  rcsb_primary_citation{
   pdbx_database_id_PubMed
  }
 }
}

Y se pueden convertir a una URL como https://data.rcsb.org/graphql?query={entry(entry_id:%229ANT%22){polymer_entities{rcsb_polymer_entity{pdbx_description}rcsb_entity_source_organism{scientific_name},rcsb_polymer_entity_container_identifiers{entry_id,auth_asym_ids,reference_sequence_identifiers{database_accession,database_name}}}struct{title}rcsb_primary_citation{pdbx_database_id_PubMed}}}

 

Si haces click verás que en Chrome o Firefox no obtienes el resultado esperado, con un error:

Invalid character found in the request target. The valid characters are defined in RFC 7230 and RFC 3986

Creo que se debe a un problema del servidor Tomcat del PDB. Sin embargo, sí funciona con wget y así obtienes la salida JSON que puedes procesar con tu lenguaje favorito:

wget 'https://data.rcsb.org/graphql?query={entry(entry_id:%229ANT%22){polymer_entities{rcsb_polymer_entity{pdbx_description}rcsb_entity_source_organism{scientific_name},rcsb_polymer_entity_container_identifiers{entry_id,auth_asym_ids,reference_sequence_identifiers{database_accession,database_name}}}struct{title}rcsb_primary_citation{pdbx_database_id_PubMed}}}' -O-

Hasta pronto,

Bruno

 

 

 

Recursos en línea sobre el lenguaje Perl

Hola,  

la bioinformática se escribe en muchos lenguajes de programación. Aunque seguramente es ahora más minoritario, por el empuje de lenguajes como Python o R en nichos como el aprendizaje automático o la estadística, uno de los lenguajes que ha tenido mucho peso en nuestra disciplina es Perl. De ahí el nombre del blog.

Para mi la principal fortaleza de Perl es lo eficiente que es para manipular ficheros de texto y para describir y ejecutar tareas complejas que llamen a otras herramientas. Y por supuesto sus enormes colecciones de módulos (core & CPAN), que puedes instalar con https://metacpan.org/pod/App::cpanminus

En cualquier caso, si estás aprendiendo a programar en Perl, o necesitas entender programas o módulos escritos por otras personas, hay buenos recursos en línea:

• https://perlweekly.com tiene un boletín semanal con actualidad, novedades y noticias por separado de Perl y el language hermano Raku
• https://www.reddit.com/r/perl es un foro general en inglés con una amplia comunidad detrás
• https://www.perlmonks.org es el foro específico para preguntas y dudas sobre código y sintaxis. Otra opción similar es https://stackoverflow.com/questions/tagged/perl 
• https://perlenespanol.com si prefieres un foro en español
• https://bioconda.github.io/search.html?q=perl para ver recetas Perl en bioconda
  
• https://github.com/Perl/perl5 si quieres cotillear el código fuente del intérprete del lenguaje y su evolución (leer esto), o contribuir
  

Hasta luego,

Bruno

 

Introduction to Bioinformatics for plant breeders

Hi, Najla Ksouri and me recently put together teaching material about Bioinformatics resources and best-practices for plant-breeders.

The material can be found at https://eead-csic-compbio.github.io/bioinformatics

It comprises four sessions and the goal is for you to explore some bioinformatics resources and databases which can be useful in plant breeding and genomics, using both standalone and Web-based tools. The last two sessions are about reproducible analysis practices and software benchmarks.

session	title	required time	URL
1	Annotation of coding sequences	4h	session 1
2	Analysis of non-coding sequences	4h	session 2
3	Reproducible analysis practices	2h	session 3
4	Benchmarks	2h	session 4

 

If you spot errors please send pull requests, hope this helps some of you out there,

Bruno