3.8: Dividir, Aplicar, Combinar

Aunque tocamos brevemente las funciones de escritura, la discusión hasta ahora ha sido en gran parte sobre los diversos tipos de datos utilizados por R y la sintaxis y funciones para trabajar con ellos. Estos incluyen vectores de diversos tipos, listas, marcos de datos, factores, indexación, reemplazo y reciclaje de vectores.

¿En algún momento no discutiremos analizar un conjunto de datos significativo? Sí, y ese punto es ahora. Continuaremos trabajando con los datos de expresión génica del capítulo 32, “Carácter y datos categóricos”, después de que se hayan procesado para incluir las columnas individuales apropiadas usando las funciones de la biblioteca stringr. Aquí está el código del último capítulo que pre-procesó los datos de entrada:

A medida que continuemos, es probable que queramos volver a ejecutar nuestro script de análisis de diferentes maneras, sin tener que volver a ejecutar todo el código anterior, ya que tarda unos minutos en ejecutar este script de preprocesamiento. Una buena estrategia sería escribir este marco de datos “limpiado” en un archivo de texto usando write.table (), y volver a leerlo en un script de solo análisis con read.table (). Sin embargo, solo por variedad, usemos la función save (), que almacena cualquier objeto R (como una lista, vector o marco de datos) en un archivo en un formato binario comprimido.

La extensión de archivo tradicional para dicho archivo de objeto R es .Rdata. Dichos archivos de datos a menudo se comparten entre los programadores de R.

Las primeras líneas del marco de datos guardado son todas lecturas para el mismo gen. Hay ~360.000 líneas en el marco de datos completo, con lecturas para más de 11 mil genes. (El script anterior, que procesa el archivo expr_long_coded.txt, se puede encontrar en el archivo Expr_Preprocess.r.) La siguiente figura muestra los resultados de print (head (expr_long_split)):

Comencemos un nuevo script de análisis que cargue este marco de datos antes de comenzar nuestro análisis del mismo. Para este script también podemos necesitar la biblioteca stringr, y vamos a incluir también la biblioteca dplyr (que suponemos que se ha instalado por separado con install.packages (“dplyr”) en la consola interactiva).

La función load () crea una variable con el mismo nombre que se utilizó en save (), y por conveniencia hemos creado una copia del marco de datos con un nombre de variable más corto, expr. Si bien este tipo de archivos de datos R son convenientes para los usuarios de R, el beneficio de un enfoque write.table () sería que otros usuarios e idiomas tendrían acceso más directo a los datos (por ejemplo, en Python o incluso Microsoft Excel).

Este conjunto de datos representa un análisis multifactorial de la expresión génica ^[1] en dos genotipos de una planta, donde se han aplicado a cada uno un tratamiento control (agua) y un tratamiento químico (pesticida). Además, se probaron tres tipos de tejido (A, B y C, para hoja, tallo y raíz, respectivamente) y se probaron dos o tres repeticiones para cada combinación para determinar el poder estadístico. Estos datos representan lecturas de micromatrices y ya se han normalizado y por lo tanto están listos para su análisis. ^[2] (Estos datos también son anonimizados para esta publicación a petición de los investigadores.)

Por ahora, nos centraremos en las variables de genotipo y tratamiento. Ignorando otras variables, esto nos da alrededor de 24 lecturas de expresión para cada gen.

Este tipo de diseño es bastante común en los estudios de expresión génica. No es sorprendente que la expresión génica difiera entre dos genotipos, y que la expresión génica difiera entre dos tratamientos. Pero los genotipos y tratamientos reaccionan frecuentemente de formas interesantes. En este caso, quizás se sabe que el genotipo L4 crece más grande que el genotipo C6, y se sabe que el tratamiento químico es dañino para las plantas. Curiosamente, parece que el tratamiento químico es menos dañino para el genotipo L4.

La pregunta natural es: ¿qué genes podrían estar contribuyendo a la mayor durabilidad de L4 bajo exposición química? No podemos simplemente buscar las diferencias en la expresión media entre genotipos (con, digamos, t.test ()), ya que eso capturaría muchos genes involucrados en el tamaño de la planta. Tampoco podemos buscar aquellos con diferencias de expresión medias que sean grandes entre tratamientos, porque sabemos que muchos genes reaccionarán a un tratamiento químico pase lo que pase. Incluso ambas señales juntas (lado izquierdo de la imagen de abajo) no son tan interesantes, ya que representan genes que se ven afectados por el químico y están involucrados con el tamaño de la planta, pero pueden no tener nada que ver con características únicas de L4 en el tratamiento químico. Lo que realmente nos interesa son los genes con una respuesta diferencial a la sustancia química entre dos de los dos genotipos (lado derecho de la imagen de abajo).

