0.1 Acerca del curso

El tutorial a continuación fue creado especialmente para guiar el trabajo práctico del curso pre-congreso ISME Latin America 2019: Análisis de datos bioinformáticos para metagenomas y amplicones usando R. A realizarse el próximo 9 y 10 de septiembre en la Universidad Técnica Federico Santa María, Valparaíso, Chile.

Ir a la página de inicio del curso

El curso cuenta con 6 unidades. Ahora, usted se encuentra en la unidad: Introducción a R: Manipulación de datos y visualización

Otras unidades del curso son:

Análisis de secuencias de 16S con DADA2

Introducción a phyloseq y a análisis de diversidad

Búsqueda de genes de interés en datos de metagenómica shotgun

Visualización y curación de genomas ensamblados desde metagenomas (MAGs)

Redes de co-ocurrencia de microorganismos


CASTRO LAB

Centro de Bioinformática y Biología Integrativa | Universidad Andrés Bello

Santiago, Chile

1 Introducción a R

1.1 ¿Por qué R?

Tus análisis serán una serie de comandos escritos (R script), lo que es muy conveniente porque:

  • Conlleva a tener el control y completo entendimiento de qué estás haciendo (adiós a la caja negra en tus análisis in silico).
  • ¡Reproducibilidad! Si necesitas repetir el análisis, agregar más datos o corregir un error, sólo necesitas correr tu script nuevamente y tus tests estadísticos y gráficos se actualizarán automáticamente.
  • Tus análisis son claros y transparentes, cualquier colega o tú en el futuro pueden leerlo para encontrar errores y/o hacer mejoras.
  • Puedes comentar tu script. Agregar comentarios, explicaciones o relatos en tu código facilitará entender, paso a paso, el por qué de tus análisis a tus colegas, e incluso a ti mísmo en el futuro.

Otras importantes ventajas de usar R son:

  • ¡Escalabilidad! R está diseñado para el análisis de datos. Las habilidades que vayas aprendiendo en R se pueden escalar fácilmente con el tamaño de tus datos (set de datos con cientos a miles, incluso millones de líneas).
  • ¡Gráficos de alta calidad! R tiene innumerables funcionalidades para todo tipo de gráficas a tu disposición para lograr una efectiva visualización de tus resultados. Visita The R Graph Gallery para ver ejemplos de lo que puedes llegar a hacer.
  • ¡Amplia documentación! R cuenta con una basta documentación y la web está minada de tutoriales, sólo debes buscar. Por ejemplo, Cookbook for R que provee de soluciones a tareas y problemas básicos en análisis de datos.

Aprender R puede ser difícil y tomar tiempo en un principio, pero ¡no te preocupes, no estás sólo! R cuenta con una amplia comunidad de usuarios dispuestos a ayudar a través de mailing-lists y websites como Stack Overflow y RStudio community. Sin embargo, lo más probable es que tus dudas ya se encuentren resueltas en la web, asi que primero haz una búsqueda en Google usando palabras clave o copia y pega en el buscador el mensaje de error que te aparezca en la consola de R. Consejo: haz tus búsquedas y consultas en inglés, tendrás acceso a más y mejor información.

1.2 Conociendo RStudio

¡Todo lo que necesitas en una ventana! Como puedes ver en la imagen a continuación, RStudio se divide en 4 paneles principales:

  1. Panel superior izquierdo: editor de texto para escribir comandos/instrucciónes en R (R script), aquí puedes agregar comentarios.
  2. Panel inferior izquierdo: consola de R, donde se ejecutan las líneas de comando en tu script.
  3. Panel superior derecho: aquí puedes visualizar los datos presentes en la memoria de R.
  4. Panel inferior derecho: aquí encuentras 5 pestañas Files/Plots/Packages/Help/Viewer en las que puedes navegar por tus archivos, visualizar tus gráficos en tiempo real, administrar paquetes de R, y pedir ayuda.

RStudio RStudio interfaz

Básicamente, lo que haremos es escribir instrucciones (a lo que llamamos comandos) en el lenguaje de programación R, y luego le indicaremos al computador que ejecute la instrucción (a lo que llamamos ejecutar o correr). En RStudio puedes escribir los comandos en el editor de texto y luego ejecutarlos, una línea de código a la vez, presionando Ctrl + Enter (cmd + return en Mac). Cuando ejecutas comandos, es importante poner atención a la consola. Si hubo algún problema con la ejecución de la instrucción, aparecerá un mensaje de error en la consola. Cuando el computador finaliza la instrucción, un nuevo > aparece en la consola indicando que a finalizado de ejecutar el último comando y está listo para correr el siguiente.

1.3 Set de datos para trabajar

Para ésta sección vamos a trabajar con una tabla de datos que contiene información acerca de índices de diversidad de microorganismos presentes en 87 muestras de piel de 3 especies de ballena recolectadas en 3 zonas de Chile (Megaptera novaeangliae, Balaenoptera musculus, Balaenoptera physalus; Estrecho Magallanes, Chiloé, Reserva Nacional Pinguino de Humboldt).

  • Descarga la tabla de datos IR_table1.csv AQUÍ.

  • Usa la función read.table() para leer o cargar la tabla de datos a la memoria de R. La forma de cargar nuevos datos a la memoria de R, es asignar dicho set de datos a un objeto usando <- (objeto <- datos).

data <- read.table(file = "~/Dropbox/CastroLab_database/workshops_data/IR_table1.csv", sep = ",", header = TRUE)
# usa la ruta correspondiente en tu computadora hasta el archivo 'IR_table1.csv'
# ahora los datos se encuentran guardados en el objeto "data"

file = "", sep = "" y header = TRUE/FALSE son argumentos de la función read.table. Cada vez que llamamos una función debemos usar argumentos para indicar datos de entrada y/o preferencias:

  • file = "~/Dropbox/CastroLab_database/workshops_data/IR_table1.csv" para indicar archivo del cual queremos copiar datos y cargarlos a la memoria de R.
  • sep = "," para indicar el separador de columnas de la tabla.
  • header = TRUE para indicar que la tabla “IR_table1.csv” sí (“TRUE”) contiene headers o títulos de columna, de lo contrario usaríamos “FALSE”.

Puedes conocer la descripción de cualquier función y sus argumentos escribiendo “?” + “nombre de la función” en la consola de R (e.g., ?read.table).

Cuando cargamos una tabla de datos usando la función read.table(), ésta pasa a conformar lo que se conoce como un data frame en la memoria de R. Un data frame es una representación de los datos en forma de tabla, donde las columnas son vectores, todos del mismo largo (igual número de filas). Un vector es el tipo de dato más básico en R, está compuesto por una serie de valores los que pueden ser números o caracteres.

  • Usa la función View() para ver el contenido del objeto data.
View(data)
  • Prueba las siguientes funciones para inspeccionar nuestro nuevo data frame “data”:
