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

proyección de variantes genómicas entre genomas

Cuando se acumulan diferentes versiones del mismo genoma, como pasa con la cebada, a menudo necesitaremos proyectar anotaciones de una versión a otra. Esta operación se llama lift-over en la literatura en inglés y tiene sus complicaciones, como se ven en la figura:

fuente: https://doi.org/10.12688/f1000research.14148.2

En una entrada anterior explicaba cómo hacerlo para genes, por ejemplo con el software LiftOff. Sin embargo, a veces lo que queremos mapear son SNPs, que se habían definido sobre una versión del genoma, sobre la siguiente. 

Una manera, para genomas que tengan precalculados alineamientos en UCSC o Ensembl (chain files), es usar el software BCFtools/liftover, que se puede descargar como binario o compilar, y requiere bcftools 1.20 o superior. Puedes leer más sobre esta opción en https://doi.org/10.1093/bioinformatics/btae038 y https://github.com/freeseek/score. Una importante limitación es que solamente hay chain files pare ciertas especies. Por ejemplo, para plantas puedes consultar https://ftp.ebi.ac.uk/ensemblgenomes/pub/plants/current/assembly_chain

Para cualquier pareja de genomas podemos usar una estrategia que usábamos en Ensembl Plants, consiste en cortar la secuencia flanqueante de cada SNP en el genoma1 y mapearla sobre el genoma2 con BWA mem. Esta estrategia tiene como limitación que se pierde una fracción de las variantes originales, aquellas cuyas secuencias no mapeen bien en genoma2, o que estén en regiones repetidas, pero eso no es necesariamente malo. La ventaja que tiene es que no necesitas calcular alineamientos de dos genomas completos, lo cual es complejo y puede requerir grandes cantidades de RAM. Además en todo momento controlas lo que estás haciendo y si algo sale mal lo puedes ver y tratar de corregir. Esta estrategia se describe paso a paso en: https://github.com/eead-csic-compbio/eead-csic-compbio.github.io 

Como resultado produce texto separado por tabuladores (TSV) cómo este (ver fichero completo):

1	51976	-	LR890096.1	77101	C	G
1	51988	-	LR890096.1	77089	C	G
1	51995	-	LR890096.1	77082	G	C
1	52015	-	LR890096.1	77062	C	G
1	263632	+	LR890096.1	148230	G	G
1	263634	+	LR890096.1	148232	A	A
1	263635	+	LR890096.1	148233	A	A
1	263637	+	LR890096.1	148235	T	T
1	263638	+	LR890096.1	148236	G	G
1	263646	+	LR890096.1	148244	C	C
1	263654	+	LR890096.1	148252	C	C
1	263699	+	LR890096.1	148297	C	C
1	263706	+	LR890096.1	148304	A	A
1	270084	+	LR890096.1	154681	C	C
1	270087	+	LR890096.1	154684	G	G

Un control de calidad posible es comprobar que la base de ambos genomas es la misma, aunque a veces estará un el reverso complementario, como se ve en el ejemplo para dos regiones de los cromomas 1 (genoma1) y LR890096.1 (genoma2).

Hasta pronto,

Bruno