Probamos miniprot para mapear proteínas sobre genomas

Hola, 

hoy escribo a mi regreso del X Congreso Nacional de Mejora Genética de Plantas, donde se habló y mucho de herramientas de genómica computacional.

Justo esos días me enteré de la liberación de las primeras versiones de miniprot, un programa de Heng Li, el creador de minimap, del que ya hablamos aquí comparándolo con BLASTN. 

Esto me recordó que hace unos años, mientras Carlos Cantalapiedra empezaba a desarrollar BARLEYMAP, nos preguntamos qué programas había disponibles para mapear secuencias de genes, tránscritos y proteínas sobre genomas. Para los dos primeros tipos de secuencias encontramos GMAP, que adoptamos para nuestro nuevo software, pero para el tercero no encontramos ninguno que nos gustara del todo más allá de BLASTX y spaln. Justo para eso es miniprot. Lo he probado con una proteína de cebada:

#installation
git clone https://github.com/lh3/miniprot.git
cd miniprot/
make

# index barley genome
./miniprot -t 8 -d GCA_904849725.1_MorexV3.mpi GCA_904849725.1_MorexV3.fna

# example: map protein HORVU3Hr1G095240
./miniprot --gff GCA_904849725.1_MorexV3.fna HvOs2/HORVU3Hr1G095240.pep.fa

Cómo veis pedí la salida en formato GFF y me la ha devuelto precedida del mismo resultado en formato PAF, que incluye CIGARs y en la columna 10 el número nucleótidos alineados (477 en el ejemplo) y :

##gff-version 3
##PAF	transcript:HORVU3Hr1G095240.2	159	0	159	...
chr3H_LR890098.1	miniprot	mRNA	577160564	577200549	...
chr3H_LR890098.1	miniprot	CDS	577200365	577200549	...
chr3H_LR890098.1	miniprot	CDS	577162212	577162293	...
chr3H_LR890098.1	miniprot	CDS	577161755	577161804	...
chr3H_LR890098.1	miniprot	CDS	577161575	577161671	...
chr3H_LR890098.1	miniprot	CDS	577160567	577160629	...
chr3H_LR890098.1	miniprot	stop_codon	577160564	577160566	...

Lo que coincide con los resultados de BARLEYMAP cuando busco la secuencia de nucleótidos (CDS) correspondiente:

HORVU3Hr1G095240.2	chr3H_LR890098.1 577160564 577200549

 

El manual está en https://lh3.github.io/miniprot/miniprot.html , si encontráis errores podéis comunicarlos en https://github.com/lh3/miniprot/issues

 

Hasta pronto,

Bruno

Cortar secuencias desde un fichero GFF y un FASTA del genoma

Hola,

hoy cuento cómo cortar subsecuencias de un genoma correspondientes a elementos anotados en un fichero GFF asociado, que tienen este aspecto, con columnas separadas por tabuladores. A diferencia del formato BED, las coordinadas aquí empiezan a contar en 1:

1  proveedor  gene  16399  20144  .  +  .  ID=500;...
1  proveedor  mRNA  16399  20144  .  +  .  ID=500-01;Parent=500;...
1  proveedor  exon  16399  16976  .  +  .  Parent=500-01;Name=500-01-E1;...
1  proveedor  CDS   16599  16976  .  +  0  ID=CDS:500-01;Parent=500-01;...
1  proveedor  exon  17383  17474  .  +  .  Parent=500-01;Name=O500-01-E2;...
1  proveedor  CDS   17383  17474  .  +  0  ID=CDS:500-01;Parent=500-01;...
...

Como se puede ver, cada gen es un intervalo dentro de un cromosoma o contig. En este ejemplo, adaptado de Ensembl Plants, el gen con identificador 500 está en el segmento 16399-20144 del cromosoma 1 en la hebra directa (+).

Para cortar la secuencia de los genes en un fichero GFF podemos hacer lo siguiente, con ayuda de bedtools getfasta, que deberás instalar previamente. Los identificadores de cromosomas/contigs deben coincidir en ambos ficheros:

# por acortar las URLs
ENSEMBL="http://ftp.ebi.ac.uk/ensemblgenomes/pub/release-51/plants"

# descargar fichero FASTA del genoma, o del chr en este caso
wget ${ENSEMBL}/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa.gz .
gunzip Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa.gz

# descargar y descomprimir fichero GFF
wget ${ENSEMBL}/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3.gz .
gunzip Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3.gz

# extraer solamente los elementos de tipo "gene"
perl -lane 'print if($F[2] eq "gene")' Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3 > \ 
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.gff3

# cortar los genes y guardarlos en un fichero FASTA
bedtools getfasta -fi Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    -bed Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.gff3 \
    -fo Oryza_sativa.IRGSP-1.0.51.chromosome.1.gene.fa

Si necesitamos cortar las secuencias codificantes (CDS) necesitamos instalar gffread desde https://github.com/gpertea/gffread (mirar artículo aquí). Con este programa puedes extraer los CDS como nucleótidos o aminoácidos, pero tiene muchas más opciones:

# secuencia de exones codificantes concatenados, una por tránscrito/mRNA
/path/to/gffread-0.12.7.Linux_x86_64/gffread -x cds.fna \
    -g Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3

# traducción de exones codificantes concatenados, una por tránscrito/mRNA
/path/to/gffread-0.12.7.Linux_x86_64/gffread -y cds.faa \
    -g Oryza_sativa.IRGSP-1.0.dna.chromosome.1.fa \
    Oryza_sativa.IRGSP-1.0.51.chromosome.1.gff3

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