Trabajar con tablas o marcos de datos con R

Última actualización:

17 de mayo de 2026

PRÁCTICA 4

Objetivos

Conocer el uso de datasets y dataformat.
Ver algún ejemplo de archivo en formato JSON.
Manipular marcos de datos sin usar hojas de cálculo.

Datos

Descargaremos un genoma de referencia del género Trypanosoma, así como su anotación funcional, en formato GFF3.

if [ ! -e  trypanosoma_annotated.txt ]; then
   ../../../bin/datasets summary genome taxon Trypanosoma --annotated --as-json-lines | \
   ../../../dataformat tsv genome --fields \
      accession,organism-name,assminfo-release-date,assminfo-level,assmstats-atgc-count,annotinfo-featcount-gene-total,assminfo-type,assminfo-status > trypanosoma_annotated.txt
fi

La línea larga que sigue a --fields especifica qué información sobre los genomas queremos ver en la tabla. Para saber qué campos se pueden seleccionar, puedes ejecutar la orden dataformat tsv genome --help. Veamos esa tabla en el informe HTML en un formato adecuado.

library('knitr')
genomas <- read.table('trypanosoma_annotated.txt',
                      header = TRUE, sep = '\t')
kable(genomas)

Assembly.Accession	Organism.Name	Assembly.Release.Date	Assembly.Level	Assembly.Stats.ATGC.Count	Annotation.Count.Gene.Total	Assembly.Type	Assembly.Status
GCF_000002445.2	Trypanosoma brucei brucei TREU927	2005-12-14	Chromosome	27984184	10221	haploid	current
GCA_000002445.1	Trypanosoma brucei brucei TREU927	2005-12-14	Chromosome	27984184	10234	haploid	current
GCA_003543875.1	Trypanosoma brucei equiperdum	2018-09-11	Chromosome	26988976	10293	haploid	current
GCF_000210295.1	Trypanosoma brucei gambiense DAL972	2009-10-14	Chromosome	22110721	9930	haploid	current
GCA_000210295.1	Trypanosoma brucei gambiense DAL972	2009-10-14	Chromosome	22110721	9961	haploid	current
GCA_000227395.2	Trypanosoma congolense IL3000	2011-08-12	Scaffold	18802389	6348	haploid	current
GCF_003719485.1	Trypanosoma conorhini	2018-11-08	Scaffold	21137671	10321	haploid	current
GCA_003719485.1	Trypanosoma conorhini	2018-11-08	Scaffold	21137671	10321	haploid	current
GCF_000209065.1	Trypanosoma cruzi	2005-08-02	Scaffold	89612356	25210	unresolved-diploid	current
GCA_000188675.2	Trypanosoma cruzi	2012-10-02	Contig	38589511	10876	haploid	current
GCA_000209065.1	Trypanosoma cruzi	2005-08-02	Scaffold	89612356	25210	unresolved-diploid	current
GCA_003177095.1	Trypanosoma cruzi	2018-05-30	Contig	87058484	29302	unresolved-diploid	current
GCA_003177105.1	Trypanosoma cruzi	2018-05-30	Contig	53271887	19115	haploid	current
GCA_003719155.1	Trypanosoma cruzi	2018-11-05	Scaffold	20565523	34344	diploid	current
GCA_003719455.1	Trypanosoma cruzi	2018-11-08	Scaffold	23852162	13503	haploid	current
GCA_013358655.1	Trypanosoma cruzi	2020-06-18	Scaffold	40788336	14352	haploid	current
GCA_015033625.1	Trypanosoma cruzi	2020-11-03	Chromosome	45527284	18708	haploid	current
GCA_015033655.1	Trypanosoma cruzi	2020-11-03	Chromosome	47195002	17651	haploid	current
GCA_002219105.2	Trypanosoma cruzi cruzi	2017-09-15	Contig	50928259	22870	haploid	current
GCA_000496795.1	Trypanosoma cruzi Dm28c	2013-11-13	Contig	27304309	11398	haploid	current
GCA_000300495.1	Trypanosoma cruzi marinkellei	2012-10-01	Contig	34226154	10127	haploid	current
GCF_001457755.1	Trypanosoma equiperdum	2016-12-15	Scaffold	26113307	9129	haploid	current
GCA_001457755.2	Trypanosoma equiperdum	2016-12-15	Scaffold	26113307	9135	haploid	current
GCF_000691245.1	Trypanosoma grayi	2014-05-15	Scaffold	20791447	10676	haploid	current
GCA_000691245.1	Trypanosoma grayi	2014-05-15	Scaffold	20791447	10676	haploid	current
GCF_022059095.1	Trypanosoma melophagium	2022-02-08	Contig	23303718	9756	haploid	current
GCA_022059095.1	Trypanosoma melophagium	2022-02-08	Contig	23303718	9768	haploid	current
GCF_003719475.1	Trypanosoma rangeli	2018-11-08	Scaffold	20977846	10589	haploid	current
GCA_003719475.1	Trypanosoma rangeli	2018-11-08	Scaffold	20977846	10589	haploid	current
GCA_000492115.1	Trypanosoma rangeli SC58	2013-10-30	Contig	15111138	7475	haploid	current
GCF_002087225.1	Trypanosoma theileri	2017-04-12	Scaffold	25688935	11312	haploid	current
GCA_002087225.1	Trypanosoma theileri	2017-04-12	Scaffold	25688935	11312	haploid	current
GCA_021307395.1	Trypanosoma vivax	2021-12-28	Contig	67823889	17806	haploid	current
GCA_019425475.1	Trypanosoma vivax	2021-07-30	Contig	23865041	10424	haploid	current
GCA_000227375.1	Trypanosoma vivax Y486	2011-08-11	Contig	24779298	4122	haploid	current

En mi caso, no me puedo descargar el paquete “kableExtra”. Elimino del código original propuesto en el guion de la práctica del Aula Virtual la parte de library(‘kableExtra’), y la sustituyo por library(‘knitr’). “kableExtra” es un paquete cosmético. La función kable la tiene que sacar R desde el paquete ‘knitr’. También eliminamos la función kable_styling.

Para conocer qué tipo de archivo ejecutamos el siguiente código:

class(genomas)

[1] "data.frame"

Ahora, sabiendo que es un dataframe, observamos el número de filas y columnas empleando el siguiente código:

dim(genomas)

[1] 35  8