dim(data) # muestra número de filas y columnas, respectivamente
## [1] 87 30
nrow(data) # muestra número de filas
## [1] 87
ncol(data) # muestra número de columnas
## [1] 30
head(data) # muestra las primeras 6 filas
##    sample_ID           geo_loc_name                species observed
## 1 SRR6442697 Estrecho de Magallanes Megaptera novaeangliae       31
## 2 SRR6442698 Estrecho de Magallanes Megaptera novaeangliae       33
## 3 SRR6442699 Estrecho de Magallanes Megaptera novaeangliae       43
## 4 SRR6442700 Estrecho de Magallanes Megaptera novaeangliae       29
## 5 SRR6442701 Estrecho de Magallanes Megaptera novaeangliae       26
## 6 SRR6442702 Estrecho de Magallanes Megaptera novaeangliae       21
##    shannon richness_0 richness_20 richness_50 richness_80
## 1 2.121629         31          31          31          31
## 2 1.499137         33          33          33          33
## 3 2.249373         43          43          43          43
## 4 1.306574         29          29          29          29
## 5 1.077438         26          26          26          26
## 6 1.131204         21          21          21          21
##   diversities_inverse_simpson diversities_gini_simpson diversities_shannon
## 1                    6.112104                0.8363902            2.121629
## 2                    2.505873                0.6009375            1.499137
## 3                    7.388570                0.8646558            2.249373
## 4                    2.644111                0.6218011            1.306574
## 5                    1.873331                0.4661914            1.077438
## 6                    2.740500                0.6351030            1.131204
##   diversities_fisher diversities_coverage evenness_camargo evenness_pielou
## 1           3.907415                    3       0.02929889       0.6178324
## 2           4.844428                    1       0.01815948       0.4287528
## 3           5.368153                    3       0.03391424       0.5980465
## 4           3.094775                    1       0.01389960       0.3880188
## 5           2.956294                    1       0.01143026       0.3306954
## 6           2.010698                    2       0.01179407       0.3715539
##   evenness_simpson evenness_evar evenness_bulla dominance_dbp
## 1      0.028296778    0.08554199     0.06038316     0.3022210
## 2      0.011601264    0.11572176     0.05579809     0.6138276
## 3      0.034206345    0.07397473     0.05992808     0.2424974
## 4      0.012241256    0.05943919     0.04011332     0.5250103
## 5      0.008672827    0.08300059     0.04374973     0.7143736
## 6      0.012687499    0.05585996     0.02383196     0.4903917
##   dominance_dmn dominance_absolute dominance_relative dominance_simpson
## 1     0.4948605               3293          0.3022210         0.1636098
## 2     0.7063907               2699          0.6138276         0.3990625
## 3     0.4386486               3919          0.2424974         0.1353442
## 4     0.8384308              19071          0.5250103         0.3781989
## 5     0.8534960              13936          0.7143736         0.5338086
## 6     0.7767577              33864          0.4903917         0.3648970
##   dominance_core_abundance dominance_gini rarity_log_modulo_skewness
## 1                0.5595631      0.9707011                   2.055226
## 2                0.8749147      0.9818405                   2.060118
## 3                0.5033723      0.9660858                   2.059857
## 4                0.9925120      0.9861004                   2.057839
## 5                0.9824175      0.9885697                   2.056901
## 6                0.9981030      0.9882059                   2.049704
##   rarity_low_abundance rarity_noncore_abundance rarity_rare_abundance
## 1          0.009177680              0.076174743           0.076174743
## 2          0.007505117              0.024789629           0.024789629
## 3          0.005692717              0.282965163           0.282965163
## 4          0.002890571              0.000220234           0.000220234
## 5          0.004459709              0.001127742           0.001127742
## 6          0.003301716              0.000086900           0.000086900
tail(data) # muestra las últimas 6 filas
##     sample_ID           geo_loc_name                species observed
## 82 SRR6442787 Estrecho de Magallanes Megaptera novaeangliae       34
## 83 SRR6442788 Estrecho de Magallanes Megaptera novaeangliae       49
## 84 SRR6442789 Estrecho de Magallanes Megaptera novaeangliae       28
## 85 SRR6442790 Estrecho de Magallanes Megaptera novaeangliae       27
## 86 SRR6442792 Estrecho de Magallanes Megaptera novaeangliae       24
## 87 SRR6442794 Estrecho de Magallanes Megaptera novaeangliae       22
##      shannon richness_0 richness_20 richness_50 richness_80
## 82 1.2165375         34          34          34          34
## 83 1.4508829         49          49          49          49
## 84 0.8779605         28          28          28          28
## 85 1.2337859         27          27          27          27
## 86 1.6592270         24          24          24          24
## 87 1.0682126         22          22          22          22
##    diversities_inverse_simpson diversities_gini_simpson
## 82                    2.370590                0.5781641
## 83                    2.723141                0.6327769
## 84                    2.132005                0.5309578
## 85                    2.365953                0.5773373
## 86                    3.274365                0.6945973
## 87                    2.195027                0.5444248
##    diversities_shannon diversities_fisher diversities_coverage
## 82           1.2165375           3.789528                    1
## 83           1.4508829           5.808398                    1
## 84           0.8779605           3.261216                    1
## 85           1.2337859           3.494975                    1
## 86           1.6592270           3.350736                    1
## 87           1.0682126           2.388742                    1
##    evenness_camargo evenness_pielou evenness_simpson evenness_evar
## 82       0.01240594       0.3449839      0.010974954    0.07523218
## 83       0.01636409       0.3728032      0.012607132    0.10265990
## 84       0.01006692       0.2634774      0.009870391    0.10082495
## 85       0.01243090       0.3743468      0.010953485    0.09603370
## 86       0.01921725       0.5220890      0.015159097    0.10629033
## 87       0.01062587       0.3455833      0.010162163    0.06442340
##    evenness_bulla dominance_dbp dominance_dmn dominance_absolute
## 82     0.02951427     0.6173890     0.7687722              18434
## 83     0.05934702     0.5651963     0.7234154              15132
## 84     0.02545193     0.5161512     0.9659794               9012
## 85     0.03539774     0.6129155     0.8023506               4850
## 86     0.05534884     0.5129630     0.6608796               2216
## 87     0.03291487     0.6099778     0.8939806              14562
##    dominance_relative dominance_simpson dominance_core_abundance
## 82          0.6173890         0.4218359                0.9939380
## 83          0.5651963         0.3672231                0.8883203
## 84          0.5161512         0.4690422                0.9822451
## 85          0.6129155         0.4226627                0.9962088
## 86          0.5129630         0.3054027                0.8837963
## 87          0.6099778         0.4555752                0.9965652
##    dominance_gini rarity_log_modulo_skewness rarity_low_abundance
## 82      0.9875941                   2.052933          0.003181727
## 83      0.9836359                   2.057605          0.015164531
## 84      0.9899331                   2.057562          0.006815578
## 85      0.9875691                   2.059455          0.007456085
## 86      0.9807828                   2.058916          0.004166667
## 87      0.9893741                   2.054947          0.002303858
##    rarity_noncore_abundance rarity_rare_abundance
## 82              0.000837297           0.000837297
## 83              0.014044000           0.014044000
## 84              0.014604811           0.014604811
## 85              0.000884620           0.000884620
## 86              0.001157407           0.001157407
## 87              0.000376995           0.000376995
summary(data) # calcula estadísticas básicas para cada columna
##       sample_ID                                 geo_loc_name
##  SRR6442697: 1   Chiloe                               :28   
##  SRR6442698: 1   Estrecho de Magallanes               :41   
##  SRR6442699: 1   Reserva Nacional Pinguino de Humboldt:18   
##  SRR6442700: 1                                              
##  SRR6442701: 1                                              
##  SRR6442702: 1                                              
##  (Other)   :81                                              
##                    species      observed       shannon      
##  Balaenoptera musculus :26   Min.   :13.0   Min.   :0.1675  
##  Balaenoptera physalus : 7   1st Qu.:28.0   1st Qu.:1.1493  
##  Megaptera novaeangliae:54   Median :36.0   Median :1.4509  
##                              Mean   :41.3   Mean   :1.4610  
##                              3rd Qu.:53.0   3rd Qu.:1.7976  
##                              Max.   :84.0   Max.   :3.1401  
##                                                             
##    richness_0    richness_20    richness_50    richness_80  
##  Min.   :13.0   Min.   :13.0   Min.   :13.0   Min.   :13.0  
##  1st Qu.:28.0   1st Qu.:28.0   1st Qu.:28.0   1st Qu.:28.0  
##  Median :36.0   Median :36.0   Median :36.0   Median :36.0  
##  Mean   :41.3   Mean   :41.3   Mean   :41.3   Mean   :41.3  
##  3rd Qu.:53.0   3rd Qu.:53.0   3rd Qu.:53.0   3rd Qu.:53.0  
##  Max.   :84.0   Max.   :84.0   Max.   :84.0   Max.   :84.0  
##                                                             
##  diversities_inverse_simpson diversities_gini_simpson diversities_shannon
##  Min.   : 1.052              Min.   :0.04979          Min.   :0.1675     
##  1st Qu.: 2.164              1st Qu.:0.53769          1st Qu.:1.1493     
##  Median : 2.814              Median :0.64461          Median :1.4509     
##  Mean   : 3.478              Mean   :0.63208          Mean   :1.4610     
##  3rd Qu.: 4.137              3rd Qu.:0.75821          3rd Qu.:1.7976     
##  Max.   :14.576              Max.   :0.93139          Max.   :3.1401     
##                                                                          
##  diversities_fisher diversities_coverage evenness_camargo  
##  Min.   : 1.620     Min.   :1.000        Min.   :0.005441  
##  1st Qu.: 3.286     1st Qu.:1.000        1st Qu.:0.012368  
##  Median : 4.056     Median :2.000        Median :0.016613  
##  Mean   : 5.018     Mean   :1.701        Mean   :0.018888  
##  3rd Qu.: 6.263     3rd Qu.:2.000        3rd Qu.:0.023079  
##  Max.   :13.163     Max.   :5.000        Max.   :0.083105  
##                                                            
##  evenness_pielou   evenness_simpson   evenness_evar     evenness_bulla   
##  Min.   :0.04711   Min.   :0.004872   Min.   :0.05064   Min.   :0.01145  
##  1st Qu.:0.32381   1st Qu.:0.010016   1st Qu.:0.06995   1st Qu.:0.03536  
##  Median :0.38802   Median :0.013027   Median :0.09603   Median :0.04923  
##  Mean   :0.40233   Mean   :0.016100   Mean   :0.10164   Mean   :0.05415  
##  3rd Qu.:0.49853   3rd Qu.:0.019151   3rd Qu.:0.12267   3rd Qu.:0.06665  
##  Max.   :0.74418   Max.   :0.067482   Max.   :0.25815   Max.   :0.19846  
##                                                                          
##  dominance_dbp    dominance_dmn    dominance_absolute dominance_relative
##  Min.   :0.1358   Min.   :0.2602   Min.   :  382      Min.   :0.1358    
##  1st Qu.:0.3772   1st Qu.:0.6092   1st Qu.: 4856      1st Qu.:0.3772    
##  Median :0.4841   Median :0.7482   Median :13936      Median :0.4841    
##  Mean   :0.5042   Mean   :0.7221   Mean   :16988      Mean   :0.5042    
##  3rd Qu.:0.6167   3rd Qu.:0.8478   3rd Qu.:20533      3rd Qu.:0.6167    
##  Max.   :0.9747   Max.   :0.9932   Max.   :80901      Max.   :0.9747    
##                                                                         
##  dominance_simpson dominance_core_abundance dominance_gini  
##  Min.   :0.06861   Min.   :0.0063           Min.   :0.9169  
##  1st Qu.:0.24179   1st Qu.:0.4664           1st Qu.:0.9769  
##  Median :0.35539   Median :0.9189           Median :0.9834  
##  Mean   :0.36792   Mean   :0.7313           Mean   :0.9811  
##  3rd Qu.:0.46231   3rd Qu.:0.9736           3rd Qu.:0.9876  
##  Max.   :0.95021   Max.   :0.9988           Max.   :0.9946  
##                                                             
##  rarity_log_modulo_skewness rarity_low_abundance rarity_noncore_abundance
##  Min.   :1.976              Min.   :0.001072     Min.   :0.0000869       
##  1st Qu.:2.054              1st Qu.:0.004414     1st Qu.:0.0017489       
##  Median :2.058              Median :0.007645     Median :0.0062500       
##  Mean   :2.053              Mean   :0.010067     Mean   :0.0630320       
##  3rd Qu.:2.060              3rd Qu.:0.013359     3rd Qu.:0.0252855       
##  Max.   :2.061              Max.   :0.039322     Max.   :0.9751045       
##                                                                          
##  rarity_rare_abundance
##  Min.   :0.0000869    
##  1st Qu.:0.0017489    
##  Median :0.0062500    
##  Mean   :0.0630320    
##  3rd Qu.:0.0252855    
##  Max.   :0.9751045    
## 

