contenido GC de un fichero FASTA

Hola,

una pregunta habitual cuando analizas un fichero de nucleótidos, por ejemplo un ensamblaje de un genoma, es qué porcentaje GC tiene.  Asumiendo que el fichero está en formato FASTA, podemos obtener fácilmente ese valor con un mini-programa (one-liner) escrito en lenguaje perl. Por ejemplo, para el genoma comprimido de Brachypodium distachyon obtenido de Ensembl Plants, podríamos obtenerlo así:

zcat Brachypodium_distachyon.Brachypodium_distachyon_v3.0.dna.toplevel.fa.gz | \
   perl -lne 'if(!/^>/){ $SQ=uc($_); while($SQ =~ /([ACTG])/g){ $stat{$1}++; $tot++ } } 
   END{ printf("%%GC=%1.1f\n",100*($stat{"G"}+$stat{"C"})/$tot);  
      foreach $nt (keys(%stat)){ print "$nt\t$stat{$nt}" } }'

%GC=46.4
A	72549289
T	72561114
C	62839311
G	62789747

Si quieres calcular el %GC solamente para ciertas regiones del genoma entonces puedes codificarlas en un fichero BED y usar bedtools nuc, como se explica en https://www.biostars.org/p/47047

Hasta pronto,

Bruno

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