Este resultado indica que el dataframe presenta 35 filas y 8 columnas.

Desafío

¿Qué diferencia hay entre las líneas 7 y 8 de la tabla anterior?

La información proviene del mismo organismo pero el identificador es diferente: en la línea 7 es GCF y en la 8 GCA. Esto indica que provienen de diferentes bases de datos. Vemos que cambia el nombre del ensamblaje. El cambio en el nombre por una cuestión técnica de la base de datos.

Busca una función de R que te permita identificar las líneas repetidas en un marco de datos.

Mediante la función duplicated podemos identificar las líneas repetidas. Para conocer más acerca de esta función he usado “help(”duplicated”)“. Basicamente esta función sirve para determinar qué elementos de un vector o de un data frame son elementos duplicados de elementos con subíndices más pequeños y devuelve un vector lógico que indica qué elementos (filas) son duplicados.

duplicated(genomas)

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Al ejecutar la función duplicated veo que en ningún caso hay duplicados.

Ordenar una tabla

La función order() produce los índices ordenados de los elementos de un vector. Consulta su ayuda y experimenta con ella. Vamos a usar esta función para ordenar las filas del marco de datos genomas por fecha, que es la tercera columna, llamada Assembly.Release.Date.

He ejecutado “help(”order”)” y he observado que esta función devuelve una permutación que reorganiza su primer argumento en orden ascendente o descendente, rompiendo los empates con argumentos adicionales.

# Observa que asignamos el resultado del ordenamiento al mismo objeto que
# estamos ordenando, para guardar el resultado sin tener que generar otro objeto.
genomas <- genomas[order(genomas$Assembly.Release.Date), ]
kable(genomas)

	Assembly.Accession	Organism.Name	Assembly.Release.Date	Assembly.Level	Assembly.Stats.ATGC.Count	Annotation.Count.Gene.Total	Assembly.Type	Assembly.Status
9	GCF_000209065.1	Trypanosoma cruzi	2005-08-02	Scaffold	89612356	25210	unresolved-diploid	current
11	GCA_000209065.1	Trypanosoma cruzi	2005-08-02	Scaffold	89612356	25210	unresolved-diploid	current
1	GCF_000002445.2	Trypanosoma brucei brucei TREU927	2005-12-14	Chromosome	27984184	10221	haploid	current
2	GCA_000002445.1	Trypanosoma brucei brucei TREU927	2005-12-14	Chromosome	27984184	10234	haploid	current
4	GCF_000210295.1	Trypanosoma brucei gambiense DAL972	2009-10-14	Chromosome	22110721	9930	haploid	current
5	GCA_000210295.1	Trypanosoma brucei gambiense DAL972	2009-10-14	Chromosome	22110721	9961	haploid	current
35	GCA_000227375.1	Trypanosoma vivax Y486	2011-08-11	Contig	24779298	4122	haploid	current
6	GCA_000227395.2	Trypanosoma congolense IL3000	2011-08-12	Scaffold	18802389	6348	haploid	current
21	GCA_000300495.1	Trypanosoma cruzi marinkellei	2012-10-01	Contig	34226154	10127	haploid	current
10	GCA_000188675.2	Trypanosoma cruzi	2012-10-02	Contig	38589511	10876	haploid	current
30	GCA_000492115.1	Trypanosoma rangeli SC58	2013-10-30	Contig	15111138	7475	haploid	current
20	GCA_000496795.1	Trypanosoma cruzi Dm28c	2013-11-13	Contig	27304309	11398	haploid	current
24	GCF_000691245.1	Trypanosoma grayi	2014-05-15	Scaffold	20791447	10676	haploid	current
25	GCA_000691245.1	Trypanosoma grayi	2014-05-15	Scaffold	20791447	10676	haploid	current
22	GCF_001457755.1	Trypanosoma equiperdum	2016-12-15	Scaffold	26113307	9129	haploid	current
23	GCA_001457755.2	Trypanosoma equiperdum	2016-12-15	Scaffold	26113307	9135	haploid	current
31	GCF_002087225.1	Trypanosoma theileri	2017-04-12	Scaffold	25688935	11312	haploid	current
32	GCA_002087225.1	Trypanosoma theileri	2017-04-12	Scaffold	25688935	11312	haploid	current
19	GCA_002219105.2	Trypanosoma cruzi cruzi	2017-09-15	Contig	50928259	22870	haploid	current
12	GCA_003177095.1	Trypanosoma cruzi	2018-05-30	Contig	87058484	29302	unresolved-diploid	current
13	GCA_003177105.1	Trypanosoma cruzi	2018-05-30	Contig	53271887	19115	haploid	current
3	GCA_003543875.1	Trypanosoma brucei equiperdum	2018-09-11	Chromosome	26988976	10293	haploid	current
14	GCA_003719155.1	Trypanosoma cruzi	2018-11-05	Scaffold	20565523	34344	diploid	current
7	GCF_003719485.1	Trypanosoma conorhini	2018-11-08	Scaffold	21137671	10321	haploid	current
8	GCA_003719485.1	Trypanosoma conorhini	2018-11-08	Scaffold	21137671	10321	haploid	current
15	GCA_003719455.1	Trypanosoma cruzi	2018-11-08	Scaffold	23852162	13503	haploid	current
28	GCF_003719475.1	Trypanosoma rangeli	2018-11-08	Scaffold	20977846	10589	haploid	current
29	GCA_003719475.1	Trypanosoma rangeli	2018-11-08	Scaffold	20977846	10589	haploid	current
16	GCA_013358655.1	Trypanosoma cruzi	2020-06-18	Scaffold	40788336	14352	haploid	current
17	GCA_015033625.1	Trypanosoma cruzi	2020-11-03	Chromosome	45527284	18708	haploid	current
18	GCA_015033655.1	Trypanosoma cruzi	2020-11-03	Chromosome	47195002	17651	haploid	current
34	GCA_019425475.1	Trypanosoma vivax	2021-07-30	Contig	23865041	10424	haploid	current
33	GCA_021307395.1	Trypanosoma vivax	2021-12-28	Contig	67823889	17806	haploid	current
26	GCF_022059095.1	Trypanosoma melophagium	2022-02-08	Contig	23303718	9756	haploid	current
27	GCA_022059095.1	Trypanosoma melophagium	2022-02-08	Contig	23303718	9768	haploid	current