1.4 Extraer información de tablas

Nuestra tabla de datos de estudio (data) consta de filas y columnas (2 dimensiones), si queremos extraer algunos datos de interés, debemos especificar las “coordinadas” de los datos que queremos obtener. Primero el número(s) de fila, seguido por el número(s) de columna. Existen diferentes formas de especificar coordenadas, que nos llevaran a obtener datos de diferente clase o tipo.

  • Primer elemento en la primera columna del data frame, como vector:
data[1, 1]
  • Primer elemento en la sexta columna, como vector:
data[1, 6]
  • Primera columna del data frame, como vector:
data[, 1]
  • Primera columna del data frame, como data frame:
data[1]
  • Primeros tres elementos en la séptima columna, como vector:
data[1:3, 7] # filas 1 a 3 de la columna 7
  • La tercera fila del data frame, como data frame:
data[3, ]
  • Equivalente a la función head():
data[1:6, ] 
  • También puedes excluir datos usando el símbolo - (“todo menos…”):
data[, -1]    # todo excepto la primera columna
data[-c(7:87), ] # equivalente a head(data)

La función c() se utiliza para indicar una serie de valores o asignar una serie de valores a un vector (e.g., peso_kg <- c(60,55,64,80,74)).

  • También puedes extraer datos usando el nombre de las columnas o headers:
data["species"]    # resultado como data.frame
data[, "species"]    # resultado como vector
data[["species"]]    # resultado como vector
data$species    # resultado como vector

Puedes revisar los headers disponibles en tu data frame usando las funciones colnames(data) o View(data).

RStudio cuenta con una muy útil función de autocompletado, presiona tab (tabulador) para obtener nombres completos y correctos de funciones, columnas (headers), etc.

1.5 Factores

  • La función str() muestra la estructura de un objeto e información acerca de la clase y contenido de cada columna:
str(data)
## 'data.frame':    87 obs. of  30 variables:
##  $ sample_ID                  : Factor w/ 87 levels "SRR6442697","SRR6442698",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ geo_loc_name               : Factor w/ 3 levels "Chiloe","Estrecho de Magallanes",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ species                    : Factor w/ 3 levels "Balaenoptera musculus",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ observed                   : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ shannon                    : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ richness_0                 : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_20                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_50                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_80                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ diversities_inverse_simpson: num  6.11 2.51 7.39 2.64 1.87 ...
##  $ diversities_gini_simpson   : num  0.836 0.601 0.865 0.622 0.466 ...
##  $ diversities_shannon        : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ diversities_fisher         : num  3.91 4.84 5.37 3.09 2.96 ...
##  $ diversities_coverage       : int  3 1 3 1 1 2 2 1 2 1 ...
##  $ evenness_camargo           : num  0.0293 0.0182 0.0339 0.0139 0.0114 ...
##  $ evenness_pielou            : num  0.618 0.429 0.598 0.388 0.331 ...
##  $ evenness_simpson           : num  0.0283 0.0116 0.03421 0.01224 0.00867 ...
##  $ evenness_evar              : num  0.0855 0.1157 0.074 0.0594 0.083 ...
##  $ evenness_bulla             : num  0.0604 0.0558 0.0599 0.0401 0.0437 ...
##  $ dominance_dbp              : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_dmn              : num  0.495 0.706 0.439 0.838 0.853 ...
##  $ dominance_absolute         : int  3293 2699 3919 19071 13936 33864 5957 30937 10163 35049 ...
##  $ dominance_relative         : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_simpson          : num  0.164 0.399 0.135 0.378 0.534 ...
##  $ dominance_core_abundance   : num  0.56 0.875 0.503 0.993 0.982 ...
##  $ dominance_gini             : num  0.971 0.982 0.966 0.986 0.989 ...
##  $ rarity_log_modulo_skewness : num  2.06 2.06 2.06 2.06 2.06 ...
##  $ rarity_low_abundance       : num  0.00918 0.00751 0.00569 0.00289 0.00446 ...
##  $ rarity_noncore_abundance   : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
##  $ rarity_rare_abundance      : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...

Como podrás ver en el output de str(data) que las columnas sample_ID, geo_loc_name y species son de una clase llamada Factor. Factores representan datos categóricos. Son guardados en la memoria de R como números enteros (integers), los que pueden estar ordenados o desordenados.

Los factores contienen un set de valores pre-definidos, conocidos como levels. Por defecto, R ordena los levels en orden alfabético. Por ejemplo, en nuestro objeto data la columna species es un Factor con 3 levels:

levels(data$species)
## [1] "Balaenoptera musculus"  "Balaenoptera physalus" 
## [3] "Megaptera novaeangliae"
nlevels(data$species)
## [1] 3

R asigna 1 al level “Balaenoptera musculus”, 2 al level “Balaenoptera physalus” y 3 al level “Megaptera novaeangliae”.

Algunas veces, el orden de los factores no importa, pero otras veces vamos a requerir especificar el orden porque es importante para el análisis o visualización de los datos. Una forma de re-ordenar los levels del factor species es:

levels(data$species)    # orden actual
## [1] "Balaenoptera musculus"  "Balaenoptera physalus" 
## [3] "Megaptera novaeangliae"
data$species <- factor(data$species, levels = c("Megaptera novaeangliae", "Balaenoptera musculus", "Balaenoptera physalus"))
levels(data$species)    # después de re-ordenar
## [1] "Megaptera novaeangliae" "Balaenoptera musculus" 
## [3] "Balaenoptera physalus"

1.5.1 Transformar factores

  • Para transformar un factor a un vector de caracteres, puedes usar la función as.character():
as.character(data$species)
  • Transformar o convertir factores cuyos niveles son números (e.g., años) a un vector numérico es un poco más complejo. La función as.numeric() muestra los números enteros asignados a cada level, no los niveles en si. Una manera de evitarlo es convertir los factores a caracteres, y luego a números:
year_fct <- factor(c(1990, 1983, 1977, 1998, 1990))
as.numeric(year_fct)    # equivocado! sin mensajes de error...
## [1] 3 2 1 4 3
as.numeric(as.character(year_fct))    # funciona!
## [1] 1990 1983 1977 1998 1990
# otra forma es usar la función levels()
as.numeric(levels(year_fct))[year_fct]    # funciona!
## [1] 1990 1983 1977 1998 1990

1.5.2 Renombrar factores

¿Necesitas renombrar algún elemento en tus datos? Supongamos que queremos cambiar el nombre de la especie “Megaptera novaeangliae” por el nombre común “ballena jorobada”.

species <- data$species    # copiamos la columna "species" en un objeto aparte para no alterar nuestro set de datos original
head(species)
## [1] Megaptera novaeangliae Megaptera novaeangliae Megaptera novaeangliae
## [4] Megaptera novaeangliae Megaptera novaeangliae Megaptera novaeangliae
## 3 Levels: Megaptera novaeangliae ... Balaenoptera physalus
levels(species)    # identifica la posición del level que quieres renombrar (1)
## [1] "Megaptera novaeangliae" "Balaenoptera musculus" 
## [3] "Balaenoptera physalus"
levels(species)[1] <- "ballena jorobada"
levels(species)
## [1] "ballena jorobada"      "Balaenoptera musculus" "Balaenoptera physalus"
head(species)
## [1] ballena jorobada ballena jorobada ballena jorobada ballena jorobada
## [5] ballena jorobada ballena jorobada
## 3 Levels: ballena jorobada ... Balaenoptera physalus
# también puedes hacerlo para los otros dos levels
levels(species)[2] <- "ballena azul"
levels(species)[3] <- "ballena de aleta"
levels(species)
## [1] "ballena jorobada" "ballena azul"     "ballena de aleta"

