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
https://gist.github.com/darencard/9497e151882c3ff366335040e20b6714
ResponderEliminarHola su Blog me ha servido bastante, me podría ayudar... Tengo una base de datos con archivos Fasta, quisiera Filtrarlos, dejar solo secuencias de 12 a 30NS. Que programa me recomendaría usar
ResponderEliminarPuedes modificar la receta 4 de https://eead-csic-compbio.github.io/perl_bioinformatica/node90.html
ResponderEliminar