Casualmente, la fecha estaba en formato “%Y-%m-%d”, en el que el orden alfabético coincide con el orden temporal. Al crear genomas, read.table() había interpretado la tercera columna como texto. Para que las fechas sean entendidas como tales, podemos usar la función as.Date(): Si le decimos que es una fecha con un formato concreto, ya podemos incluso podemos operar sumas y restas de fechas (nos dirá las diferencias entre ambas fechas por ejemplo).

genomas$Assembly.Release.Date <- as.Date(genomas$Assembly.Release.Date,
                                         '%Y-%m-%d')

Una vez convertida en “fecha”, la columna “Assembly.Release.Date” puede usarse correctamente, por ejemplo, como eje de abscisas en una gráfica:

# Aquí usamos "cumsum()" para obtener la suma acumulada del número de bases
# secuenciadas.
plot(x = genomas$Assembly.Release.Date,
     y = cumsum(genomas$Assembly.Stats.ATGC.Count),
     type = 'b',
     xlab = 'Fecha de publicación',
     ylab = 'Longitud total secuenciada (bp)',
     main = 'Genomas anotados de Trypanosoma')

Archivos JSON y GFF

Habiendo visto las opciones, decidimos descargar el genoma y la anotación funcional de una de ellas, por ejemplo GCA_022059095.1. Como tengo descargado el programa dataset dentro de la carpeta bin he tenido que indicar el path dentro del siguiente código.

if [ ! -e ncbi_data/data/GCA_022059095.1/genomic.gff ]; then
    if [ ! -e ncbi_dataset.zip ]; then
        /home/alicia/bin/datasets download genome accession GCA_022059095.1 --include genome,gff3
    fi
    unzip -o ncbi_dataset.zip
fi

Archive:  ncbi_dataset.zip
  inflating: README.md               
  inflating: ncbi_dataset/data/assembly_data_report.jsonl  
  inflating: ncbi_dataset/data/GCA_022059095.1/GCA_022059095.1_T.mel.1_genomic.fna  
  inflating: ncbi_dataset/data/GCA_022059095.1/genomic.gff  
  inflating: ncbi_dataset/data/dataset_catalog.json  
  inflating: md5sum.txt

En la carpeta ncbi_dataset/data encontrarás un archivo en formato JSON con información sobre los datos descargados, y un archivo en formato “jsonlines”. El formato “jsonlines” es una versión de JSON en la que cada registro ocupa una sola línea. Esto dificulta la visualización, pero agiliza el análisis de los datos en la línea de órdenes. Para ver mejor el contenido del archivo con la extensión jsonl, podríamos utilizar el programa jq, o bien utilizar la función prettify() del paquete jsonlite de R:

Por el mismo motivo que antes, me daba un error con la dirección que habia propuesta, utilizando en la consola “file.choose()” he podido buscar la ubicación del documento assembly_data_report.jsonl y la he copiado en el código de bajo tras la función readLines.

library('jsonlite')
prettify(readLines("/home/alicia/cuaderno/data/ncbi_dataset/data/assembly_data_report.jsonl"))

{
    "assemblyInfo": {
        "assemblyLevel": "Contig",
        "assemblyName": "T.mel.1",
        "assemblyType": "haploid",
        "submitter": "University of Edinburgh",
        "refseqCategory": "reference genome",
        "bioprojectLineage": [
            {
                "bioprojects": [
                    {
                        "accession": "PRJNA786535",
                        "title": "Trypanosoma melophagium isolate:St. Kilda Genome sequencing and assembly"
                    }
                ]
            }
        ],
        "sequencingTech": "Oxford Nanopore MinION; DNBseq",
        "blastUrl": "https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastSearch&PROG_DEF=blastn&BLAST_SPEC=GDH_GCA_022059095.1",
        "biosample": {
            "accession": "SAMN23701209",
            "lastUpdated": "2021-12-29T05:45:01.878",
            "publicationDate": "2021-12-29T00:00:00.000",
            "submissionDate": "2021-12-06T10:53:03.556",
            "sampleIds": [
                {
                    "label": "Sample name",
                    "value": "T. melophagium"
                },
                {
                    "db": "SRA",
                    "value": "SRS11278131"
                }
            ],
            "description": {
                "title": "MIGS Eukaryotic sample from Trypanosoma melophagium",
                "organism": {
                    "taxId": 715481,
                    "organismName": "Trypanosoma melophagium"
                },
                "comment": "Keywords: GSC:MIxS;MIGS:6.0"
            },
            "owner": {
                "name": "University of Edinburgh",
                "contacts": [
                    {

                    }
                ]
            },
            "models": [
                "MIGS.eu",
                "MIGS/MIMS/MIMARKS.host-associated"
            ],
            "bioprojects": [
                {
                    "accession": "PRJNA786535"
                }
            ],
            "package": "MIGS.eu.host-associated.6.0",
            "attributes": [
                {
                    "name": "isolate",
                    "value": "St. Kilda"
                },
                {
                    "name": "isolation_source",
                    "value": "blood"
                },
                {
                    "name": "ecotype",
                    "value": "St. Kilda"
                },
                {
                    "name": "collection_date",
                    "value": "2020-09-17"
                },
                {
                    "name": "env_broad_scale",
                    "value": "Cultured"
                },
                {
                    "name": "env_local_scale",
                    "value": "Cultured"
                },
                {
                    "name": "env_medium",
                    "value": "Cultured"
                },
                {
                    "name": "geo_loc_name",
                    "value": "United Kingdom: Scotland, St. Kilda"
                },
                {
                    "name": "host",
                    "value": "Ovis aries"
                },
                {
                    "name": "isol_growth_condt",
                    "value": "This study"
                },
                {
                    "name": "lat_lon",
                    "value": "57.8287 N 8.6352 W"
                }
            ],
            "status": {
                "status": "live",
                "when": "2021-12-29T05:45:01.878"
            },
            "collectionDate": "2020-09-17",
            "ecotype": "St. Kilda",
            "geoLocName": "United Kingdom: Scotland, St. Kilda",
            "host": "Ovis aries",
            "isolate": "St. Kilda",
            "isolationSource": "blood",
            "latLon": "57.8287 N 8.6352 W"
        },
        "assemblyStatus": "current",
        "pairedAssembly": {
            "accession": "GCF_022059095.1",
            "status": "current",
            "annotationName": "Annotation submitted by University of Edinburgh"
        },
        "bioprojectAccession": "PRJNA786535",
        "assemblyMethod": "Wtdbg2 v. 13/05/2021",
        "releaseDate": "2022-02-08"
    },
    "assemblyStats": {
        "totalSequenceLength": "23303718",
        "totalUngappedLength": "23303718",
        "numberOfContigs": 64,
        "contigN50": 505851,
        "contigL50": 15,
        "numberOfScaffolds": 64,
        "scaffoldN50": 505851,
        "scaffoldL50": 15,
        "numberOfComponentSequences": 64,
        "gcCount": "9608043",
        "gcPercent": 41,
        "genomeCoverage": "225",
        "atgcCount": "23303718"
    },
    "annotationInfo": {
        "name": "Annotation submitted by University of Edinburgh",
        "provider": "University of Edinburgh",
        "releaseDate": "2022-02-08",
        "stats": {
            "geneCounts": {
                "total": 9768,
                "proteinCoding": 9768
            }
        }
    },
    "wgsInfo": {
        "wgsProjectAccession": "JAJSEB01",
        "masterWgsUrl": "https://www.ncbi.nlm.nih.gov/nuccore/JAJSEB000000000.1",
        "wgsContigsUrl": "https://www.ncbi.nlm.nih.gov/Traces/wgs/JAJSEB01"
    },
    "currentAccession": "GCA_022059095.1",
    "accession": "GCA_022059095.1",
    "pairedAccession": "GCF_022059095.1",
    "sourceDatabase": "SOURCE_DATABASE_GENBANK",
    "organism": {
        "taxId": 715481,
        "organismName": "Trypanosoma melophagium",
        "infraspecificNames": {
            "ecotype": "St. Kilda",
            "isolate": "St. Kilda"
        }
    }
}

