Dividir una sola columna de pares clave-valor en varias columnas

Dos formatos de archivo ampliamente utilizados en bioinformática, VCF y GTF, tienen columnas únicas que contienen información de anotaciones. Esto hace que sea un poco inconveniente trabajar con ellos en R cuando se usan marcos de datos porque los valores deben desempaquetarse, es decir, dividirse. Además, esto viola una de las condiciones de los datos ordenados, que es que cada celda es un valor único. En esta publicación, utilizaremos herramientas de la tidyverse para dividir los valores en varias columnas para facilitar el trabajo con los datos en R.

Para comenzar, instale el tidyverse si aún no lo has hecho.

if(!require("tidyverse"))
library(tidyverse)

Tengo un pequeño paquete llamado importbio que se puede usar para cargar un archivo VCF y GTF. Puedes instalarlo usando el remotes paquete.

if(!require("remotes"))

if(!require("importbio"))

Dividir la columna de información en un archivo VCF

Cargaremos un pequeño archivo VCF usando importvcf.

library(importbio)
my_vcf <- importvcf("https://raw.githubusercontent.com/davetang/learning_vcf_file/master/eg/Pfeiffer.vcf")

my_vcf %>%
  select(info) %>%
  head()
# A tibble: 6 × 1
  info                                                                                                         
  <chr>                                                                                                        
1 AC=2;AF=1.00;AN=2;DB;DP=11;FS=0.000;HRun=0;HaplotypeScore=41.3338;MQ0=0;MQ=61.94;QD=23.51;set=variant        
2 AC=1;AF=0.50;AN=2;BaseQRankSum=1.455;DB;DP=21;Dels=0.00;FS=1.984;HRun=0;HaplotypeScore=0.0000;MQ0=0;MQ=60.00…
3 AC=1;AF=0.50;AN=2;BaseQRankSum=1.934;DP=48;Dels=0.00;FS=4.452;HRun=0;HaplotypeScore=0.5784;MQ0=0;MQ=59.13;MQ…
4 AC=1;AF=0.50;AN=2;BaseQRankSum=-4.517;DB;DP=29;Dels=0.00;FS=1.485;HRun=0;HaplotypeScore=0.0000;MQ0=0;MQ=56.9…
5 AC=1;AF=0.50;AN=2;BaseQRankSum=0.199;DB;DP=33;Dels=0.00;FS=0.000;HRun=1;HaplotypeScore=1.8893;MQ0=0;MQ=60.00…
6 AC=1;AF=0.50;AN=2;BaseQRankSum=-0.259;DB;DP=12;Dels=0.00;FS=0.000;HRun=1;HaplotypeScore=0.0000;MQ0=0;MQ=53.2…

Tenga en cuenta que el info columna está repleta de todo tipo de información para cada variante. También tenga en cuenta el formato consistente de la info columna: cada anotación está separada por un punto y coma (;) y las anotaciones se almacenan como pares clave-valor con un signo igual en el medio.

En primer lugar, utilizaremos separate_rows para crear una nueva fila para cada anotación usando ; como separador/delimitador; tenga en cuenta que he incluido \s* después ;, que es una expresión regular para especificar un único espacio en blanco que aparece 0 o más veces. Al incluir la expresión regular, espacios en blanco después ; será eliminado, lo cual es bueno porque no queremos espacios en blanco en nuestros datos. Además mutate la llamada se usa antes de llamar separate_rows y se utiliza para eliminar un punto y coma final, si existe.

my_vcf %>%
  mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) %>%
  separate_rows(info, sep = ";\s*") %>%
  head()
# A tibble: 6 × 10
  vid              chrom    pos id         ref   alt   qual   filter info     type 
  <chr>            <fct>  <int> <chr>      <chr> <chr> <chr>  <chr>  <chr>    <chr>
1 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AC=2     ins  
2 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AF=1.00  ins  
3 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AN=2     ins  
4 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   DB       ins  
5 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   DP=11    ins  
6 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   FS=0.000 ins  

El siguiente paso es dividir los pares clave-valor y usaremos el separate función para separar los pares en dos columnas, a las que llamaremos key y value, utilizando el signo igual como separador/delimitador. A veces, a una clave le falta un valor y, en estos casos, el valor será NA.

my_vcf %>%
  mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) %>%
  separate_rows(info, sep = ";\s*") %>%
  separate(info, c('key', 'value'), sep = "=") %>%
  head(10)
# A tibble: 10 × 11
   vid              chrom    pos id         ref   alt   qual   filter key            value   type 
   <chr>            <fct>  <int> <chr>      <chr> <chr> <chr>  <chr>  <chr>          <chr>   <chr>
 1 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AC             2       ins  
 2 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AF             1.00    ins  
 3 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   AN             2       ins  
 4 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   DB             NA      ins  
 5 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   DP             11      ins  
 6 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   FS             0.000   ins  
 7 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   HRun           0       ins  
 8 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   HaplotypeScore 41.3338 ins  
 9 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   MQ0            0       ins  
10 1_866511_C_CCCCT 1     866511 rs60722469 C     CCCCT 258.62 PASS   MQ             61.94   ins  

