Mostrando entradas con la etiqueta gffread. Mostrar todas las entradas
Mostrando entradas con la etiqueta gffread. Mostrar todas las entradas

15 de septiembre de 2021

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