Volviendo a lo que nos interesa, vamos a cargar en un marco de datos de una sesión de R la anotación funcional del genoma escogido. Aunque existen funciones específicas para leer archivos GFF, usaremos la función básica read.table(), con alguna precaución. Ocurre lo mismo, pongo la ubicación del archivo con file.choose().

# Como no tiene cabecera, le damos los nombres de las columnas:
gff <- read.table("/home/alicia/cuaderno/data/ncbi_dataset/data/GCA_022059095.1/genomic.gff", sep = '\t',
                  col.names = c('seqid', 'source', 'type', 'start', 'end',
                                'score', 'strand', 'phase', 'attributes'),
                  quote = '')

Al leer el archivo GFF con read.table(), debemos especificar que el separador de columnas es el tabulador y que no hay texto entrecomillado en el archivo. Esto último se consigue con la opción quote = '', y sirve para evitar que todo lo que se encuentre entre dos apóstrofes (’) se considere parte de un mismo campo de información, aunque sean cientos de líneas.

Para comprobar que el archivo ha sido leído completo y sin errores, conviene revisar las dimensiones del marco de datos (dim(gff)) o incluso visualizar el principio (head(gff)) y el final (tail(gff)) del objeto, que es demasiado grande para mostrarlo en pantalla completo.

dim(gff)

[1] 42497     9

head(gff)

              seqid  source   type start    end score strand phase
1 JAJSEB010000004.1 Genbank region     1 117947     .      +     .
2 JAJSEB010000004.1 Genbank   gene   246    935     .      -     .
3 JAJSEB010000004.1 Genbank   mRNA   246    935     .      -     .
4 JAJSEB010000004.1 Genbank   exon   801    935     .      -     .
5 JAJSEB010000004.1 Genbank   exon   246    614     .      -     .
6 JAJSEB010000004.1 Genbank    CDS   801    935     .      -     0
                                                                                                                                                                                                                                                                           attributes
1 ID=JAJSEB010000004.1:1..117947;Dbxref=taxon:715481;collection-date=2020-09-17;country=United Kingdom: Scotland%2C St. Kilda;ecotype=St. Kilda;gbkey=Src;genome=genomic;isolate=St. Kilda;isolation-source=blood;lat-lon=57.8287 N 8.6352 W;mol_type=genomic DNA;nat-host=Ovis aries
2                                                                                                                                                                                ID=gene-LSM04_009031;Name=LSM04_009031;gbkey=Gene;gene_biotype=protein_coding;locus_tag=LSM04_009031
3                                                    ID=rna-gnl|WGS:JAJSEB|jg9720.t1;Parent=gene-LSM04_009031;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
4                                      ID=exon-gnl|WGS:JAJSEB|jg9720.t1-1;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
5                                      ID=exon-gnl|WGS:JAJSEB|jg9720.t1-2;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
6                              ID=cds-KAH9601282.1;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;Dbxref=NCBI_GP:KAH9601282.1;Name=KAH9601282.1;gbkey=CDS;locus_tag=LSM04_009031;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase;protein_id=KAH9601282.1

tail(gff)

                  seqid  source type start   end score strand phase
42492 JAJSEB010000064.1 Genbank exon 61805 62122     .      -     .
42493 JAJSEB010000064.1 Genbank  CDS 61805 62122     .      -     0
42494 JAJSEB010000064.1 Genbank gene 67444 67818     .      +     .
42495 JAJSEB010000064.1 Genbank mRNA 67444 67818     .      +     .
42496 JAJSEB010000064.1 Genbank exon 67444 67818     .      +     .
42497 JAJSEB010000064.1 Genbank  CDS 67444 67818     .      +     0
                                                                                                                                                                                                                                 attributes
42492           ID=exon-gnl|WGS:JAJSEB|jg153.t1-1;Parent=rna-gnl|WGS:JAJSEB|jg153.t1;gbkey=mRNA;locus_tag=LSM04_007206;orig_protein_id=gnl|WGS:JAJSEB|jg153.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg153.t1;product=hypothetical protein
42493 ID=cds-KAH9577183.1;Parent=rna-gnl|WGS:JAJSEB|jg153.t1;Dbxref=NCBI_GP:KAH9577183.1;Name=KAH9577183.1;gbkey=CDS;locus_tag=LSM04_007206;orig_transcript_id=gnl|WGS:JAJSEB|jg153.t1;product=hypothetical protein;protein_id=KAH9577183.1
42494                                                                                                                                  ID=gene-LSM04_000875;Name=LSM04_000875;gbkey=Gene;gene_biotype=protein_coding;locus_tag=LSM04_000875
42495                        ID=rna-gnl|WGS:JAJSEB|jg154.t1;Parent=gene-LSM04_000875;gbkey=mRNA;locus_tag=LSM04_000875;orig_protein_id=gnl|WGS:JAJSEB|jg154.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg154.t1;product=hypothetical protein
42496           ID=exon-gnl|WGS:JAJSEB|jg154.t1-1;Parent=rna-gnl|WGS:JAJSEB|jg154.t1;gbkey=mRNA;locus_tag=LSM04_000875;orig_protein_id=gnl|WGS:JAJSEB|jg154.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg154.t1;product=hypothetical protein
42497 ID=cds-KAH9577164.1;Parent=rna-gnl|WGS:JAJSEB|jg154.t1;Dbxref=NCBI_GP:KAH9577164.1;Name=KAH9577164.1;gbkey=CDS;locus_tag=LSM04_000875;orig_transcript_id=gnl|WGS:JAJSEB|jg154.t1;product=hypothetical protein;protein_id=KAH9577164.1