¿Significa eso que podemos buscar genes donde la diferencia media sea grande en un genotipo, pero no en el otro? ¿O dentro de un tratamiento, pero no ese otro? No, porque muchos de estos interesantes escenarios se perderían:

Para las mediciones de expresión de cada gen, necesitamos un modelo clásico de ANOVA (análisis de varianza). Después de contabilizar la varianza explicada por la diferencia de genotipo (a la que podemos asignar un valor de p que indique si existe una diferencia) así como la varianza explicada por el tratamiento (que también tendrá un valor p), es la variación sobrante—para la interacción/ respuesta diferencial—significativa?

Una prueba de un solo gen, fórmulas

Para comenzar, ejecutaremos un ANOVA para un solo gen. ^[3] Primero, necesitaremos aislar un submarco de datos que consiste en lecturas para un solo ID. Podemos usar la función unique () para obtener un vector de identificadores únicos, y luego el operador %in% para seleccionar filas que coincidan solo con la primera de ellas.

Este marco de datos contiene filas para un solo gen y columnas para expresión, genotipo, tratamiento y otros (podríamos verificar esto con una simple print () si es necesario). Se desea ejecutar un ANOVA para estas variables, basado en un modelo lineal que prediga valores de expresión a partir del genotipo, el tratamiento y la interacción de estos términos.

El siguiente paso es construir un modelo lineal simple usando la función lm (), que toma como único parámetro requerido una “fórmula” que describa cómo deben relacionarse los valores predictores con los valores predichos. Observe el formato de la fórmula: response_vector ~ predictor_vector1 + predictor_vector2 + predictor_vector1: predictor_vector2, indicando que deseamos determinar cómo predecir valores en el vector de respuesta usando elementos de los vectores predictores.

Si los vectores o factores de interés existen dentro de un marco de datos (como es el caso aquí), entonces la función lm () puede funcionar solo con nombres de columna y un argumento data =.

Entonces podemos imprimir el modelo lineal devuelto con print (lm1), así como su estructura con str (lm1), para ver que contiene bastante información como una lista con nombre:

Antes de continuar, ¿cuál es esta expresión ~ genotipo + tratamiento + genotipo: tratamiento “fórmula”? Como es habitual, podemos investigar un poco más asignándola a una variable y usando class () y str ().

El resultado indica que la clase es de tipo “formula” y que tiene un atributo “.Environment”. Es un poco opaco, pero un tipo de fórmula en R es un contenedor para un vector de caracteres de nombres de variables y una sintaxis para relacionarlos entre sí. El atributo environment especifica dónde deben buscarse esas variables por defecto, aunque las funciones son libres de ignorarlas (como en el caso de data = for lm (), que especifica que las variables deben considerarse nombres de columna en el marco de datos dado). Las fórmulas ni siquiera necesitan especificar variables que realmente existen. Considere la siguiente fórmula y la función all.vars (), que inspecciona una fórmula y devuelve un vector de caracteres de los nombres de variables únicos que aparecen en la fórmula.

En fin, volvamos a ejecutar la prueba estadística para el gen único que hemos aislado. El resultado lm1 contiene el modelo que predice expresiones de los otros términos (donde este “modelo” es una lista con un conjunto particular de elementos). Para obtener los valores p asociados, tendremos que ejecutarlo a través de la función anova () de R.

Al imprimir el resultado se muestra que los valores de p (etiquetados “Pr (>F)”) para el genotipo, el tratamiento y los términos de interacción son 0.97, 0.41 y 0.56, respectivamente. Si queremos extraer los valores p individualmente, primero tendremos que inspeccionar su estructura con str (), revelando que el resultado es tanto una lista como un marco de datos, no es sorprendente porque un marco de datos es un tipo de lista. Los tres valores p se almacenan en el nombre “Pr (>F)” de la lista.