1.5.3 Argumento stringsAsFactors

Por defecto, al importar un data frame en R, las columnas que contienen caracteres (i.e. texto) son convertidas en factores. Dependiendo de qué queramos hacer con los datos, en algún caso podríamos necesitar que la columna se mantenga como caracter. Para ésto, la función read.table() tiene disponible el argumento stringsAsFactors que puede ser definido como “FALSE” (stringsAsFactors = FALSE).

  • Compara la diferencia entre la tabla de datos de estudio leída como factor vs. caracter:
data <- read.table("~/Dropbox/CastroLab_database/workshops_data/IR_table1.csv", sep = ",", header = TRUE, stringsAsFactors = TRUE)
str(data)
## 'data.frame':    87 obs. of  30 variables:
##  $ sample_ID                  : Factor w/ 87 levels "SRR6442697","SRR6442698",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ geo_loc_name               : Factor w/ 3 levels "Chiloe","Estrecho de Magallanes",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ species                    : Factor w/ 3 levels "Balaenoptera musculus",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ observed                   : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ shannon                    : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ richness_0                 : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_20                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_50                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_80                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ diversities_inverse_simpson: num  6.11 2.51 7.39 2.64 1.87 ...
##  $ diversities_gini_simpson   : num  0.836 0.601 0.865 0.622 0.466 ...
##  $ diversities_shannon        : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ diversities_fisher         : num  3.91 4.84 5.37 3.09 2.96 ...
##  $ diversities_coverage       : int  3 1 3 1 1 2 2 1 2 1 ...
##  $ evenness_camargo           : num  0.0293 0.0182 0.0339 0.0139 0.0114 ...
##  $ evenness_pielou            : num  0.618 0.429 0.598 0.388 0.331 ...
##  $ evenness_simpson           : num  0.0283 0.0116 0.03421 0.01224 0.00867 ...
##  $ evenness_evar              : num  0.0855 0.1157 0.074 0.0594 0.083 ...
##  $ evenness_bulla             : num  0.0604 0.0558 0.0599 0.0401 0.0437 ...
##  $ dominance_dbp              : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_dmn              : num  0.495 0.706 0.439 0.838 0.853 ...
##  $ dominance_absolute         : int  3293 2699 3919 19071 13936 33864 5957 30937 10163 35049 ...
##  $ dominance_relative         : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_simpson          : num  0.164 0.399 0.135 0.378 0.534 ...
##  $ dominance_core_abundance   : num  0.56 0.875 0.503 0.993 0.982 ...
##  $ dominance_gini             : num  0.971 0.982 0.966 0.986 0.989 ...
##  $ rarity_log_modulo_skewness : num  2.06 2.06 2.06 2.06 2.06 ...
##  $ rarity_low_abundance       : num  0.00918 0.00751 0.00569 0.00289 0.00446 ...
##  $ rarity_noncore_abundance   : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
##  $ rarity_rare_abundance      : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
data <- read.table("~/Dropbox/CastroLab_database/workshops_data/IR_table1.csv", sep = ",", header = TRUE, stringsAsFactors = FALSE)
str(data)
## 'data.frame':    87 obs. of  30 variables:
##  $ sample_ID                  : chr  "SRR6442697" "SRR6442698" "SRR6442699" "SRR6442700" ...
##  $ geo_loc_name               : chr  "Estrecho de Magallanes" "Estrecho de Magallanes" "Estrecho de Magallanes" "Estrecho de Magallanes" ...
##  $ species                    : chr  "Megaptera novaeangliae" "Megaptera novaeangliae" "Megaptera novaeangliae" "Megaptera novaeangliae" ...
##  $ observed                   : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ shannon                    : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ richness_0                 : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_20                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_50                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_80                : int  31 33 43 29 26 21 37 32 28 30 ...
##  $ diversities_inverse_simpson: num  6.11 2.51 7.39 2.64 1.87 ...
##  $ diversities_gini_simpson   : num  0.836 0.601 0.865 0.622 0.466 ...
##  $ diversities_shannon        : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ diversities_fisher         : num  3.91 4.84 5.37 3.09 2.96 ...
##  $ diversities_coverage       : int  3 1 3 1 1 2 2 1 2 1 ...
##  $ evenness_camargo           : num  0.0293 0.0182 0.0339 0.0139 0.0114 ...
##  $ evenness_pielou            : num  0.618 0.429 0.598 0.388 0.331 ...
##  $ evenness_simpson           : num  0.0283 0.0116 0.03421 0.01224 0.00867 ...
##  $ evenness_evar              : num  0.0855 0.1157 0.074 0.0594 0.083 ...
##  $ evenness_bulla             : num  0.0604 0.0558 0.0599 0.0401 0.0437 ...
##  $ dominance_dbp              : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_dmn              : num  0.495 0.706 0.439 0.838 0.853 ...
##  $ dominance_absolute         : int  3293 2699 3919 19071 13936 33864 5957 30937 10163 35049 ...
##  $ dominance_relative         : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_simpson          : num  0.164 0.399 0.135 0.378 0.534 ...
##  $ dominance_core_abundance   : num  0.56 0.875 0.503 0.993 0.982 ...
##  $ dominance_gini             : num  0.971 0.982 0.966 0.986 0.989 ...
##  $ rarity_log_modulo_skewness : num  2.06 2.06 2.06 2.06 2.06 ...
##  $ rarity_low_abundance       : num  0.00918 0.00751 0.00569 0.00289 0.00446 ...
##  $ rarity_noncore_abundance   : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
##  $ rarity_rare_abundance      : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...

2 Manipulación de datos en R

Las funciones que hemos estado usando hasta ahora, como read.table y str(), vienen incorporadas en R. Pero existen incontables funciones en R a las que puedes acceder instalando lo que se conoce como paquetes (packages). Antes de usar un paquete por primera vez, primero debes instalarlo y luego debes importarlo en cada sesión de R en la que quieras usarlo.

Hay dos formas de instalar paquetes en R, a través de la consola y a través de la opción Install en la pestaña Packages de RStudio:

  1. Para instalar e importar el paquete tidyverse usando la consola, escribe:
install.packages("tidyverse")
library(tidyverse)
  1. Para instalar e importar el paquete tidyverse usando RStudio, dirígete a la pestaña Packages -> Install -> escribe el nombre del paquete y haz clic en “Install”. Una vez instalado, para importarlo sólo debes identificarlo y seleccionarlo en la lista de paquetes.

install_packages Instalación de paquetes en RStudio

tidyverse es un paquete que incluye la instalación varios paquetes ´útiles para la manipulaci´ón y an´álisis de datos, tales como dplyr, tidyr, ggplot2, etc. Anteriormente practicamos como extraer información de tablas usando [ ], a continuación vamos a utilizar los paquetes dplyr y tidyr para manipular nuestra tabla de datos más fácil, rápido y con más funcionalidades.

  • Vamos a leer nuevamente nuestra tabla de datos, ésta vez usando la función read_csv() del paquete tidyverse:
data <- read_csv("~/Dropbox/CastroLab_database/workshops_data/IR_table1.csv")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   sample_ID = col_character(),
##   geo_loc_name = col_character(),
##   species = col_character()
## )
## See spec(...) for full column specifications.
# Usa la función str() para inspeccionar los datos
str(data)
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 87 obs. of  30 variables:
##  $ sample_ID                  : chr  "SRR6442697" "SRR6442698" "SRR6442699" "SRR6442700" ...
##  $ geo_loc_name               : chr  "Estrecho de Magallanes" "Estrecho de Magallanes" "Estrecho de Magallanes" "Estrecho de Magallanes" ...
##  $ species                    : chr  "Megaptera novaeangliae" "Megaptera novaeangliae" "Megaptera novaeangliae" "Megaptera novaeangliae" ...
##  $ observed                   : num  31 33 43 29 26 21 37 32 28 30 ...
##  $ shannon                    : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ richness_0                 : num  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_20                : num  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_50                : num  31 33 43 29 26 21 37 32 28 30 ...
##  $ richness_80                : num  31 33 43 29 26 21 37 32 28 30 ...
##  $ diversities_inverse_simpson: num  6.11 2.51 7.39 2.64 1.87 ...
##  $ diversities_gini_simpson   : num  0.836 0.601 0.865 0.622 0.466 ...
##  $ diversities_shannon        : num  2.12 1.5 2.25 1.31 1.08 ...
##  $ diversities_fisher         : num  3.91 4.84 5.37 3.09 2.96 ...
##  $ diversities_coverage       : num  3 1 3 1 1 2 2 1 2 1 ...
##  $ evenness_camargo           : num  0.0293 0.0182 0.0339 0.0139 0.0114 ...
##  $ evenness_pielou            : num  0.618 0.429 0.598 0.388 0.331 ...
##  $ evenness_simpson           : num  0.0283 0.0116 0.03421 0.01224 0.00867 ...
##  $ evenness_evar              : num  0.0855 0.1157 0.074 0.0594 0.083 ...
##  $ evenness_bulla             : num  0.0604 0.0558 0.0599 0.0401 0.0437 ...
##  $ dominance_dbp              : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_dmn              : num  0.495 0.706 0.439 0.838 0.853 ...
##  $ dominance_absolute         : num  3293 2699 3919 19071 13936 ...
##  $ dominance_relative         : num  0.302 0.614 0.242 0.525 0.714 ...
##  $ dominance_simpson          : num  0.164 0.399 0.135 0.378 0.534 ...
##  $ dominance_core_abundance   : num  0.56 0.875 0.503 0.993 0.982 ...
##  $ dominance_gini             : num  0.971 0.982 0.966 0.986 0.989 ...
##  $ rarity_log_modulo_skewness : num  2.06 2.06 2.06 2.06 2.06 ...
##  $ rarity_low_abundance       : num  0.00918 0.00751 0.00569 0.00289 0.00446 ...
##  $ rarity_noncore_abundance   : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
##  $ rarity_rare_abundance      : num  0.07617 0.02479 0.28297 0.00022 0.00113 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   sample_ID = col_character(),
##   ..   geo_loc_name = col_character(),
##   ..   species = col_character(),
##   ..   observed = col_double(),
##   ..   shannon = col_double(),
##   ..   richness_0 = col_double(),
##   ..   richness_20 = col_double(),
##   ..   richness_50 = col_double(),
##   ..   richness_80 = col_double(),
##   ..   diversities_inverse_simpson = col_double(),
##   ..   diversities_gini_simpson = col_double(),
##   ..   diversities_shannon = col_double(),
##   ..   diversities_fisher = col_double(),
##   ..   diversities_coverage = col_double(),
##   ..   evenness_camargo = col_double(),
##   ..   evenness_pielou = col_double(),
##   ..   evenness_simpson = col_double(),
##   ..   evenness_evar = col_double(),
##   ..   evenness_bulla = col_double(),
##   ..   dominance_dbp = col_double(),
##   ..   dominance_dmn = col_double(),
##   ..   dominance_absolute = col_double(),
##   ..   dominance_relative = col_double(),
##   ..   dominance_simpson = col_double(),
##   ..   dominance_core_abundance = col_double(),
##   ..   dominance_gini = col_double(),
##   ..   rarity_log_modulo_skewness = col_double(),
##   ..   rarity_low_abundance = col_double(),
##   ..   rarity_noncore_abundance = col_double(),
##   ..   rarity_rare_abundance = col_double()
##   .. )

Como podrás notar, la función read_csv() hace algunos cambios al cargar la tabla con respecto a lo que revisamos anteriormente usando read.table(). Las diferencias son:

  1. Al leer la tabla, muestra un resumen del tipo de dato de cada columna, y sólo muestra las primeras filas y tantas columnas como se puedan visualizar en la pantalla.
  2. Las columnas de clase character (caracteres) no son convertidas en factores.

A continuación vamos a aprender algunas de las funciones más comunes de dplyr:

  • select(): extraer columnas.
  • filter(): extraer filas según condiciones.
  • mutate(): crear nuevas columnas usando la información de otras columnas.
  • group_by() y summarize(): cálcula estadísticas en datos agrupados.
  • arrange(): ordena resultados.
  • count(): conteo de datos.

2.1 Seleccionar columnas y filtrar filas

  • Para seleccionar columnas de una tabla de datos o data frame, vamos a usar la función select(). Luego, para seleccionar filas de acuerdo a cierto criterio (filtrar), vamos a usar la función filter().
select(data, sample_ID, geo_loc_name, species, observed) # los argumentos son, primero el objeto que contiene el data frame, y luego los títulos de las columnas a extraer
## # A tibble: 87 x 4
##    sample_ID  geo_loc_name           species                observed
##    <chr>      <chr>                  <chr>                     <dbl>
##  1 SRR6442697 Estrecho de Magallanes Megaptera novaeangliae       31
##  2 SRR6442698 Estrecho de Magallanes Megaptera novaeangliae       33
##  3 SRR6442699 Estrecho de Magallanes Megaptera novaeangliae       43
##  4 SRR6442700 Estrecho de Magallanes Megaptera novaeangliae       29
##  5 SRR6442701 Estrecho de Magallanes Megaptera novaeangliae       26
##  6 SRR6442702 Estrecho de Magallanes Megaptera novaeangliae       21
##  7 SRR6442703 Estrecho de Magallanes Megaptera novaeangliae       37
##  8 SRR6442704 Estrecho de Magallanes Megaptera novaeangliae       32
##  9 SRR6442705 Estrecho de Magallanes Megaptera novaeangliae       28
## 10 SRR6442706 Estrecho de Magallanes Megaptera novaeangliae       30
## # … with 77 more rows
dplyr::filter(data, geo_loc_name == "Chiloe") # los argumentos son, primero el objeto que contiene el data frame, y luego la columna con el criterio de filtro
## # A tibble: 28 x 30
##    sample_ID geo_loc_name species observed shannon richness_0 richness_20
##    <chr>     <chr>        <chr>      <dbl>   <dbl>      <dbl>       <dbl>
##  1 SRR64427… Chiloe       Megapt…       63   1.28          63          63
##  2 SRR64427… Chiloe       Megapt…       14   0.732         14          14
##  3 SRR64427… Chiloe       Balaen…       18   1.41          18          18
##  4 SRR64427… Chiloe       Balaen…       19   1.68          19          19
##  5 SRR64427… Chiloe       Balaen…       36   2.31          36          36
##  6 SRR64427… Chiloe       Balaen…       63   1.31          63          63
##  7 SRR64427… Chiloe       Balaen…       29   1.73          29          29
##  8 SRR64427… Chiloe       Balaen…       38   1.65          38          38
##  9 SRR64427… Chiloe       Balaen…       81   2.05          81          81
## 10 SRR64427… Chiloe       Balaen…       78   0.725         78          78
## # … with 18 more rows, and 23 more variables: richness_50 <dbl>,
## #   richness_80 <dbl>, diversities_inverse_simpson <dbl>,
## #   diversities_gini_simpson <dbl>, diversities_shannon <dbl>,
## #   diversities_fisher <dbl>, diversities_coverage <dbl>,
## #   evenness_camargo <dbl>, evenness_pielou <dbl>, evenness_simpson <dbl>,
## #   evenness_evar <dbl>, evenness_bulla <dbl>, dominance_dbp <dbl>,
## #   dominance_dmn <dbl>, dominance_absolute <dbl>,
## #   dominance_relative <dbl>, dominance_simpson <dbl>,
## #   dominance_core_abundance <dbl>, dominance_gini <dbl>,
## #   rarity_log_modulo_skewness <dbl>, rarity_low_abundance <dbl>,
## #   rarity_noncore_abundance <dbl>, rarity_rare_abundance <dbl>

Es importante tener en cuenta que funciones provenientes de diferentes paquetes de R pueden tener el mismo nombre. Esto es altamente probable considerando que solo el repositorio CRAN contiene más de 10.000 paquetes. Por lo tanto, en ocasiones es necesario agregar el paquete fuente de la función a utilizar, así: dplyr::filter(). Con el fin de evitar errores asociados a la ejecución de una función no diferente a la deseada, pero que tiene el mismo nombre.

Para saber si el nombre función a utilizar se repite en otros paquetes, escribe en la consola el signo de interrogación seguido por el nombre de la función, así: ?filter Si el nombre se repite en diferentes paquetes, verás una lista en la pestaña Help de RStudio, sino verás la página de ayuda de la función directamente.

¿Quieres seleccionar y filtrar al mismo tiempo? Claro! Hay formas de hacer varias operaciones consecutivas en una misma instrucción, de esta manera evitamos tener que guardar objetos “intermedios” innecesariamente.

2.2 Funciones anidadas y pipes

  • Vamos a anidar funciones (i.e. una función dentro de otra):
nrow(data)
## [1] 87
ncol(data)
## [1] 30
data_div <- select(filter(data, shannon > 1.5), sample_ID, geo_loc_name, species, observed, shannon)
head(data_div)
## # A tibble: 6 x 5
##   sample_ID  geo_loc_name           species                observed shannon
##   <chr>      <chr>                  <chr>                     <dbl>   <dbl>
## 1 SRR6442697 Estrecho de Magallanes Megaptera novaeangliae       31    2.12
## 2 SRR6442699 Estrecho de Magallanes Megaptera novaeangliae       43    2.25
## 3 SRR6442703 Estrecho de Magallanes Megaptera novaeangliae       37    1.86
## 4 SRR6442708 Estrecho de Magallanes Megaptera novaeangliae       27    1.91
## 5 SRR6442709 Estrecho de Magallanes Megaptera novaeangliae       36    1.88
## 6 SRR6442712 Estrecho de Magallanes Megaptera novaeangliae       36    1.94
nrow(data_div) # 39 muestras tienen un índice de shannon mayor a 1.5
## [1] 39
ncol(data_div)
## [1] 5