Seleccionar filas y columnas

En el objeto gff hay 5 valores diferentes de la variable type: “CDS”, “exon”, “gene”, “mRNA” y “region”. Puedes comprobarlo con unique(gff$type), o con table(gff$type). Podemos separar la tabla en cinco marcos de datos, que correspondan a esos cinco tipos de anotación funcional.

unique(gff$type)

[1] "region" "gene"   "mRNA"   "exon"   "CDS"

table(gff$type)


   CDS   exon   gene   mRNA region 
 11304  11304   9768  10057     64

# Usamos un vector booleano como índice de las filas.
# En las tablas nuevas, no necesitamos la tercera columna
# y la eliminamos con el "-3".

gene   <- gff[gff$type == 'gene',   -3]
mrna   <- gff[gff$type == 'mRNA',   -3]
exon   <- gff[gff$type == 'exon',   -3]
cds    <- gff[gff$type == 'CDS',    -3]
region <- gff[gff$type == 'region', -3]

# Alinear las órdenes nos ayuda a detectar posibles errores

Desafío

Identifica en cada uno de estos marcos de datos qué columnas son invariables. Es decir, las columnas que tienen el mismo valor en todas las líneas. Y elimina esas columnas. Visualmente en todos vemos que se debería eliminar la columna de source, score y phase, ya que en source es GenBank para todos, y en score y phase solo hay puntos. Sin embargo, las tablas tienen muchas filas, así que conviene usar comandos para averiguar esto de forma más sencilla. Podemos usar los siguientes comandos:

unique(gene$score)

[1] "."

unique(gene$phase)

[1] "."

table(gene$score)


   . 
9768

Para eliminarlo podríamos haberlo hecho a través de asignación de NULL:

gene$source <- NULL
gene$score  <- NULL
gene$phase  <- NULL

Agrupar, operar y resumir

Habiendo separado los diferentes tipos de anotaciones, ahora resulta fácil obtener estadísticas de cada uno. Podemos querer una tabla que nos diga para cada tipo de anotación, cuántos elementos de ese tipo hay en el genoma, qué longitud media tienen y qué desviación estándar. Hay diferentes maneras de conseguirlo. La más rudimentaria implicaría obtener cada estadística individualmente y después juntarlas en una nueva tabla. Por ejemplo, para regiones codificantes (“CDS”):

cds$longitud <- cds$end - cds$start + 1
resumen <- data.frame(
   tipo       = 'CDS',
   número     = nrow(cds),
   long.Media = mean(cds$longitud),
   long.SD    = sd(cds$longitud)
)

Desafío

¿Sabrías completar el marco de datos resumen con los datos de los otros cuatro tipos de anotación funcional? Primero calculamos la longitud de cada dato:

gene$longitud   <- gene$end - gene$start + 1
mrna$longitud   <- mrna$end - mrna$start + 1
exon$longitud   <- exon$end - exon$start + 1
region$longitud <- region$end - region$start + 1

Ahora lo juntamos todo en una tabla:

resumen_completo <- rbind(
  resumen, 
  data.frame(tipo = 'gene',   número = nrow(gene),   long.Media = mean(gene$longitud),   long.SD = sd(gene$longitud)),
  data.frame(tipo = 'mRNA',   número = nrow(mrna),   long.Media = mean(mrna$longitud),   long.SD = sd(mrna$longitud)),
  data.frame(tipo = 'exon',   número = nrow(exon),   long.Media = mean(exon$longitud),   long.SD = sd(exon$longitud)),
  data.frame(tipo = 'region', número = nrow(region), long.Media = mean(region$longitud), long.SD = sd(region$longitud))
)


print(resumen_completo)

    tipo número long.Media    long.SD
1    CDS  11304   1358.661   1377.899
2   gene   9768   1553.822   1466.322
3   mRNA  10057   1549.331   1466.240
4   exon  11304   1358.661   1377.899
5 region     64 364120.594 314801.149

Este tipo de operación es bastante habitual, y se puede hacer sin necesidad de crear marcos de datos parciales. La función split() divide las líneas de un marco de datos en tantos marcos de datos como valores diferentes (niveles) haya en el factor que usemos para definir los grupos. El resultado es un solo objeto: una lista de marcos de datos:

# Esto es una lista de marcos de datos:
porTipo <- split(gff, gff$type)

# Estas funciones nos permiten inspeccionar algunas características
# del nuevo objeto:
class(porTipo)

[1] "list"

length(porTipo)

[1] 5

names(porTipo)

[1] "CDS"    "exon"   "gene"   "mRNA"   "region"

sapply(porTipo, class)

         CDS         exon         gene         mRNA       region 
"data.frame" "data.frame" "data.frame" "data.frame" "data.frame"

sapply(porTipo, nrow)

   CDS   exon   gene   mRNA region 
 11304  11304   9768  10057     64

La función sapply(), toma al menos dos argumentos: una lista (o vector) y una función. Devuelve un vector (si es posible) con los resultados de aplicar la función a cada elemento de la lista. Con esto, el resumen que queríamos se puede obtener así:

# Primero calculamos las longitudes de todas las anotaciones en gff:
gff$longitud <- gff$end - gff$start + 1

# Volvemos a separar las filas por tipo, ahora que tenemos las longitudes:
porTipo <- split(gff, gff$type)

resumen <- data.frame(
   tipo = names(porTipo),
   número = sapply(porTipo, nrow),
   long.Media = sapply(porTipo, function(x) mean(x$longitud)),
   long.SD    = sapply(porTipo, function(x)   sd(x$longitud))
)

