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 

 

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