Podemos así extraer el vector de valores p ya que pvals1 <- anova1 The resulting list will have the same three values as above, and the names of the elements in the output list will be inherited from the input list (so sample_irqs" title="Rendido por Quicklatex.com" height="2836" width="2358" style="vertical-align: -4px;"> s1 contendrá el primer rango o 5.44, para esta muestra aleatoria de todos modos).

Este es un ejemplo de la poderosa estrategia de “split-apply-combine” para el análisis de datos. En general, esta estrategia implica dividir un conjunto de datos en pedazos, aplicar una operación o función a cada uno, y luego combinar los resultados en un único conjunto de salida de alguna manera.

Cuando la entrada y la salida son ambas listas, esta operación también se conoce como “mapa”. De hecho, este paradigma subyace a algunos frameworks de supercomputadoras de procesamiento paralelo como MapReduce de Google y la versión de código abierto de Hadoop (aunque estos programas no están escritos en R).

Debido a que las listas son un tipo de datos tan flexible, lapply () es bastante útil. Aún así, en algunos casos, como este, preferiríamos que la salida fuera un vector en lugar de una lista. Esto es lo suficientemente común como para garantizar una función de conversión llamada unlist () que extrae todos los elementos y subelementos de una lista para producir un vector.

Ser advertido: debido a que las listas son más flexibles que los vectores, si la lista dada a unlist () no es simple, los elementos pueden estar presentes en el vector en un orden impar y serán coaccionados a un tipo de datos común.

La función lapply () también puede tomar como entrada un vector en lugar de una lista, y no es más que una de las muchas funciones “apply” en R. Otros ejemplos incluyen apply (), que aplica una función a cada fila o columna de una matriz, ^[4] y sapply () (para “simple apply”), que detecta si el tipo de retorno debe ser un vector o matriz dependiendo de la función aplicada. Es fácil confundirse con esta diversidad, así que para comenzar, lo mejor es enfocarse en lapply () dada su naturaleza relativamente sencilla (pero útil).

Recopilación de parámetros con...

La página de ayuda para lapply () es interesante. Especifica que la función toma tres parámetros: X (la lista a aplicar), FUN (la función a aplicar), y... , que la página de ayuda detalla como “argumentos opcionales a FUN”.

El... parámetro es común en R, si un poco misterioso a veces. Este “parámetro” recoge cualquier número de parámetros nombrados suministrados a la función en una sola unidad. Esta colección puede ser utilizada entonces por la función; a menudo se pasan a funciones llamadas por la función.

Para ilustrar este parámetro, modifiquemos nuestra función ipercentile_range () para tomar dos parámetros opcionales, un percentil inferior y superior a partir de los cuales se debe calcular el rango. Si bien el rango intercuartílico se define como el rango en los datos del percentil 25 al 75, nuestra función podrá calcular el rango desde cualquier percentil inferior o superior.

Ahora podemos llamar a nuestra función como ipercentile_range (muestra to get the default interquartile range of the first sample, or as ipercentile_range(samples" title="Rendido por Quicklatex.com" height="40" width="586" style="vertical-align: -4px;"> s2, inferior = 0.05, superior = 0.95), para obtener el rango del percentil 90 medio. De igual manera, tenemos la capacidad de suministrar los parámetros opcionales a través del... de lapply (), que a su vez los suministrará a cada llamada de la función ipercentile_range durante el paso de aplicación.

Esta sintaxis puede parecer un poco incómoda, pero algo así es necesario si esperamos pasar funciones que toman múltiples parámetros como parámetros a otras funciones. Esto también debería ayudar a iluminar una de las características comúnmente incomprendidas de R.

Las funciones están en todas partes

Aunque en su mayoría ocultas a la vista, la mayoría (si no todas) las operaciones en R son en realidad llamadas a funciones. A veces los nombres de las funciones pueden contener caracteres no alfanuméricos. En estos casos, podemos encerrar el nombre de la función en backticks, como en iqr1 <- `ipercentile_range` (samples$s1). (En el capítulo 27, “Variables y datos”, aprendimos que cualquier nombre de variable se puede incluir en backticks de esta manera, aunque rara vez se recomienda).

¿Y algo así como el operador de suma humilde, +? Funciona como una función: dadas dos entradas (vectores del lado izquierdo y del lado derecho), devuelve una salida (la suma elemento por elemento).

De hecho, + realmente es una función que toma dos parámetros, pero para llamarlo así necesitamos usar backticks.

Esto podría considerarse algo sorprendente: la versión “legible por humanos” es convertida en una llamada de función correspondiente detrás de escena por el intérprete. Esta última, la representación funcional es más fácil de usar para el intérprete que la representación amigable con el ser humano. Este es un ejemplo de azúcar sintáctico: ajustes a un lenguaje que lo hacen “más dulce” para los programadores humanos. R contiene bastante azúcar sintáctica. ¡Incluso el operador de asignación <- es una llamada de función! La línea a <- c (1, 2, 3) es en realidad azúcar sintáctica para `<-` (“a”, c (1, 2, 3)).

Debido a que las funciones son solo un tipo de datos y los nombres de las funciones son variables, podemos reasignar fácilmente nuevas funciones a nombres antiguos. Quizás una compañera ha dejado abierta su sesión R, así que decidimos redefinir + para significar “al poder de”:

Ok, tal vez eso sea un poco demasiado malo. Otra característica de R es que los nombres de función que comienzan y terminan con% y toman dos parámetros señales de que la función puede ser utilizada como una función de infijo vía azúcar sintáctica.

Marco de datos azucarado Split-Apply-Combinar

Aunque las unidades básicas de datos en R son vectores, y las listas se utilizan para almacenar datos heterogéneos, los marcos de datos son importantes como mecanismo de almacenamiento de datos tabulares. Debido a que lapply () puede aplicar una función a cada elemento de una lista, también puede aplicar una función a cada columna de un marco de datos. Este hecho se puede utilizar, por ejemplo, para extraer todas las columnas numéricas de un marco de datos con una aplicación de lapply () y is.numeric () (este último devuelve TRUE cuando su entrada es de tipo numérico; dejaremos los detalles como un ejercicio).

Para la mayoría de los datos tabulares, sin embargo, no nos interesa aplicar independientemente una función a cada columna. Más bien, más a menudo deseamos aplicar una función a diferentes conjuntos de filas agrupadas por una o más de las columnas. (Esto es exactamente lo que queremos hacer con la función sub_df_to_pvals () que escribimos para nuestro análisis de expresión génica). El paquete dplyr (install.packages (“dplyr”)) proporciona esta habilidad y es potente y fácil de usar, una vez que nos acostumbramos a su azúcar sintáctica especializada.

Inicialmente, cubriremos dos funciones en el paquete dplyr, group_by () y do (): group_by () agrega metadatos (como atributos) a un marco de datos indicando qué columnas categóricas definen grupos de filas y devuelve dicho marco de datos “agrupado”, mientras que do () aplica una función dada a cada grupo de un marco de datos agrupado y devuelve un marco de datos de resultados con información de agrupación incluida.

Para ilustrar, primero crearemos un marco de datos simple en el que trabajar, representando muestras de peces de uno de dos lagos diferentes: Green Lake o Detroit Lake.

La impresión muy bien formateada:

Primero, agruparemos el marco de datos por la columna de especies usando group_by (). Esta función toma un marco de datos (o matriz con nombre de columna o similar) y una lista de nombres de columna (sin comillas) que definen los grupos como parámetros adicionales. Al igual que con el uso del parámetro data = para lm (), la función busca dentro del marco de datos las columnas especificadas.

¿Este marco de datos es diferente del original? Sí, pero sobre todo en los atributos/metadatos. Cuando se imprime, obtenemos cierta información sobre la “fuente” (donde se almacenan los datos, dplyr puede acceder a datos de bases de datos remotas) y los grupos:

Prácticamente, a diferencia de los marcos de datos regulares, los marcos de datos “agrupados” imprimen solo las primeras filas y columnas, incluso si tienen muchos miles de filas de largo. A veces esta es una característica indeseable; afortunadamente, ejecutar data.frame () en un marco de datos agrupado devuelve una versión desagrupada. La clase () de un marco de datos agrupado devuelve “data.frame” así como “tbl”, “tbl_df” y “grouped_df”. Porque un marco de datos agrupado también es un marco de datos regular (¡y también es una lista!) , todavía podemos hacer todas las operaciones de indexación elegantes sobre ellas cubiertas en capítulos anteriores.

La función do () aplica una función a cada grupo de filas de un marco de datos agrupado usando una estrategia split-apply-combine. La función aplicada debe tomar como primer parámetro un marco de datos, y debe devolver un marco de datos. Vamos a escribir una función que, dado un marco de datos o subtrama de datos (grupo) de los datos de peces, devuelve un marco de datos de una sola fila con columnas para mean_weight y sd_weight. Tanto las funciones mean () como sd () pueden tomar un argumento na.rm para eliminar cualquier valor NA potencial en sus entradas, así que quizás nuestra función debería tomar un parámetro opcional similar.

Esta función toma como entrada un marco de datos con una columna para el peso y devuelve un marco de datos con estadísticas de resumen. Podemos ejecutarlo en el marco de datos original:

Alternativamente, podemos llamar a do () en el marco de datos agrupado, diciéndole que ejecute mean_sd_weight () en cada grupo (sub-trama de datos). La sintaxis para do () difiere ligeramente de la de lapply () en que especificamos a. para el argumento posicional que representa el submarco de datos.

El resultado del do () es otro marco de datos agrupado, compuesto por las filas devueltas por la función aplicada. Observe que se han agregado las columnas de agrupación, aunque no las especificamos en el ret_df dentro de la función.

Al desarrollar funciones que funcionan con do (), es posible que se encuentre con un error como Error: Los resultados no son marcos de datos en las posiciones: 1, 2. Este mensaje de error indica que la función no está devolviendo un tipo de trama de datos, que se requiere para do (). Para diagnosticar problemas como este, puede agregar sentencias print () dentro de la función para inspeccionar el contenido de las variables a medida que se aplica la función.

También podemos agrupar marcos de datos por múltiples columnas, lo que resulta en un solo grupo por combinación de entradas de columna. ^[5]

Los valores de NA para algunas desviaciones estándar son el resultado de llamar sd () en un vector de longitud uno (porque solo hubo una medición de trucha por lago).

Aunque la función aplicada debe tomar un marco de datos y devolver un marco de datos, no hay restricciones sobre la naturaleza del marco de datos devuelto. Aquí nuestra función devuelve un marco de datos de una sola fila, pero podría devolver varias filas que serían cosidas juntas en el paso de combinación. Como ejemplo, aquí hay una función que, dada una trama de datos, calcula la media () de la columna de peso, y resta esa media de todas las entradas, devolviendo la trama de datos modificada (la llamada “normalización media” de los datos).

¡Y entonces podemos fácilmente normalizar los datos por grupo!

En la salida anterior, -1.15 y 1.15 son las desviaciones de la media del grupo de truchas, y las otras son desviaciones de la media para el grupo bajo.

Más Azúcar, Parámetros Opcionales, Resumir

Algo a tener en cuenta sobre el call to do () es que difiere sintácticamente de la llamada a lapply (). En do (), especificamos no solo la función a aplicar a cada grupo, sino también cómo se llamará esa función, usando. para denotar el marco de datos del grupo de entrada. Esto es algo más claro cuando queremos especificar argumentos opcionales a la función aplicada. En este caso, es posible que queramos especificar que los valores NA deben eliminarse estableciendo remove_nas = TRUE en cada llamada a mean_sd_weight ().

Hablando de azúcar sintáctico, el paquete magrittr (que se instala y carga junto con dplyr, aunque escrito por un autor diferente) proporciona una interesante función de infijo, %>%. Considere el patrón común anterior; después de la creación de un marco de datos (fish), lo ejecutamos a través de una función para crear un resultado intermediario (fish_by_species de group_by ()), ejecutarlo a través de otra función para obtener otro resultado (stats_by_species de mean_sd_weight ()), y así sucesivamente. Si bien este proceso es claro, podríamos evitar la creación de múltiples líneas de código si estuviéramos dispuestos a anidar algunas de las llamadas a funciones.

Esta línea de código ciertamente sufre algunos problemas de legibilidad. La función de infijo %>% suministra su lado izquierdo a su lado derecho (ubicado por. en el lado derecho), devolviendo el resultado de la llamada a la función derecha. Así podemos reescribir lo anterior de la siguiente manera.

Observe la similitud entre esto y el encadenamiento de métodos en Python y tuberías stdin/stdout en la línea de comandos. (También es posible omitir el. si sería el primer parámetro, como en peces %>% group_by (species) %>% do (mean_sd_weight (.)) .)

Demostramos que las funciones aplicadas por do () pueden devolver marcos de datos de varias filas (en el ejemplo de normalización media), pero nuestra función mean_sd_weight () solo devuelve un marco de datos de una sola fila con columnas hechas de estadísticas de resumen simples. Para este tipo de necesidad simplificada, el paquete dplyr proporciona una función especializada llamada summary (), tomando un marco de datos agrupado y otros parámetros que indican cómo se deben calcular dichas estadísticas de resumen grupales. Este ejemplo produce la misma salida que el anterior, sin la necesidad de la función intermediaria mean_sd_weight ().

El paquete dplyr proporciona bastantes otras características y funciones para explorar, pero incluso la simple combinación de group_by () y do () representan un paradigma poderoso. Debido a que R es un lenguaje tan flexible, es probable que el tipo de azúcar sintáctico avanzado que usa dplyr se vuelva más común. Aunque se requiere un cierto esfuerzo para comprender estos lenguajes específicos de dominio (DSL) y cómo se relacionan con R “normal”, el tiempo empleado suele valer la pena.

Ejercicios

1. El propósito principal de lapply () es ejecutar una función en cada elemento de una lista y cotejar los resultados en una lista. Con algo de creatividad, también se puede utilizar para ejecutar una llamada de función idéntica una gran cantidad de veces. En este ejercicio veremos cómo se distribuyen los valores p bajo el modelo “nulo”, cuando se comparan dos muestras aleatorias de la misma distribución.

   Comienza escribiendo una función null_pval () que genere dos muestras aleatorias con rnorm (100, media = 0, sd = 1), las compara con t.test (), y devuelve el valor p. La función debe tomar un solo parámetro, digamos, x, pero en realidad no hacer ningún uso de ella.

   A continuación, genere una lista de 10,000 números con 10k_nums_list <- as.list (seq (1,10000)), y llame a 10k_pvals_list <- lapply (10k_nums, null_pval). Convierte esto en un vector con 10k_pvals_vec <- unlist (10k_pvals_list) e inspecciona la distribución con hist (10k_pvals_vec).

   ¿Qué revela esta prueba? ¿Qué hace el código que escribiste y por qué funciona? ¿Por qué la función necesita tomar un parámetro x que ni siquiera se usa? ¿Qué sucede si cambias una de las muestras aleatorias para usar rnorm (100, media = 0.1, sd = 1)?

2. Si df es un marco de datos, el uso de la indexación [] lo trata como una lista de elementos (columnas). Escribe una función llamada numeric_cols () que tome un marco de datos como parámetro y devuelva una versión del marco de datos con solo las columnas numéricas guardadas. Es posible que desee hacer uso de lapply (), unlist (), como.numeric (), y el hecho de que las listas (y los marcos de datos cuando se tratan como listas) se pueden indexar por vector lógico.

   Como ejemplo, si df1 <- data.frame (id = c (“PRQ”, “XL2", “BB4"), val = c (23, 45.6, 62)), entonces print (numeric_cols (df1)) debe imprimir un marco de datos con solo la columna val. Si df2 <- data.frame (srn = c (461, 514), name = c (“Mel”, “Ben”), age = c (27, 24)), entonces print (numeric_cols (df2)) debe imprimir un marco de datos con solo columnas srn y age.

3. Escribe una función llamada subset_rows () que tome dos parámetros: primero, un marco de datos df, y segundo, un entero n. La función debe devolver un marco de datos que consiste en n filas aleatorias de df (puede encontrar la función sample () de uso).

   El marco de datos del iris es un conjunto de datos R incorporado comúnmente referenciado, que describe mediciones de pétalos para una variedad de especies de iris. Aquí está la salida de print (head (iris)):
   Use group_by () para agrupar el marco de datos del iris por la columna Species, y luego do () junto con su función subset_rows () para generar un marco de datos que consta de 10 random hileras por especie.

4. Curiosamente, no hay nada que detenga una función llamada por do () (o lapply ()) de sí misma llamando a do () (y/o lapply ()). Escriba un conjunto de funciones que compare cada especie en el marco de datos del iris con todas las demás especies, reportando la diferencia media en Pétal.Ancho. (En esta impresión la diferencia de medias compara especiesA con especiesB.)
   Para lograr esto, querrá escribir una función compare_petal_widths () que tome dos sub-marcos de datos, uno que contenga datos para la especie A (sub_dfa) y el otro para la especie B (sub_dfb), y devuelva un marco de datos que contenga el nombre de la A especie, el nombre de la especie B y la diferencia de la media de pétalos. Ancho. Puede probar su función extrayendo manualmente marcos de datos que representan setosa y versicolor.

   A continuación, escribe una función one_vs_all_by_species () que nuevamente tome dos parámetros; el primero será un marco de subdatos que represente una sola especie (sub_df), pero el segundo será un marco de datos con datos para todas las especies (all_df). Esta función debe agrupar all_df por especie, y usar do () para llamar compare_petal_widths () en cada grupo enviando a lo largo de sub_df como parámetro secundario. Se debe devolver el resultado.

   Finalmente, una función all_vs_all_by_species () puede tomar una sola trama de datos df, agruparla por especie y llamar a one_vs_all_by_species () en cada grupo, enviando a lo largo de df como parámetro secundario, devolviendo el resultado. Al final, todo lo que se necesita es llamar a esta función en el iris.

Expresión Génica, Terminada

Con la discusión de split-apply-combine y dplyr en nuestro haber, volvamos a la tarea de crear y analizar un modelo lineal para cada ID en el conjunto de datos de expresión génica. Como recordatorio, habíamos dejado de haber leído en el conjunto de datos “limpiados”, extraer un marco de subdatos que representaba un único ID y escribir una función que toma dicho subcuadro de datos y devuelve un marco de datos de una sola fila de valores p. (Ahora debería quedar claro por qué pasamos por la molestia de asegurarnos de que nuestra función tome un marco de datos como parámetro y devuelva uno también).

Ahora, podemos usar group_by () en el marco de datos expr para agrupar por la columna id, y do () para aplicar la función sub_df_to_pvals_df () a cada grupo. Sin embargo, en lugar de trabajar en todo el conjunto de datos, creemos un expr10 para contener un marco de datos que represente mediciones para 10 IDs; si estamos satisfechos con los resultados, siempre podemos analizar la tabla expr completa (aunque el conjunto de datos completo tarda solo un par de minutos en analizarlo) .

El resultado es una tabla bien organizada de valores p para cada gen en el conjunto de datos:

Hay un tema más importante a considerar para un análisis como este: la corrección de pruebas múltiples. Supongamos por un momento que ninguno de los ~11.000 genes se expresa diferencialmente de ninguna manera. Debido a que los valores de p se distribuyen uniformemente entre cero y uno bajo la hipótesis nula (sin diferencia), solo para la columna de genotipos podríamos esperar ~11,000 * 0.05 = 550 aparentemente (pero erróneamente) diferencias significativas. Alterando el escenario, supongamos que hubo alrededor de 50 genes que realmente se expresaron diferencialmente en el experimento: estos probablemente tendrían valores de p pequeños, pero es difícil separarlos de los otros 550 valores aparentemente significativos.

Hay una variedad de soluciones a este problema. Primero, podríamos considerar un umbral de significación mucho más pequeño. ¿Qué umbral debería ser? Bueno, si quisiéramos que el número de valores de p aparentemente —pero no realmente— significativos fuera pequeño (en promedio), digamos, 0.05, podríamos resolver para 11,000 * α = 0.05, sugiriendo un corte de α = 0.000004545. Esto se conoce como corrección Bonferroni. El problema es que esta es una corrección conservadora, y es probable que incluso nuestros 50 genes significativos no pasen ese pequeño umbral.

La razón por la que la corrección de Bonferroni es tan conservadora es que hemos especificado que sólo aceptaremos que el número promedio de falsos positivos sea 0.05, lo cual es bastante pequeño. Alternativamente podríamos decir que aceptaríamos 10 falsos positivos (y resolveríamos por 11,000 * α = 10); si el número de resultados guardados es de 100, eso es solo una tasa de descubrimiento falso (FDR) del 10%, y así podemos esperar que el 90% de los resultados guardados sean reales (¡pero no sabremos cuáles son cuáles!). Por otro lado, si el número de resultados guardados termina siendo de 15, esa es una tasa de falsos descubrimientos del 75%, lo que indica que podemos esperar que la mayoría de los resultados guardados sean falsos positivos.

En lugar de especificar el número de falsos positivos que estamos dispuestos a aceptar, en su lugar podemos especificar la tasa y tratar con cualquier número de resultados guardados que salgan. Hay varios de estos métodos de “control de FDR” disponibles, y algunos hacen más suposiciones sobre los datos que otros. Por ejemplo, los métodos de Benjamini y Hochberg (descritos en 1995 y conocidos por las siglas “BH”) y Benjamini, Hochberg y Yekutieli (de 2001, conocidos como “BY”) difieren en la cantidad de independencia que se supone que tienen las muchas pruebas. El método BY asume menos independencia entre las pruebas, pero generalmente da como resultado más resultados para una tasa de FDR dada. (Consulte a su estadístico local o libro de texto de estadísticas para obtener más detalles.)

La función p.fit () nos permite ejecutar los métodos de corrección Bonferroni, BH o BY. Se necesitan dos argumentos: primero, un vector de valores p para ajustar, y, segundo, un método = parámetro que se puede establecer en “bonferroni”, “BH” o “BY”. (Hay otras posibilidades, pero estas tres son las más utilizadas). Se devuelve un vector de la misma longitud de valores de FDR ajustados; seleccionar todas las entradas con valores menores que Q equivale a establecer un límite de tasa de FDR de Q.

En el análisis final de todos los genes, produciremos una columna By-ajustada para los valores de interacción p, y luego seleccionaremos del conjunto de datos solo aquellas filas donde la tasa de FDR sea menor a 0.05. (El script de análisis completo se encuentra en el archivo Expr_Analyze.r.)

Desafortunadamente (o afortunadamente, dependiendo de la cantidad de datos que se espera analizar más a fondo), solo se identifica que un gen tiene una interacción significativa después de la corrección BY en este análisis. (Si se devolvieran 100, esperaríamos que 5 de ellos fueran falsos positivos debido al corte de 0.05 FDR utilizado).

Para los marcos de datos agrupados, print () omitirá la impresión de algunas columnas si no caben en la ventana de visualización. Para ver todas las columnas (incluyendo nuestra nueva columna Interaction_by), podemos convertir la tabla agrupada de nuevo en un marco de datos regular e imprimir la cabeza () de la misma con print (head (data.frame (pvals_interaction_sig))). Para los curiosos, el valor Interaction_by de este ID es 0.048.

Ejercicios

1. Gastamos una buena cantidad de trabajo asegurando que el marco de datos generado y devuelto por sub_df_to_pvals_df () se generó mediante programación. Aprovecha este hecho haciendo de la fórmula un parámetro secundario de la función, es decir, como función (sub_df, form), de manera que dentro de la función el modelo lineal se pueda producir como lm1 <- lm (form, data = sub_df).

   A continuación, ejecute el análisis con group_by () y do (), pero especifique alguna otra fórmula en la llamada do () (tal vez algo así como expresión ~ genotipo + tratamiento + tejido).

2. En un ejercicio anterior de este capítulo, escribimos una función que tomó una parte del marco de datos de CO2 y devolvió un marco de datos con un valor p para una comparación de los valores de absorción de CO2 entre plantas refrigeradas y no refrigeradas. Ahora, use group_by () y do () para ejecutar esta prueba para cada grupo conc, produciendo un marco de datos de valores p. Use p.fit () con corrección de Bonferroni para agregar una columna adicional de valores corregidos múltiples.
3. Si bien las funciones llamadas por do () deben tomar un marco de datos como su primer parámetro y deben devolver un marco de datos, también son libres de realizar cualquier otra acción, como generar una trama o escribir un archivo. Utilice group_by (), do () y write.table () para escribir el contenido del marco de datos de CO2 en siete archivos de texto, uno por cada nivel conc. Denles el nombre conc_95.txt, conc_175.txt, y así sucesivamente. Es posible que necesite paste () o str_c () (de la biblioteca stringr) para generar los nombres de archivo.

---

1. Aquí usaremos el término coloquial “gen” para referirnos a un locus genético que produce ARN mensajero y su expresión para ser representada por la cantidad de ARNm (de las diversas isoformas) presente en la célula en el momento del muestreo.
   
2. Las micromatrices son una tecnología más antigua para medir la expresión génica que han sido suplantadas en gran medida por los métodos directos de RNA-seq. Para ambas tecnologías, los niveles generales de expresión varían entre muestras (por ejemplo, debido a la eficiencia de la medición) y deben modificarse para que los valores de diferentes muestras sean comparables; los datos de micromatrices a menudo también se transforman en preparación para el análisis estadístico. Ambos pasos se han realizado para estos datos; hemos evitado la discusión de estos pasos preparatorios ya que son específicos de la tecnología utilizada y un área de investigación activa para RNA-seq. Generalmente, los análisis de RNA-seq implican trabajo tanto en la línea de comandos para el preprocesamiento de datos como en el uso de funciones de paquetes R especializados.
   
3. Hay muchas formas interesantes en las que podríamos analizar estos datos. Por el bien de enfocarnos en R en lugar de métodos estadísticos, nos apegaremos a un modelo ANOVA 2x2 por gen relativamente simple.
   
4. Debido a que la función apply () opera sobre matrices, convertirá silenciosamente cualquier marco de datos que se le haya dado en una matriz. Debido a que las matrices pueden contener solo un tipo (mientras que los marcos de datos pueden contener diferentes tipos en diferentes columnas), usar apply () en un marco de datos primero forzará una autoconversión para que todas las columnas sean del mismo tipo. Como tal, casi nunca es correcto usar apply () en un marco de datos con columnas de diferentes tipos.
   
5. Debido a que group_by () y funciones similares toman nombres de columna “desnudos” o “sin comillas”, pueden ser difíciles de llamar dado un vector de caracteres. Para funciones en los paquetes dplyr y tidyr (cubiertas más adelante), las versiones con _ como parte del nombre de la función soportan esto. Por ejemplo, fish_by_species_lake <- group_by (fish, species, lake) puede ser reemplazado por fish_by_species_lake <- group_by_ (fish, c (“species”, “lake”)).