# Aprovechamos para mostrar nombres de columnas más legibles:
kable(resumen,
      col.names = c('Tipo', 'Número', 'Longitud Media',
                    'Desviación estándar'),
      row.names = FALSE)

Tipo	Número	Longitud Media	Desviación estándar
CDS	11304	1358.661	1377.899
exon	11304	1358.661	1377.899
gene	9768	1553.822	1466.322
mRNA	10057	1549.331	1466.240
region	64	364120.594	314801.149

   kable(resumen)

	tipo	número	long.Media	long.SD
CDS	CDS	11304	1358.661	1377.899
exon	exon	11304	1358.661	1377.899
gene	gene	9768	1553.822	1466.322
mRNA	mRNA	10057	1549.331	1466.240
region	region	64	364120.594	314801.149

Las operaciones de agrupar y resumir se pueden hacer también con funciones del paquete dplyr, que forma parte de un conjunto de paquetes diseñados para funcionar armoniosamente entre ellos. Son los paquetes del tidyverse, muy populares en ciencia de datos. Casi se pueden considerar un dialecto de R. Las funciones (o verbos) más interesantes de dplyr son group_by() y summarise().

library('dplyr')


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

# Observa el uso de "tuberías" ("%>%"):
resumen <- gff %>% group_by(type) %>% summarise(
   número     = n(),
   Long.Media = mean(longitud),
   Long.SD    = sd(longitud)
)
kable(resumen,
      col.names = c('Tipo', 'Número', 'Longitud Media',
                    'Desviación estándar'),
      row.names = FALSE) %>%
   kable(full_width = FALSE)

x
\|Tipo \| Número\| Longitud Media\| Desviación estándar\|
\|:——\|——:\|————–:\|——————-:\|
\|CDS \| 11304\| 1358.661\| 1377.899\|
\|exon \| 11304\| 1358.661\| 1377.899\|
\|gene \| 9768\| 1553.822\| 1466.322\|
\|mRNA \| 10057\| 1549.331\| 1466.240\|
\|region \| 64\| 364120.594\| 314801.149\|

Desafío

Crea una tabla con el número total de genes en cada cromosoma. ¿Puedes contar cuántos hay codificados en cada sentido (“strand”) y en cada cromosoma?

library(dplyr)
genes_por_cromosomaysentido <- gff %>%
  dplyr::filter(type == 'gene') %>%      
  group_by(seqid, strand) %>%             
  summarise(total_genes = n())

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by seqid and strand.
ℹ Output is grouped by seqid.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(seqid, strand))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

print(genes_por_cromosomaysentido)

# A tibble: 116 × 3
# Groups:   seqid [64]
   seqid             strand total_genes
   <chr>             <chr>        <int>
 1 JAJSEB010000001.1 +              376
 2 JAJSEB010000001.1 -              201
 3 JAJSEB010000002.1 -                3
 4 JAJSEB010000003.1 -               39
 5 JAJSEB010000004.1 +               34
 6 JAJSEB010000004.1 -               14
 7 JAJSEB010000005.1 +               32
 8 JAJSEB010000005.1 -              114
 9 JAJSEB010000006.1 +              208
10 JAJSEB010000006.1 -               10
# ℹ 106 more rows

En el código filtramos por genes y agrupamos por sequid, que es como se denominan los diferentes cromosomas en el objeto gff y por strand, que es el sentido. El resultado nos da el numero de genes por cromosoma y por sentido. ## Comparar o cruzar datos

Imagina que un estudio señala los cromosomas siguientes como posiblemente asociados a la resistencia a un fármaco. Quieres saber qué genes hay en estos cinco cromosomas:

JAJSEB010000048.1
JAJSEB010000012.1
JAJSEB010000020.1
JAJSEB010000016.1
JAJSEB010000033.1

Antes hemos usado el operador lógico == para crear un vector booleano con el que seleccionar algunas líneas. En este caso usaremos el operador %in%, que compara conjuntos o vectores de caracteres:

candidatos <- c("JAJSEB010000048.1",
                "JAJSEB010000012.1",
                "JAJSEB010000020.1",
                "JAJSEB010000016.1",
                "JAJSEB010000033.1")
interesantes <- gff[gff$seqid %in% candidatos, ]
dim(interesantes)

[1] 4792   10

Desafío

¿Puedes limitar la selección a anotaciones de tipo “gene”? Utilizamos el paquete dplyr como antes:

library(dplyr)

interesantes2 <- gff %>%
  dplyr::filter(seqid %in% candidatos, type == "gene")

# Comprobamos las nuevas dimensiones
dim(interesantes2)

[1] 1097   10

Ahora el objeto interesantes2 es más pequeño, porque se ha limitado la selección a anotaciones de tipo gene.

El estudio avanza y reduce el tamaño de las regiones candidatas, así que ahora necesitas identificar los genes presentes en exactamente las regiones siguientes:

JAJSEB010000048.1:184638-253478
JAJSEB010000012.1:133444-348860
JAJSEB010000020.1:73393-149614
JAJSEB010000016.1:161247-272009
JAJSEB010000033.1:123135-170687

Aunque existen programas específicos para comparar segmentos cromosómicos (bedtools), continuaremos practicando con marcos de datos de R. Ahora ya no podemos usar %in%, porque de cada cromosoma queremos unos segmentos diferentes. Una estrategia es crear cinco filtros booleanos, y después combinarlos con la operación lógica “o”, es decir |. En cada filtro seleccionaremos las anotaciones que estén completamente contenidas en el intervalo indicado.

f1 <- gff$seqid == 'JAJSEB010000048.1' & gff$start >= 184638 & gff$end <= 253478
f2 <- gff$seqid == 'JAJSEB010000012.1' & gff$start >= 133444 & gff$end <= 348860
f3 <- gff$seqid == 'JAJSEB010000020.1' & gff$start >= 73393  & gff$end <= 149614
f4 <- gff$seqid == 'JAJSEB010000016.1' & gff$start >= 161247 & gff$end <= 272009
f5 <- gff$seqid == 'JAJSEB010000033.1' & gff$start >= 123135 & gff$end <= 170687

f <- f1 | f2 | f3 | f4 | f5

interesantes <- gff[f, ]

Desafío