El estado actual de la transformación produce una nueva fila para cada anotación de variante y dos columnas que contienen la clave y el valor. Si queremos nuestros datos en formato ancho donde cada anotación es una columna, podemos usar el pivot_wider función.

En el siguiente código, he usado las primeras ocho columnas (id_cols = vid:filter) para especificar las columnas que identifican de forma única cada variante, es decir, la misma variante tendrá los mismos valores en estas columnas. Especificamos nuestros nombres de columna de la key columna y los valores para estas celdas vendrán de la value columna.

my_vcf %>%
  mutate(info = sub(pattern = ";$", replacement = "", x = .data$info)) %>%
  separate_rows(info, sep = ";\s*") %>%
  separate(info, c('key', 'value'), sep = "=") %>%
  distinct() %>% # remove duplicated annotations, if any
  pivot_wider(id_cols = vid:filter, names_from = key, values_from = value) %>%
  head(10)

Ahora cada fila es una sola variante y cada columna es una variable.

# A tibble: 10 × 28
   vid     chrom    pos id    ref   alt   qual  filter AC    AF    AN    DB    DP    FS    HRun  HaplotypeScore
   <chr>   <fct>  <int> <chr> <chr> <chr> <chr> <chr>  <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>         
 1 1_8665… 1     866511 rs60… C     CCCCT 258.… PASS   2     1.00  2     NA    11    0.000 0     41.3338       
 2 1_8793… 1     879317 rs75… C     T     150.… PASS   1     0.50  2     NA    21    1.984 0     0.0000        
 3 1_8794… 1     879482 .     G     C     484.… PASS   1     0.50  2     NA    48    4.452 0     0.5784        
 4 1_8803… 1     880390 rs37… C     A     288.… PASS   1     0.50  2     NA    29    1.485 0     0.0000        
 5 1_8816… 1     881627 rs22… G     A     486.… PASS   1     0.50  2     NA    33    0.000 1     1.8893        
 6 1_8840… 1     884091 rs75… C     G     65.46 PASS   1     0.50  2     NA    12    0.000 1     0.0000        
 7 1_8841… 1     884101 rs49… A     C     85.81 PASS   1     0.50  2     NA    12    0.000 3     0.0000        
 8 1_8924… 1     892460 rs41… G     C     1736… PASS   1     0.50  2     NA    152   0.000 0     0.0000        
 9 1_8977… 1     897730 rs75… C     T     225.… PASS   1     0.50  2     NA    21    6.419 1     1.8410        
10 1_9092… 1     909238 rs38… G     C     247.… PASS   1     0.50  2     NA    19    0.000 1     0.0000        
# … with 12 more variables: MQ0 <chr>, MQ <chr>, QD <chr>, set <chr>, BaseQRankSum <chr>, Dels <chr>,
#   MQRankSum <chr>, ReadPosRankSum <chr>, DS <chr>, GENE <chr>, INHERITANCE <chr>, MIM <chr>

Dividir la columna de grupo en un archivo GTF

El GTF también tiene una columna repleta de pares clave-valor.

my_gtf <- read_tsv(
  file = "https://github.com/davetang/importbio/raw/master/inst/extdata/gencode.v38.annotation.sample.gtf.gz",
  comment = "#",
  show_col_types = FALSE,
  col_names = c('chr', 'src', 'feat', 'start', 'end', 'score', 'strand', 'frame', 'group'))

my_gtf %>%
  select(group) %>%
  head()
# A tibble: 6 × 1
  group                                                                                                        
  <chr>                                                                                                        
1 "gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; lev…
2 "gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pse…
3 "gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pse…
4 "gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pse…
5 "gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pse…
6 "gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pse…

Podemos usar la misma estrategia (pero con algunos pasos de formato adicionales) para dividir la columna.

my_gtf %>%
  mutate(group = sub(pattern = ";$", replacement = "", x = .data$group)) %>%
  mutate(group = gsub(pattern = '"', replacement = "", x = .data$group)) %>%
  separate_rows(group, sep = ";\s*") %>%
  separate(group, c('key', 'value'), sep = "\s") %>%
  distinct() %>% # remove duplicated annotations, if any
  pivot_wider(id_cols = chr:frame, names_from = key, values_from = value) -> my_gtf_split

head(my_gtf_split, 10)
# A tibble: 10 × 24
   chr   src    feat       start   end score strand frame gene_id gene_type gene_name level hgnc_id havana_gene
   <chr> <chr>  <chr>      <dbl> <dbl> <chr> <chr>  <chr> <list>  <list>    <list>    <lis> <list>  <list>     
 1 chr1  HAVANA gene       11869 14409 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 2 chr1  HAVANA transcript 11869 14409 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 3 chr1  HAVANA exon       11869 12227 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 4 chr1  HAVANA exon       12613 12721 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 5 chr1  HAVANA exon       13221 14409 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 6 chr1  HAVANA transcript 12010 13670 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 7 chr1  HAVANA exon       12010 12057 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 8 chr1  HAVANA exon       12179 12227 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
 9 chr1  HAVANA exon       12613 12697 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
