Comprime y consulta ficheros FASTA con agc

Hola,  hoy voy a escribir sobre el software AGC (Assembled Genomes Compressor), escrito en C++, que sirve para comprimir ensamblajes genómicos en múltiples ficheros FASTA (un pangenoma) y luego extraer de manera eficiente secuencias de nucleótidos arbitrarias o contigs completos. El artículo que lo describe está en https://doi.org/10.1093/bioinformatics/btad097 y el código en https://github.com/refresh-bio/agc. 

AGC comprime los ensamblajes/ficheros FASTA en varias etapas usando por defecto k=31:

Principales etapas de la compresión: (a) selección de k-meros distribuidos de manera uniforme en genoma de referencia cada 60Kb, splitters;  (b)  compresión del genoma de referencia tras partirlo en grupos y segmentos; (c y d)  compresión del resto de genomas. Los segmentos flanqueados por los mismos splitters se comprimen juntos; los de referencia con zstd y los demás con LZSS para marcar las diferencias respecto al de referencia correspondiente para comprimir solamente esas partes con zstd. Figura tomada de https://doi.org/10.1093/bioinformatics/btad097

En nuestro caso, probamos AGC con 20 genomas de cebada del pangenoma V2, cada uno de unos 4GB, obteniendo una estructura de datos comprimida de 3.4GB con el siguiente comando:

$ agc create <ref.fa> <genoma2.fa> ... <genoma20.fa> assemblies.agc

# nota: se pueden añadir nuevos genomas más tarde con agc append 

Comprobamos el contenido del archivo comprimido:

$ agc listset assemblies.agc
Akashinriki
B1K-04-12
Barke
Chiba
Du_Li_Huang
GoldenPromise
HOR_10350
HOR_13821
HOR_13942
HOR_21599
HOR_3081
HOR_3365
HOR_7552
HOR_8148
HOR_9043
Hockett
Igri
MorexV3
OUN333
Planet
 

Finalmente, lo más interesante, podemos obtener secuencias de interés por medio de sus coordenadas:

$ agc getctg assemblies.agc chr1H@Barke:1-20 

>chr1H sampleName=Barke:1-20
ATGCTATTAGTCACTAATTT

Toda la documentación y más ejemplos están en https://github.com/refresh-bio/agc, que lo disfrutéis,

Bruno

Notas de la 3ª asamblea general de la conexiónBCB del CSIC

Hola, comparto aquí mis notas de algunas charlas que escuché en la asamblea de la conexiónBCB los días 21 y 22 de abril en Madrid.

 

Francesc Montardit (Chesco) habló sobre sus progresos en "Genomic footprint of olive 

orchard management on a plant groundcover", todavía sin publicar.

 

Manolo Fernández Pérez presentó PhyloCNN, que analiza imágenes de árboles filogenéticos 

mediante redes convolucionales.

 

Pablo Herrera Nieto (CNB-CSIC) está diseñando y evaluando complejos de proteínas: 

 

Fabian Jetzinger (I2SYSBIO) explica su trabajo prepublicado con long reads en torno al 

problema de unir o no diferentes muestras antes de definir isoformas (join, call [more novelty] vs 

call, join [less novelty, more conservative]). Ven diferentes comportamientos con diferentes 

herramientas (IsoQuant, FLAIR) que creen se deben a umbrales relativos a la expresión total 

de un experimento. Recomienda join & call con FLAIR antes de lanzarse a estudiar nuevas 

isoformas raras.

 

Juan Martín Menor de Gaspar (CBGP-INIA) habla de multiomic integration en programas de 

mejora a largo plazo donde se imputan datos perdidos, en concreto missing covariances entre 

indivíduos en base a su pedigrí [covcomb], antes de hacer GS. 

 

Julen Santiago presenta "Development of an LLM and Retrieval Augmented Arch (RAG)-based 

Conversational Assistant for the Conexión BCB Platform" with https://weaviate.io para construir 

un knowledge graph e inmersiones especializadas (S-BioBERT) para calcular distancias. Para 

validar preparan un conjunto de consultas y respuestas correctas. Le preguntan si ha probado 

si el sistema alucina, pidiéndole cosas imposibles, y si está suficientemente protegida la 

privacidad de los datos. 

 

Carlos O Sorzano (CNB) habla de su trabajo con cryoEM y hetSIREN 