¿Puedes cambiar los filtros para que incluyan también genes parcialmente contenidos en los intervalos indicados?

f1 <- gff$seqid == 'JAJSEB010000048.1' & gff$start <= 253478 & gff$end >= 184638
f2 <- gff$seqid == 'JAJSEB010000012.1' & gff$start <= 348860 & gff$end >= 133444
f3 <- gff$seqid == 'JAJSEB010000020.1' & gff$start <= 149614 & gff$end >= 73393
f4 <- gff$seqid == 'JAJSEB010000016.1' & gff$start <= 272009 & gff$end >= 161247
f5 <- gff$seqid == 'JAJSEB010000033.1' & gff$start <= 170687 & gff$end >= 123135

f <- f1 | f2 | f3 | f4 | f5

interesantes <- gff[f, ]

Se intercambia el inicio y el final del intervalo. El final del gen tiene que ser igual o mayor al inicio del intervalo y el inicio del gen tiene que ser menor o igual que el final del intervalo. ## Las llaves y la función match()

En las bases de datos relacionales se llama llave a la columna o al conjunto de columnas cuyos valores identifican unívocamente cada línea de la tabla. En los marcos de datos de R, las líneas pueden tener nombre, el cual puede usarse como llave o índice para seleccionar líneas. Pero el uso del nombre de las líneas no es recomendable. Lo más práctico es que en cualquier tabla haya una columna cuyos valores no puedan repetirse. Esto ayuda a evitar duplicidades, mantener la integridad de la tabla y relacionar unas tablas con otras.

¿Cuál es la llave en un archivo GFF? Es la columna de atributos, se puede ver utilizando el comando head(gff)

head(gff)

              seqid  source   type start    end score strand phase
1 JAJSEB010000004.1 Genbank region     1 117947     .      +     .
2 JAJSEB010000004.1 Genbank   gene   246    935     .      -     .
3 JAJSEB010000004.1 Genbank   mRNA   246    935     .      -     .
4 JAJSEB010000004.1 Genbank   exon   801    935     .      -     .
5 JAJSEB010000004.1 Genbank   exon   246    614     .      -     .
6 JAJSEB010000004.1 Genbank    CDS   801    935     .      -     0
                                                                                                                                                                                                                                                                           attributes
1 ID=JAJSEB010000004.1:1..117947;Dbxref=taxon:715481;collection-date=2020-09-17;country=United Kingdom: Scotland%2C St. Kilda;ecotype=St. Kilda;gbkey=Src;genome=genomic;isolate=St. Kilda;isolation-source=blood;lat-lon=57.8287 N 8.6352 W;mol_type=genomic DNA;nat-host=Ovis aries
2                                                                                                                                                                                ID=gene-LSM04_009031;Name=LSM04_009031;gbkey=Gene;gene_biotype=protein_coding;locus_tag=LSM04_009031
3                                                    ID=rna-gnl|WGS:JAJSEB|jg9720.t1;Parent=gene-LSM04_009031;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
4                                      ID=exon-gnl|WGS:JAJSEB|jg9720.t1-1;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
5                                      ID=exon-gnl|WGS:JAJSEB|jg9720.t1-2;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;gbkey=mRNA;locus_tag=LSM04_009031;orig_protein_id=gnl|WGS:JAJSEB|jg9720.t1.CDS1;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase
6                              ID=cds-KAH9601282.1;Parent=rna-gnl|WGS:JAJSEB|jg9720.t1;Dbxref=NCBI_GP:KAH9601282.1;Name=KAH9601282.1;gbkey=CDS;locus_tag=LSM04_009031;orig_transcript_id=gnl|WGS:JAJSEB|jg9720.t1;product=Metal-independent alpha-mannosidase;protein_id=KAH9601282.1
  longitud
1   117947
2      690
3      690
4      135
5      369
6      135

El código siguiente descarga la tabla de genomas de Trypanosoma ensamblados a nivel de cromosoma, de la European Nucleotide Archive

# Aprenderemos esta sintaxis en el próximo tema:
portal <- 'https://www.ebi.ac.uk/ena/portal/api/search?'
result <- 'result=assembly&'
query  <- 'query=tax_tree(5690)&'
fields <- 'fields=accession%2Cversion%2Cscientific_name%2Cbase_count%2Clast_updated%2Cassembly_level'
URL <- paste0(portal, result, query, fields)
genomas.ena <- read.table(URL, header = TRUE, sep = '\t', quote = '', na.strings = '')

Tanto el nuevo marco de datos genomas.ena, como el que habíamos creado antes con datos de NCBI, genomas, tienen una columna de valores únicos: “accession” o “Assembly.Accession”. Para los registros comunes entre las dos tablas, podemos transferir información de una a otra. Por ejemplo, vamos a añadir a la tabla genomas la fecha de última modificación del ensamblaje, que encontramos en genomas.ena. Para asegurarnos de que a cada ensamblaje le asignamos el valor correcto, usaremos la función match().

En genomas, “Assembly.Accession” incluye un sufijo que indica el número de versión del ensamblaje. Este sufijo no está presente en genomas.ena$accession. Para hacer comparables ambas tablas, quitaremos el sufijo de la primera tabla.

# Este es un uso avanzado de `sapply()` y `strsplit()`:
genomas$acc <- sapply(
   strsplit(genomas$Assembly.Accession, '.', fixed = TRUE),
   '[', 1
)

genomas$last_updated <- genomas.ena[match(genomas$acc,
                                          genomas.ena$accession),
                                    'last_updated']

Desafíos

¿Cómo comprobarías que al eliminar el sufijo de la versión, genomas$acc sigue siendo una llave de la tabla genomas?
Siguiendo el mismo esquema añade a genomas el “base_count” y compáralo con el “Assembly.Stats.ATGC.Count”. Si observas alguna diferencia, intenta explicarla.

Apéndice

En los archivos GFF, la columna attributes es una lista de parejas de valores “clave = valor”, separadas por punto y coma (“;”). Habiendo creado antes un marco de datos para cada tipo de anotación (genes, CDS, etc.), es de suponer que todas las líneas de un marco de datos tendrán el mismo conjunto de parejas “clave = valor”. Aquí se muestra una manera de extraer esta información y convertirla en columnas adicionales de un marco de datos.

Veamos el primer valor de attributes en el marco de datos gene:

gene[1, 'attributes']