Es importante recordar que R lee la línea de comando desde dentro hacia fuera. En éste caso, primero se hizo el filtro y luego la selección.

  • Usar funciones anidadas puede ser engorroso cuando quieres hacer muchas operaciones consecutivas, en cuyo caso es convenientes usar pipes. Pipes te permiten usar varias funciones consecutivas, de forma que el output de una función será en input de la siguiente. Los pipes en R lucen así: %>%.
data %>%
  dplyr::filter(shannon > 1.5) %>% 
  dplyr::select(sample_ID, geo_loc_name, species, observed, shannon)
## # A tibble: 39 x 5
##    sample_ID  geo_loc_name           species               observed shannon
##    <chr>      <chr>                  <chr>                    <dbl>   <dbl>
##  1 SRR6442697 Estrecho de Magallanes Megaptera novaeangli…       31    2.12
##  2 SRR6442699 Estrecho de Magallanes Megaptera novaeangli…       43    2.25
##  3 SRR6442703 Estrecho de Magallanes Megaptera novaeangli…       37    1.86
##  4 SRR6442708 Estrecho de Magallanes Megaptera novaeangli…       27    1.91
##  5 SRR6442709 Estrecho de Magallanes Megaptera novaeangli…       36    1.88
##  6 SRR6442712 Estrecho de Magallanes Megaptera novaeangli…       36    1.94
##  7 SRR6442714 Estrecho de Magallanes Megaptera novaeangli…       26    1.81
##  8 SRR6442715 Estrecho de Magallanes Megaptera novaeangli…       35    1.61
##  9 SRR6442722 Chiloe                 Balaenoptera musculus       19    1.68
## 10 SRR6442723 Chiloe                 Balaenoptera musculus       36    2.31
## # … with 29 more rows

Como %>% pasa el objeto de su izquierda como el primer argumento de la función a su derecha, no necesitamos especificar el data frame como el primer argumento de las funciones filter() y select().

2.3 Mutate: crear nuevas columnas desde información existente en otras columnas

En algunas ocasiones necesitamos generar nueva información a partir de la existente, por ejemplo, para hacer conversiones o cálculos matemáticos.

data <- data %>% 
  dplyr::mutate(log10_dom_abs = log10(dominance_absolute))
# Dale un vistazo a la tabla (View(data)), la última columna es la nueva "log10_dom_abs", también notarás que, en la pestaña "Environment" de RStudio, ahora el objeto "data" tiene "31 variables".
# Si deseas tener una vista previa de como quedará tu "data" antes de guardar cualquier cambio, puedes usar pipe para agregar la función head() al final
data %>% 
  dplyr::mutate(log10_dom_abs = log10(dominance_absolute)) %>% 
  head()
## # A tibble: 6 x 31
##   sample_ID geo_loc_name species observed shannon richness_0 richness_20
##   <chr>     <chr>        <chr>      <dbl>   <dbl>      <dbl>       <dbl>
## 1 SRR64426… Estrecho de… Megapt…       31    2.12         31          31
## 2 SRR64426… Estrecho de… Megapt…       33    1.50         33          33
## 3 SRR64426… Estrecho de… Megapt…       43    2.25         43          43
## 4 SRR64427… Estrecho de… Megapt…       29    1.31         29          29
## 5 SRR64427… Estrecho de… Megapt…       26    1.08         26          26
## 6 SRR64427… Estrecho de… Megapt…       21    1.13         21          21
## # … with 24 more variables: richness_50 <dbl>, richness_80 <dbl>,
## #   diversities_inverse_simpson <dbl>, diversities_gini_simpson <dbl>,
## #   diversities_shannon <dbl>, diversities_fisher <dbl>,
## #   diversities_coverage <dbl>, evenness_camargo <dbl>,
## #   evenness_pielou <dbl>, evenness_simpson <dbl>, evenness_evar <dbl>,
## #   evenness_bulla <dbl>, dominance_dbp <dbl>, dominance_dmn <dbl>,
## #   dominance_absolute <dbl>, dominance_relative <dbl>,
## #   dominance_simpson <dbl>, dominance_core_abundance <dbl>,
## #   dominance_gini <dbl>, rarity_log_modulo_skewness <dbl>,
## #   rarity_low_abundance <dbl>, rarity_noncore_abundance <dbl>,
## #   rarity_rare_abundance <dbl>, log10_dom_abs <dbl>

2.4 Dividir -> aplicar -> combinar

Varias operaciones de análisis de datos se pueden realizar, primero dividiendo los datos en grupos, segundo aplicando análisis a cada grupo, y tercero combinando los resultados. Hacemos esto usando las funciones group_by() y summarize() juntas. group_by() toma como argumento los nombres de la columna que contiene valores categóricos, a partir de las cuales queremos hacer algún cálculo. summarize() colapsa cada grupo en una única fila.

  • Vamos a calcular el promedio y desviación estándar (mean(); sd()) del número de taxas observadas (columna “observed”) por zona (columna “geo_loc_name”) y especie (columna “species”) de ballena:
data %>% 
  group_by(geo_loc_name, species) %>% # puedes agrupar por una o múltiples columnas
  dplyr::summarize(mean_observed = mean(observed), # una vez que los datos están agrupados, también puedes aplicar múltiples análisis al mismo tiempo y en múltiples variables
            sd_observed = sd(observed), 
            mean_shannon = mean(shannon), 
            sd_shannon = sd(shannon))
## # A tibble: 5 x 6
## # Groups:   geo_loc_name [3]
##   geo_loc_name   species  mean_observed sd_observed mean_shannon sd_shannon
##   <chr>          <chr>            <dbl>       <dbl>        <dbl>      <dbl>
## 1 Chiloe         Balaeno…          55.0       19.9          1.57      0.568
## 2 Chiloe         Megapte…          38.5       34.6          1.00      0.385
## 3 Estrecho de M… Megapte…          33.5        7.60         1.42      0.457
## 4 Reserva Nacio… Balaeno…          33.9       22.3          1.34      0.717
## 5 Reserva Nacio… Megapte…          43.3       23.6          1.52      0.564
  • Si lo necesitas, también puedes agregar un filtro antes de agrupar los datos y hacer estadística de los grupos. Por ejemplo, supongamos que decidimos no considerar aquellas muestras con menos de 20 taxas (observed < 20):
data %>% 
  dplyr::filter(observed > 20) %>% 
  group_by(geo_loc_name, species) %>% 
  dplyr::summarize(mean_observed_min20 = mean(observed)) %>% 
  print() # para ver el resultado (output) en la consola
## # A tibble: 5 x 3
## # Groups:   geo_loc_name [3]
##   geo_loc_name                      species              mean_observed_min…
##   <chr>                             <chr>                             <dbl>
## 1 Chiloe                            Balaenoptera muscul…               59.7
## 2 Chiloe                            Megaptera novaeangl…               63  
## 3 Estrecho de Magallanes            Megaptera novaeangl…               33.5
## 4 Reserva Nacional Pinguino de Hum… Balaenoptera physal…               46.8
## 5 Reserva Nacional Pinguino de Hum… Megaptera novaeangl…               49.2
  • Muchas veces es útil re-organizar los datos para una más eficiente interpretación de los resultados. Por ejemplo, si queremos ordenar los resultados por número promedio de taxa observado por grupo (“mean_observed”) en orden decreciente:
data %>% 
  group_by(geo_loc_name, species) %>% 
  dplyr::summarize(mean_observed = mean(observed), 
            sd_observed = sd(observed), 
            mean_shannon = mean(shannon), 
            sd_shannon = sd(shannon)) %>% 
  dplyr::arrange(desc(mean_observed)) # la función arrange(), por defecto, ordena los datos en orden creciente, usamos la función desc() para ordenar en orden decreciente
## # A tibble: 5 x 6
## # Groups:   geo_loc_name [3]
##   geo_loc_name   species  mean_observed sd_observed mean_shannon sd_shannon
##   <chr>          <chr>            <dbl>       <dbl>        <dbl>      <dbl>
## 1 Chiloe         Balaeno…          55.0       19.9          1.57      0.568
## 2 Reserva Nacio… Megapte…          43.3       23.6          1.52      0.564
## 3 Chiloe         Megapte…          38.5       34.6          1.00      0.385
## 4 Reserva Nacio… Balaeno…          33.9       22.3          1.34      0.717
## 5 Estrecho de M… Megapte…          33.5        7.60         1.42      0.457

2.5 Contar

La función count() nos permite conocer el número de observaciones por cada variable, o combinación de ellas, en tus datos.

data %>% 
  dplyr::count(species) # ¿cuantas muestras por especie de ballena tenemos?
data %>% 
  dplyr::count(geo_loc_name, sort = TRUE) # argumento "sort = TRUE" para ordenar (decreciente)