(https://www.nature.com/articles/s41467-025-59135-0).

 

Isabel Díaz (adjunta VRI) nos recuerda que aprovechemos los correos de programas 

internacionales del CSIC para presionar y conseguir topics en futuras convocatorias. 57% 

fondos VRI van a cooperación. Distingue infraestructuras de investigación (servicio) y 

tecnológicas (para escalar y llevar al mercado):

 

 

Enrique Bernal Delgado (IACS) habla de la European Open Science Cloud y ciencia abierta. 

Es epidemiólogo y del board of directors de EOSC-A. Pilota la propuesta ES para un nodo 

propio en  EOSC.   

 

Rocío Tuda, oficina del dato del CSIC, SGAI presenta sabio.csic.es , haya, pino, caoba y 

resume la estrategia del dato, que también incluye a https://digital.csic.es

 

Manuel Ferrer (ICP) presenta la conexión microbioma y habla de la brecha entre recursos 

genómicos de microbios y recursos biológicos (muestras, ecosistemas; muy por detrás).

 

Jesús Cerquides presenta https://aihub.csic.es y habla de algunos hitos muy recientes de la IA, 

como el benchmark BioAgent, y la tendencia a pensar que los científicos dejaremos de hacer 

ciencia y pasaremos a supervisarla, apoyándose en experiencias como ésta, donde un sistema

IA envío un abstract a un congreso que fue aceptado.  Resume el impacto que está teniendo 

ya la IA en la biología computacional y la ciencia en general:

 

Le dedica un tiempo a hablar del impacto inmediato del copilotaje en programación:

 

Nos propone una agenda en torno a esto para la conexión BCB:

 

Hasta pronto, Bruno 

 

Lecturas con SEQ y QUAL de diferente longitud en FASTQ

Hola, si alguna vez te encuentras un mensaje como éste:

[E::sam_parse1] SEQ and QUAL are of different length

Se debe a que el fichero FASTQ en cuestión contiene líneas de secuencia (SEQ) y calidad (QUAL) de diferente longitud, lo cual viola el formato y es posiblemente un error. El siguiente oneliner te sirve para encontrar las secuencias problemáticas:

$ zcat file.fastq.gz | perl -lne 'if($.%4==1){$n=$_}elsif($.%4==2){$l=length($_)}elsif($.%4==0){print $n if(length($_) != $l)}' > lista.errores.txt

Luego puedes eliminar las secuencias del fichero FASTQ como se explica por ejemplo aquí con seqtk, hasta luego,

Bruno 

Obtén taxonomías del NCBI desde el terminal

Hola, en una entrada anterior comentaba cómo aprendimos a descargar datos del NCBI desde el terminal, usando el binario datasets, en mi caso para Linux. 

Hoy os cuento cómo se puede utilizar para obtener de manera sencilla la taxonomía completa correspondiente a un taxonID, como los que utiliza la base de datos https://www.ncbi.nlm.nih.gov/taxonomy . Lo explico con un ejemplo, buscando en dicho portal la primera especie que se me ocurre:

• En la página de resultados (ver imágen) paso el ratón sobre el nombre de la especie y me fijo en el navegador a qué URL me lleva, en este caso https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=56046 

• Deduzco que el taxonID de esta especie es el 56046
• Invoco el binario datasets con este identificador, obteniendo resultados en formato JSON:

       $ datasets summary taxonomy taxon 56046

 

{"reports": [{"query":["56046"],"taxonomy":{"children":[928733,928732,928731],"classification":{"class":{"id":58019,"name":"Pinopsida"},"domain":{"id":2759,"name":"Eukaryota"},"family":{"id":3318,"name":"Pinaceae"},"genus":{"id":3319,"name":"Abies"},"kingdom":{"id":33090,"name":"Viridiplantae"},"order":{"id":1446380,"name":"Pinales"},"phylum":{"id":35493,"name":"Streptophyta"},"species":{"id":56046,"name":"Abies pinsapo"}},"curator_common_name":"Spanish fir","current_scientific_name":{"authority":"Boiss., 1838","name":"Abies pinsapo"},"current_scientific_name_is_formal":true,"genomic_moltype":"dsDNA","group_name":"seed plants","parents":[1,131567,2759,33090,35493,131221,3193,58023,78536,58024,1437180,58019,3313,2821352,1446380,3318,3319],"rank":"SPECIES","tax_id":56046}}],"total_count": 1}

Cuando nos interesa solamente un taxonID es realmente más cómodo hacerlo en la Web, como se ha explicado. Esta operación tiene más sentido en el terminal cuando queremos hacerlo para toda una lista. Por ejemplo, podemos hacer una búsqueda con BLASTN contra la colección core_nt, obteniendo el taxonID de cada resultado en la columna 13:

$ blastn -query test.fna -db core_nt -outfmt "6 std staxids" -out results.tsv 

Con un comando como el siguiente podemos repetir la llamada a datasets para ese fichero:

$ perl -lane 'foreach $i (split(";",$F[12])){ $t=$see{$i}||`datasets summary taxonomy taxon $i`; printf("%s\t%s\t\%s",$F[0],$i,$t); $see{$i}=$t }' results.tsv

Hasta pronto, Bruno 

 

footprintDB March 2026 version

Hi, we just updated the motifs, transcription factors and sites in the database footprintDB. 

This version adds:

 

1)  Motifs inferred from  protein-DNA complexes at the Protein Data Bank, added to 3d-footprint by 19/08/2025; note the complexes are also used to annotate interface residues  of all transcription factors (TFs).

 

2) New plant data at EEADannot, with TFs assigned to Plant-TFClass families (see repo).

The current contents include:

 

The footprintDB motifs have also been synced with RSAT (see repo), see you soon,

Bruno