[1] "ID=gene-LSM04_009031;Name=LSM04_009031;gbkey=Gene;gene_biotype=protein_coding;locus_tag=LSM04_009031"

La función strsplit() en R divide una cadena de caracteres en partes delimitadas por algún carácter o expresión regular. En este caso, queremos dividir esa cadena por los puntos y coma:

strsplit(gene[1, 'attributes'], ';')

[[1]]
[1] "ID=gene-LSM04_009031"        "Name=LSM04_009031"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_009031"

strsplit() ha dividido la cadena de caracteres original en una lista que contiene un único vector de cinco cadenas de caracteres. El motivo de generar una lista en lugar de un simple vector es que strsplit() puede actuar simultáneamente sobre todas las cadenas de caracteres de un vector, y necesita repartir el resultado en diferentes vectores:

strsplit(gene[1:3, 'attributes'], ';')

[[1]]
[1] "ID=gene-LSM04_009031"        "Name=LSM04_009031"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_009031"     

[[2]]
[1] "ID=gene-LSM04_007573"        "Name=LSM04_007573"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_007573"     

[[3]]
[1] "ID=gene-LSM04_000706"        "Name=LSM04_000706"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_000706"

Podemos guardar la lista de los vectores de los atributos de cada gen en un nuevo objeto, para continuar operando sobre el después:

atributos <- strsplit(gene[ , 'attributes'], ';')

# Veamos los 6 primeros elementos:
head(atributos)

[[1]]
[1] "ID=gene-LSM04_009031"        "Name=LSM04_009031"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_009031"     

[[2]]
[1] "ID=gene-LSM04_007573"        "Name=LSM04_007573"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_007573"     

[[3]]
[1] "ID=gene-LSM04_000706"        "Name=LSM04_000706"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_000706"     

[[4]]
[1] "ID=gene-LSM04_000559"        "Name=LSM04_000559"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_000559"     

[[5]]
[1] "ID=gene-LSM04_006454"        "Name=LSM04_006454"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_006454"     

[[6]]
[1] "ID=gene-LSM04_009090"        "Name=LSM04_009090"          
[3] "gbkey=Gene"                  "gene_biotype=protein_coding"
[5] "locus_tag=LSM04_009090"

# Y la longitud total de la lista:
length(atributos)

[1] 9768

Una vez tenemos una lista de los atributos de cada gen, podemos usar funciones de la familia lapply() para operar en paralelo sobre todos sus elementos. Examina la ayuda de lapply() para familiarizarte con estas funciones. A continuación, la usaremos para separar la clave del valor en todos los atributos de todos los genes. Además asignaremos la clave como nombre del valor correspondiente (en R, los elementos de un vector pueden tener nombre).

atributos2 <- lapply(
   atributos,
   # Definimos una función para aplicarla a cada elemento de "atributos".
   function(x) {
      # "x" es un vector de atributos de un solo gen
      lista <- strsplit(x, '=')
      # "lista" es una lista de vectores de 2 componentes, clave y valor
      valores <- sapply(lista, '[', 2)
      names(valores) <- sapply(lista, '[', 1)
      return(valores)
   }
)
head(atributos2)

[[1]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_009031"      "LSM04_009031"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_009031" 

[[2]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_007573"      "LSM04_007573"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_007573" 

[[3]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_000706"      "LSM04_000706"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_000706" 

[[4]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_000559"      "LSM04_000559"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_000559" 

[[5]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_006454"      "LSM04_006454"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_006454" 

[[6]]
                 ID                Name               gbkey        gene_biotype 
"gene-LSM04_009090"      "LSM04_009090"              "Gene"    "protein_coding" 
          locus_tag 
     "LSM04_009090"

El orden de los atributos se ha mantenido idéntico al orden de sus correspondientes genes. Por eso podemos ahora añadir las columnas de atributos que queramos al marco de datos gene:

gene$ID           <- sapply(atributos2, '[', 'ID')
gene$Name         <- sapply(atributos2, '[', 'Name')
gene$gene_biotype <- sapply(atributos2, '[', 'gene_biotype')
gene$locus_tag    <- sapply(atributos2, '[', 'locus_tag')

f <- names(gene) != 'attributes'
kable(gene[1:15, f])

	seqid	start	end	strand	longitud	ID	Name	gene_biotype	locus_tag
2	JAJSEB010000004.1	246	935	-	690	gene-LSM04_009031	LSM04_009031	protein_coding	LSM04_009031
8	JAJSEB010000004.1	2204	2776	+	573	gene-LSM04_007573	LSM04_007573	protein_coding	LSM04_007573
12	JAJSEB010000004.1	3011	3706	+	696	gene-LSM04_000706	LSM04_000706	protein_coding	LSM04_000706
16	JAJSEB010000004.1	4768	4986	+	219	gene-LSM04_000559	LSM04_000559	protein_coding	LSM04_000559
20	JAJSEB010000004.1	5540	8959	+	3420	gene-LSM04_006454	LSM04_006454	protein_coding	LSM04_006454
24	JAJSEB010000004.1	9198	9512	+	315	gene-LSM04_009090	LSM04_009090	protein_coding	LSM04_009090
28	JAJSEB010000004.1	9877	12936	+	3060	gene-LSM04_002876	LSM04_002876	protein_coding	LSM04_002876
32	JAJSEB010000004.1	13160	14620	+	1461	gene-LSM04_001088	LSM04_001088	protein_coding	LSM04_001088
36	JAJSEB010000004.1	14840	16024	+	1185	gene-LSM04_006263	LSM04_006263	protein_coding	LSM04_006263
40	JAJSEB010000004.1	16367	16935	+	569	gene-LSM04_007732	LSM04_007732	protein_coding	LSM04_007732
46	JAJSEB010000004.1	17016	17411	+	396	gene-LSM04_006531	LSM04_006531	protein_coding	LSM04_006531
52	JAJSEB010000004.1	17834	19039	+	1206	gene-LSM04_007917	LSM04_007917	protein_coding	LSM04_007917
56	JAJSEB010000004.1	19420	20313	+	894	gene-LSM04_009571	LSM04_009571	protein_coding	LSM04_009571
60	JAJSEB010000004.1	20625	21875	+	1251	gene-LSM04_003035	LSM04_003035	protein_coding	LSM04_003035
64	JAJSEB010000004.1	22111	22929	+	819	gene-LSM04_003619	LSM04_003619	protein_coding	LSM04_003619