10 chr1  HAVANA exon       12975 13052 .     +      .     <chr>   <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
# … with 10 more variables: transcript_id <list>, transcript_type <list>, transcript_name <list>,
#   transcript_support_level <list>, tag <list>, havana_transcript <list>, exon_number <list>, exon_id <list>,
#   ont <list>, protein_id <list>

Sin embargo, las columnas divididas son listas porque hubo algunos casos en los que había varias anotaciones con la misma clave y se necesitaba una lista para almacenar varios valores. por ejemplo el tag key se repitió más de una vez con diferentes valores únicos para algunas anotaciones.

map_lgl(my_gtf_split$tag, function(x) length(x) > 1)
 [1] FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
[18]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
[35]  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE
[52]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[69]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE
[86] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE

Podemos verificar qué columnas tienen múltiples valores.

check_column <- function(x)

my_check <- map_lgl(my_gtf_split, check_column)
names(my_check[my_check])
[1] "transcript_id"     "transcript_name"   "tag"               "havana_transcript" "exon_number"      
[6] "ont"       

Por lo tanto, a pesar de que solo un subconjunto de las columnas contiene valores múltiples, todas las columnas dinámicas se convirtieron en listas. Sin embargo, podemos volver a convertir las columnas en caracteres, lo que tiene sentido para el gene_id columna que solo contenía valores de caracteres únicos únicos en primer lugar.

my_gtf_split %>%
  mutate(gene_id = as.character(gene_id)) %>%
  head()
# A tibble: 6 × 24
  chr   src    feat       start   end score strand frame gene_id  gene_type gene_name level hgnc_id havana_gene
  <chr> <chr>  <chr>      <dbl> <dbl> <chr> <chr>  <chr> <chr>    <list>    <list>    <lis> <list>  <list>     
1 chr1  HAVANA gene       11869 14409 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
2 chr1  HAVANA transcript 11869 14409 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
3 chr1  HAVANA exon       11869 12227 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
4 chr1  HAVANA exon       12613 12721 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
5 chr1  HAVANA exon       13221 14409 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
6 chr1  HAVANA transcript 12010 13670 .     +      .     ENSG000… <chr [1]> <chr [1]> <chr> <chr>   <chr [1]>  
# … with 10 more variables: transcript_id <list>, transcript_type <list>, transcript_name <list>,
#   transcript_support_level <list>, tag <list>, havana_transcript <list>, exon_number <list>, exon_id <list>,
#   ont <list>, protein_id <list>

Pero también podemos hacer esto con el tag columna (incluso si necesitaba una lista para almacenar los valores múltiples) y las entradas con valores múltiples se convierten en código R (carácter).

my_gtf_split %>%
  mutate(tag = as.character(tag)) %>%
  select(tag) %>%
  head()
# A tibble: 6 × 1
  tag                                  
  <chr>                                
1 "NULL"                               
2 "basic"                              
3 "basic"                              
4 "basic"                              
5 "basic"                              
6 "c("basic", "Ensembl_canonical")"

Si no le importa tener código R (carácter) en sus datos, puede realizar esta transformación en todas las columnas dinámicas.

my_gtf_split %>%
  mutate(across(gene_id:protein_id, as.character))
# A tibble: 94 × 24
   chr   src    feat       start   end score strand frame gene_id gene_type gene_name level hgnc_id havana_gene
   <chr> <chr>  <chr>      <dbl> <dbl> <chr> <chr>  <chr> <chr>   <chr>     <chr>     <chr> <chr>   <chr>      
 1 chr1  HAVANA gene       11869 14409 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 2 chr1  HAVANA transcript 11869 14409 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 3 chr1  HAVANA exon       11869 12227 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 4 chr1  HAVANA exon       12613 12721 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 5 chr1  HAVANA exon       13221 14409 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 6 chr1  HAVANA transcript 12010 13670 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 7 chr1  HAVANA exon       12010 12057 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 8 chr1  HAVANA exon       12179 12227 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
 9 chr1  HAVANA exon       12613 12697 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
10 chr1  HAVANA exon       12975 13052 .     +      .     ENSG00… transcri… DDX11L1   2     HGNC:3… OTTHUMG000…
# … with 84 more rows, and 10 more variables: transcript_id <chr>, transcript_type <chr>,
#   transcript_name <chr>, transcript_support_level <chr>, tag <chr>, havana_transcript <chr>,
#   exon_number <chr>, exon_id <chr>, ont <chr>, protein_id <chr>

Resumen

Los siguientes pasos se pueden usar para dividir una columna que contiene pares clave-valor en columnas separadas:

  1. Usar separate_rows para dividir una sola columna en filas
  2. Usar separate para dividir un par clave-valor en dos columnas
  3. Usar pivot_wider para convertir la tabla de formato largo de nuevo a formato ancho

Sin embargo, a veces los datos se empaquetan en una sola columna porque, en primer lugar, no se pueden formatear bien.

Imprimir amigable, PDF y correo electrónico

Fuente del artículo

Deja un comentario