data %>% 
  dplyr::count(geo_loc_name, species) # ¿cuantas muestras de cada especie tenemos por zona?

2.6 Remodelar tablas usando spread() y gather()

Con el objetivo de explorar las relaciones entre ciertas variables de interés en nuestros datos, podemos remodelar la tabla de datos de acuerdo a éstas variables.

Supongamos que nos interesa explorar la relación entre la diversidad del microbioma de la piel (“shannon”) de las especies de ballena (“species”) y su locación geográfica.

  • Primero, necesitamos usar group_by() y summarize(), para agrupar nuestras variables de interés y crear una nueva columna con los valores de índice de Shannon promedio para cada grupo. Después, usamos la función spread() para transformar data de modo que: cada especie ahora sea una columna, cada locación geográfica ahora sea una fila, y los valores de Shannon estén en cada celda según corresponda.
data_spread <- data %>% 
  group_by(geo_loc_name, species) %>% 
  dplyr::summarize(mean_shannon = mean(shannon)) %>% 
  spread(key = species, value = mean_shannon, fill = NA)
head(data_spread)
## # A tibble: 3 x 4
## # Groups:   geo_loc_name [3]
##   geo_loc_name       `Balaenoptera mus… `Balaenoptera ph… `Megaptera novae…
##   <chr>                           <dbl>             <dbl>             <dbl>
## 1 Chiloe                           1.57             NA                 1.00
## 2 Estrecho de Magal…              NA                NA                 1.42
## 3 Reserva Nacional …              NA                 1.34              1.52
# Como el número de taxa observada varía por cada muestra, tenemos como resultado varias celdas "NA" ("missing data"). Para éstos casos, la función spread() viene con el argumento "fill".

Ahora vamos a suponer la situación contraria. Inicialmente, tenemos una tabla de datos como data_spread, en la que los nombres de las especies (“species”) son columnas, pero en vez de ello, queremos que las especies sean valores de la variable “species” (columna: “species”).

  • Para lograrlo, necesitamos reunir los nombres de las columnas (especies) y convertirlos en un set de variables:
data_gather <- data_spread %>% 
  gather(key = species, value = mean_shannon, -geo_loc_name) %>% 
  dplyr::filter(!is.na(mean_shannon))
head(data_gather)
## # A tibble: 5 x 3
## # Groups:   geo_loc_name [3]
##   geo_loc_name                          species                mean_shannon
##   <chr>                                 <chr>                         <dbl>
## 1 Chiloe                                Balaenoptera musculus          1.57
## 2 Reserva Nacional Pinguino de Humboldt Balaenoptera physalus          1.34
## 3 Chiloe                                Megaptera novaeangliae         1.00
## 4 Estrecho de Magallanes                Megaptera novaeangliae         1.42
## 5 Reserva Nacional Pinguino de Humboldt Megaptera novaeangliae         1.52

La función is.na() determina si un dato es NA (Not Available). El símbolo ! niega el el resultado. Por lo tanto, al usar !is.na(mean_shannon) estamos pidiendo por valor que no es NA en la columna mean_shannon.

2.7 Exportar datos

El trabajo hecho en una sesión de R (e.g., análisis, nuevas tablas, etc.) sólo existe en la memoria de R mientras la sesión esté abierta. Necesitas exportar los nuevos datos creados para archivarlos en tu computadora. Entonces, así como existe la función read_csv() para leer archivos CSV (comma-separated values) en R, hay una función para generar archivos CSV a partir de tablas de datos en contenidas como objetos en la memoria de R.

  • Por ejemplo, para exportar el data frame que creamos recién usando la función gather() y filter():
# Los argumentos básicos de write_csv() son primero indicar el objeto que se quiere exportar y luego la ruta y nombre del archivo de salida, incluyendo la extensión .csv
write_csv(data_gather, path = "data_gather.csv")

3 Visualización de datos en R

Para comenzar, primero debemos importar el o los paquetes requeridos. ggplot2 está incluido en el paquete tidyverse.

library(tidyverse)

En caso de que la tabla de datos de estudio no esté en la memoria de R, puedes cargarla nuevamente.

data <- read_csv("IR_table1.csv")

ggplot2 es un paquete para graficar, que facilita crear gráficas complejas a partir de datos en un data frame. Incluye varias funciones para especificar que variables graficar, como éstas son expuestas, y varias otras características visuales. ggplot2 funciona mejor con datos extensos, i.e., una columna por cada dimensión, y una fila por cada observación. Tener datos bien estructurados te ayudará a crear figuras más eficientemente.

3.1 Graficar con ggplot2

Las gráficas de ggplot son construidas paso a paso, agregando nuevos elementos cada vez. Para graficar con ggplot, vamos a usar una línea de comando templado que es útil para diferentes tipos de gráficos:

ggplot(data = , mapping = aes()) + ()
  • Usamos la función ggplot() y el argumento data para indicar a partir de qué datos se debe crear la gráfica. Luego, la función aes() (aesthetic) para seleccionar las variables a graficar y como presentarlas, e.g. ejes x e y o características como tamaño, forma, color, etc.

  • ggplot2 ofrece “geoms” para indicar la representación gráfica de los datos (puntos, líneas, barras), algunos de ellos son:

  • geom_point() para gráficos de dispersión.
  • geom_boxplot() para gráficos de caja (boxplots).
  • geom_line() para líneas de tendencia, series de tiempo, etc.

Usa + para agregar un geom a línea de comando de ggplot.

ggplot(data = data, aes(x = shannon, y = evenness_camargo)) + geom_point()

Notas

  • Toda la información o inputs que agregues en la función ggplot() son considerados por cualquier geom (i.e., éstas son configuraciones universales), lo que incluye las variables para el eje x e y indicados usando aes().
  • También puedes especificar inputs para un geom dado, independientemente de los definidos en la función ggplot().
  • El signo + se usa para añadir nuevas funciones y se usa al final de cada línea que contiene la función anterior.

3.2 Construyendo gráficos iterativamente

Construir gráficos con ggplot2 es comúnmente un proceso iterativo. Primero definimos un set de datos a utilizar, luego indicamos las variables a graficar (ejes x e y), y escogemos un geom.

ggplot(data = data, aes(x = observed, y = shannon)) + geom_point()

Seguidamente, empezamos a modificar el gráfico para extraer más información de él y hacerlo más explíito o auto-explicativo. Por ejemplo, podemos agregar transparencia (alpha) para evidenciar la sobre-posición de los puntos (overplotting). También podemos colorear los puntos (color).

ggplot(data = data, aes(x = observed, y = shannon)) + geom_point(alpha = 0.4, color = "blue")

Además, es posible “anotar” tus datos de acuerdo a cierta(s) variable(s) categórica(s), es decir, colorearlos o asignarles formas según ciertas características que tu consideres influyentes en tus datos. Por ejemplo, vamos a anotar nuestra gráfica por especie (“species”) y por locación geográfica (“geo_loc_name”).

ggplot(data = data, aes(x = observed, y = shannon)) + geom_point(alpha = 0.7, aes(color = geo_loc_name, shape = species))

3.3 Boxplot

Boxplots o gráficos de caja, son útiles para visualizar la distribución de los datos de acuerdo a una variable o condición de interés.

  • Por ejemplo, para visualizar la distribución de número de taxas observadas en el microbioma de la piel de las ballenas muestreadas en Chiloé, Estrecho Magallanes y Parque Nacional Pingüino de Humboldt:
ggplot(data = data, aes(x = geo_loc_name, y = observed)) + 
  geom_boxplot()

  • Sin embargo, ¿podría éste resultado estar sesgado por el tamaño muestreal de cada locación geográfica?. Veamos cuántas muestras tenemos por locación geográfica:
data %>% 
  dplyr::count(geo_loc_name, sort = TRUE)
  • Agregar puntos al boxplot nos da una mejor idea del número de muestras y su distribución:
ggplot(data = data, aes(x = geo_loc_name, y = observed)) + 
  geom_boxplot(alpha = 0.5) +
  geom_jitter(alpha = 0.5, color = "tomato")

¿Notas como las cajas (boxplot layer) están detrás de los puntos (jitter layer)? ¿Qué necesitas cambiar en el código de ggplot para que las cajas aparezcan en frente de los puntos? (Pista: el orden sí importa)

3.4 Graficar datos con series de tiempo

Cuando tienes datos tomados en una serie de tiempo, una buena manera de visualizarlos es a través de un gráfico de líneas. En el siguiente ejemplo vamos a utilizar una tabla de datos correspondiente a un sub-muestreo de un experimento RNAseq, que contiene las read counts de 8 transcritos después de 0, 3, 6, 12 y 24 horas de exposición a un estímulo X.

  • Descarga la tabla de datos VD_table2.tbl AQUÍ.

  • Primero vamos a cargar dicha tabla a la memoria de R:

tr_counts <- read_table2(file = "~/Dropbox/CastroLab_database/workshops_data/VD_table2.tbl", col_names = TRUE)
## Parsed with column specification:
## cols(
##   transcripts = col_character(),
##   `0hrs` = col_double(),
##   `3hrs` = col_double(),
##   `6hrs` = col_double(),
##   `12hrs` = col_double(),
##   `24hrs` = col_double()
## )
# Usamos la función read_table2() para leer el archivo `data/VD_table2.tbl` y asignarlo al nuevo objeto `tr_counts`.
# El argumento col_names = TRUE es para indicar la presencia de títulos o nombres para cada columna.
tr_counts # dale un vistazo a la nueva tabla de datos
## # A tibble: 8 x 6
##   transcripts  `0hrs` `3hrs` `6hrs` `12hrs` `24hrs`
##   <chr>         <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
## 1 Transcript_1  2.67   24.5   59.2     9.00    14.2
## 2 Transcript_2  0      87.0    5.71    7.61     0  
## 3 Transcript_3  1.84   54.9  131.     63.5     15.6
## 4 Transcript_4  0.256  16.6   27.1     1.33    18.5
## 5 Transcript_5  0       0      0       0        0  
## 6 Transcript_6  0       6.31   0       0        0  
## 7 Transcript_7  0       6.64   0       0        0  
## 8 Transcript_8  0     767.    92.7   141.       0

Generalmente, a través de la metodología de RNAseq, podemos obtener una tabla como tr_counts, que informa de las read counts de cada transcrito por tiempo muestral. Sin embargo, para crear un gráfico lineal, usando ggplot() + geom_line(), en el que el eje X represente una serie de tiempo, es necesario que los tiempos muestreales sean una variable (no nombres de columnas, como en tr_counts).

  • Entonces, recordamos lo aprendido anteriormente y utilizamos la función gather() para remodelar tr_counts, reuniendo los nombres de las columnas (tiempos) y transformándolos en un set de variables:
tr_plot <- tr_counts %>% 
  gather(key = time, value = read_counts, -transcripts)
# View(tr_plot)
  • Ahora podemos graficar las read counts de cada transcrito a través de los 5 puntos muestreales, usando ggplot() + geom_line():
ggplot(data = tr_plot, aes(x = time, y = read_counts, group = transcripts)) + 
  geom_line()

  • ¿Aún no? Podemos anotar las líneas para identificar cada transcrito:
ggplot(data = tr_plot, aes(x = time, y = read_counts, group = transcripts)) + 
  geom_line(aes(color = transcripts))

Como ya habrás notado, tanto en la gráfica como en la tabla, el rango de los valores de read counts es bastante amplio…

min(tr_plot$read_counts) # 'objeto$columna'
## [1] 0
max(tr_plot$read_counts)
## [1] 766.8096

Al graficar, en muchos cachos es conveniente considerar la escala de las observaciones, con el fin de obtener una mejor visualización de la distribución de los datos.

  • Vamos a convertir los valores de read counts al logarítmo base 10 de las read counts, de modo de reducir el rango de distribución de los valores. Para ello, modificamos la escala del eje Y usando scale_y_log10():
ggplot(data = tr_plot, aes(x = time, y = read_counts, group = transcripts)) + 
  geom_line(aes(color = transcripts)) + 
  scale_y_log10()
## Warning: Transformation introduced infinite values in continuous y-axis

¡Warning message! Obtenemos un gráfico, pero algo no está bien…

¿A qué crees que se debe el mensaje de advertencia? ¿Qué otro aspecto del gráfico crees tú es necesario modificar?

  • Solución:
ggplot(data = tr_plot, aes(x = time, y = log10(read_counts+1), group = transcripts)) + 
  geom_line(aes(color = transcripts))

¡Truco! Transforma la columna “time” en un factor ordenado, para que ggplot respete el orden deseado de tu serie de tiempo:

tr_plot$time <- factor(tr_plot$time, levels=unique(tr_plot$time))
ggplot(data = tr_plot, aes(x = time, y = log10(read_counts+1), group = transcripts)) + 
  geom_line(aes(color = transcripts))

3.5 Faceting

Faceting es una propiedad de ggplot que permite dividir un gráfico en múltiples gráficos basado en un factor incluído en el set de datos.

  • Continuando con el gráfico que hemos estado trabajando, usamos facet_wrap() para hacer un gráfico por cada transcrito:
ggplot(data = tr_plot, aes(x = time, y = log10(read_counts+1), group = transcripts)) + 
  geom_line(aes(color = transcripts)) + 
  facet_wrap(~transcripts)

3.6 ggplot2 themes

Themes permite configurar la apariencia o características estéticas de los gráficos creados con ggplot. Puedes revisar la lista completa de themes disponibles aquí.

  • Por ejemplo, gráficos con fondo blanco lucen mejor cuando son impresos. Podemos hacer que el fondo de las gráficas sea blanco usando la función theme_bw(), también podemos remover la cuadrícula:
ggplot(data = tr_plot, aes(x = time, y = log10(read_counts+1), group = transcripts)) + 
  geom_line(aes(color = transcripts)) + 
  facet_wrap(~ transcripts) + 
  theme_bw() + 
  theme(panel.grid = element_blank())

Mientras que, la función facet_wrap() organiza los gráficos en un número arbitrario de filas y columnas, también existe la función facet_grid que permite especificar como organizar los gráficos usando la siguiente nomenclatura: filas ~ columnas, también podemos usar . para indicar sólo una fila o columna (e.g., todas las gráficas en una columna: facet_grid(transcripts ~ .)).

3.7 Personalización

ggplot2 cuenta con varias funciones para personalizar gráficas, dale un vistazo al Cheat Sheet “Data Visualization” para pensar en formas de mejorar tus gráficos, también puedes inspirarte buscando ejemplos en internet, como aquí.

  • Vamos a revisar un ejemplo de como personalizar el boxplot que estuvimos trabajando anteriormente:
ggplot(data = data, aes(x = geo_loc_name, y = observed, color = geo_loc_name)) + 
  geom_boxplot(alpha = 0.5) +
  geom_jitter(alpha = 0.5) + 
  coord_cartesian(ylim = c(10, 85)) + 
  scale_y_continuous(breaks = c(10,20,30,40,50,60,70,80)) + 
  theme(axis.text.x = element_text(face = "bold", angle = 45, hjust = 1, size = 10), 
        axis.title.x = element_blank(), 
        axis.text.y = element_text(size = 9), 
        axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 10, b = 0, l = 0)), 
        legend.position = "none") + 
  labs(y = "Observed number of taxa")

3.8 Organizar y exportar gráficos

Anteriormente usamos faceting para separar un gráfico en múltiples gráficos. Sin embargo, también podemos crear una figura con múltiples gráficos diferentes (i.e., diferentes variables, set de datos, tipo de gráfica). El paquete gridExtra nos permite combinar gráficos separados en una única figura, usando la función grid.arrange().

p1 <- ggplot(data = data, aes(x = geo_loc_name, y = observed, color = geo_loc_name)) + 
  geom_boxplot(alpha = 0.5) +
  geom_jitter(alpha = 0.5) + 
  coord_cartesian(ylim = c(10, 85)) + 
  scale_y_continuous(breaks = c(10,20,30,40,50,60,70,80)) + 
  theme(axis.text.x = element_text(face = "bold", angle = 45, hjust = 1, size = 10), 
        axis.title.x = element_blank(), 
        axis.text.y = element_text(size = 9), 
        axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 10, b = 0, l = 0)), 
        legend.position = "none") + 
  labs(y = "Observed number of taxa")

p2 <- ggplot(data = data, aes(x = shannon, y = evenness_camargo, color = species, shape = geo_loc_name)) + 
  geom_point(alpha = 0.7, size = 4) + 
  coord_cartesian(xlim = c(0.0,3.2), ylim = c(0.005,0.085)) + 
  theme(axis.text.x = element_text(size = 9), 
        axis.title.x = element_text(size = 10, margin = margin(t = 10, r = 0, b = 0, l = 0)), 
        axis.text.y = element_text(size = 9), 
        axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 10, b = 0, l = 0)), 
        legend.text = element_text(size = 10), 
        legend.title = element_blank(), 
        legend.position = c(.02, .98),
        legend.justification = c("left", "top"),
        legend.box.just = "left",
        legend.margin = margin(2, 2, 2, 2)) + 
  labs(x = "Shannon", y = "Camargo evenness")

grid.arrange(p1, p2, ncol = 2, widths = c(1,2.5))

Después de crear tu gráfico o figura, puedes guardarlo en un archivo de formato de preferencia (e.g., png, pdf):

  • En RStudio, dirígete a la pestaña Plots y haz clic en Export para guardar el gráfico que estés visualizando.

  • También puedes usar la función ggsave() para, no solo guardar tus gráficos, sino que también modificar la dimensión y resolución por medio de los argumentos width, height y dpi.

# Asigna el gráfico a un 'objeto'
plot <- grid.arrange(p1, p2, ncol = 2, widths = c(1,2.5), heights = c(3,1))
# Exporta la nueva figura en fig_output/
ggsave("example_plot.pdf", plot, width = 14, height = 8, device = "pdf")

Ir a la página de inicio del curso