Necesidad

Aunque la implementación y ejecución de los experimentos de paneles celulares está ampliamente validada y automatizada, los centros de investigación se encuentran ahora con limitaciones muy importantes a la hora de aplicar análisis que permitan la identificación de biomarcadores. Las principales razones son:

• Integración de grandes volúmenes de datos: una de las características de las disciplinas Ómicas es que trabajan con grandes volúmenes de datos, por ejemplo, el genoma de un solo paciente tiene aproximadamente 20,000 genes y su proteoma unos 111,451 transcritos. Integrar todos estos datos con los generados en el laboratorio sin comprometer su integridad requiere conocimiento sobre cómo extraer información de bases de datos.

• Complejidad del análisis: identificar de forma veraz relaciones de causalidad entre eventos requiere del uso de métodos estadísticos y de Data Mining avanzados que deben de ser adaptados de forma correcta en cada uno de los análisis. Este hecho se vuelve todavía más crítico si tiene en cuenta la extrema complejidad de los datos biológicos (miles de variables, información no relevante, comparaciones múltiples…).

• Visualización de datos multidimensionales: la complejidad intrínseca de los datos biológicos requiere de métodos de visualización específicos, generalmente interactivos, que no suelen estar disponibles en herramientas gráficas básicas.

• Software: las herramientas disponibles para los análisis bioinformáticos son muy variadas, existen herramientas comerciales muy optimizadas pero con un alto coste por licencia (GeneData, Ingenuity-Pathway-Analysis QiaGen), o herramientas de open source de acceso gratuito pero poco unificadas y cuya curva de aprendizaje suele ser más compleja.

Todos los conocimientos necesarios para solventar con éxito los problemas descritos no suelen formar parte de las habilidades de los investigadores del ámbito de la biomedicina (biólogos, biotecnólogos, médicos), lo que, en la práctica, suele conllevar que el análisis sea delegado a expertos en otras áreas (bioinformáticos y estadísticos en su mayoría). Esto supone un riesgo para el avance de los proyectos de ámbito biomédico. En primer lugar, el análisis de la información se convierte en un cuello de botella, los datos generados por muchos investigadores deben de ser analizados por unos pocos, lo que ralentiza la obtención de resultados. En segundo lugar, y más importante, si la trasferencia de conocimiento entre profesionales de las distintas áreas no es suficiente o el analista no tiene formación sobre el ámbito del que proceden los datos, difícilmente se podrá extraer la información adecuada.

El objetivo de este documento es hacer accesibles a través del lenguaje de programación R algunos los análisis necesarios para la identificación de biomarcadores:

• Carga de datos

• Integración con datos Ómicos (información genómica y transcriptómica)

• Contrastes de hipótesis

• Análisis de correlación

• Análisis de pathways

• Clustering de compuestos

• Visualizar los resultados.
  
  

Biomarcadores

Varias décadas de tratamientos contra el cáncer han dejado patente que existe una alta heterogeneidad en cuanto a la efectividad que tiene un determinado fármaco oncológico en la población de pacientes. Los avances en el estudio de la biología molecular del cáncer, han sacado a la luz claras evidencias de que la respuesta de un fármaco depende en gran medida de factores genómicos propios del paciente y del tumor.

Un biomarcador se define como una característica biológica que puede cuantificarse y emplearse como indicador de la presencia o ausencia de una patología, o de la respuesta que se tendrá frente a un determinado tratamiento. Algunos ejemplos incluyen: patrones de expresión genética, niveles de proteína en sangre, mutaciones, etc. Son, por lo tanto, una herramienta para caracterizar y segmentar la población de pacientes de cara a posibles tratamientos (Figura 1).

Fig.1 Esquema sobre la estratificación de los pacientes. Imagen obtenida de 2017 The University of Texas MD Anderson Cancer Center.

Dentro del ámbito de los fármacos de oncología, son muchas las ventajas que implica disponer de un biomarcador, algunas de las principales son: identificación certera de los pacientes que pueden beneficiarse de un tratamiento, reducción de efectos secundarios y aumento del éxito en las fases clínicas. Un claro ejemplo es el fármaco Imatinib, que consigue un porcentaje de supervivencia en los 5 primeros años del 90% en aquellos enfermos de Leucemia Milode Crónica que tienen una aberración en los genes BCR-ABL (Druker BJ, et al. Five-year follow-up of patients receiving imatinib for chronic myeloid leukemia).

Aunque el anterior es un caso muy prometedor y que pone de manifiesto el potencial médico que tiene la identificación de biomarcadores asociados con enfermedades, son pocos los casos descubiertos. De hecho, la gran mayoría de tratamientos oncológicos capaces de combatir con notable eficacia el cáncer, no han sido asociados a una alteración biológica concreta que pueda ser empleada para seleccionar a aquellos pacientes en los que el tratamiento resulta efectivo. La combinación de la información generada en las diferentes disciplinas Ómicas (genómica, trascriptómica, proteómica, etc) junto con las metodologías de Data Mining se han convertido en la principal estrategia para la identificación de biomarcadores.

Cultivo celular

El término cultivo celular hace referencia al proceso por el cual células vivas se mantienen y se reproducen bajo condiciones controladas fuera de su ambiente natural. Dada la característica de las células cancerígenas para multiplicarse de forma descontrolada, es posible, a partir de unas pocas células extraídas en una biopsia médica, mantenerlas y expandirlas de forma ilimitada dentro de dispositivos diseñados para reproducir las condiciones del cuerpo humano (temperatura, humedad, nutrientes, gases…) (Figura 2).

Fig.2. Esquema de la obtención de cultivos celulares a partir de biopsias de pacientes. Imagen obtenida de Personalized In Vitro and In Vivo Cancer Models to Guide Precision Medicine, 2017 American Association for Cancer Research.

Esta forma de proceder ha permitido crear bancos de células en los que se recogen multitud tipos de tumorales (pulmón, intestino, colon, etc), haciendo posible que los investigadores dispongan de material biológico con el que experimentar independientemente de los pacientes. Un ejemplo de ello es la American Type Culture Collection (ATCC), una organización sin ánimo de lucro que reúne, almacena y distribuye hasta 3000 líneas celulares de origen animal.

Actividad in vitro

Estudiar la actividad in vitro de un fármaco consiste en cuantificar su actividad a través de experimentos empíricos con el objetivo de estimar el efecto que tiene sobre un determinado sistema. Estos estudios se engloban dentro de las fases pre-clínicas del desarrollo y su éxito es clave para maximizar los resultados que consiguen los fármacos cuando llegan a los pacientes (fases clínicas).

En el ámbito de los fármacos antitumorales, uno de los estudios más frecuentes de actividad in vitro consiste en medir la relación que existe entre la dosis del fármaco y la respuesta celular que provoca, lo que se conoce como cell-based drug rensponse. Para ello, las celular tumorales se exponen a distintas concentraciones del fármaco en estudio y, tras un determinado tiempo de exposición, se cuenta el número de células vivas. A continuación, se normaliza el número de células vivas respecto a un control máximo (células no expuestas al fármaco, señal máxima) y un control mínimo (células expuestas a un fármaco de referencia que las mata, señal mínima) convirtiendo así las cuentas en porcentajes de actividad.

\[\text{% actividad} = \frac{ \text{señal} - \text{señal}_{ \text{mínima}}}{\text{señal}_{\text{máxima}} - \text{señal}_{\text{mínima}}}\]

Finalmente, empleando los % de actividad, se ajusta una curva sigmoidea que representa la relación entre la dosis del fármaco y la respuesta conseguida.

\[\text{respuesta} = \text{límite inferior} + \frac{\text{límite superior - límite inferior}}{1 + (\frac{X}{IC50})^p}\]

donde los límites inferior y superior son las asíntotas de la curva, X la concentración del fármaco, p la pendiente de la curva e IC50 la concentración del fármaco con el que se consigue un 50% de la actividad máxima.

A partir de este modelo se pueden obtener múltiples métricas que describen la actividad del fármaco, algunas de las más empleadas son:

• IC50: concentración del fármaco con el que se consigue un 50% de la actividad máxima observada.

• Emax: actividad máxima observada, sea o no a la concentración máxima estudiada.

• AUC: área bajo la curva dosis-respuesta.
  
  

Curva Dosis Respuesta con R

Desde el punto de vista matemático/estadístico, una curva dosis respuesta es un modelo de regresión en el que la variable independiente (predictor) es la concentración y la variable dependiente (variable respuesta) es la actividad o efecto del fármaco. Existen varios paquetes en R que permiten ajustar curvas dosis respuesta, uno de ellos en drc, cuyos autores, además de crear el paquete, han publicado un artículo donde se pueden encontrar información muy detallada sobre el análisis de curvas dosis respuesta. A continuación, se muestran algunos ejemplos.

Supóngase un experimento de laboratorio en el que se exponen células cancerígenas a distintas concentraciones de un fármaco y se registra la actividad antiproliferativa. Se asume que los valores de actividad han sido normalizados respecto a una referencia máxima y una mínima.

library(ggplot2)

concentracion <- c(0.0108, 0.0488, 0.2195, 0.9877, 2.643, 4.4444, 20.0, 50.0)
actividad     <- c(-1.6941, -5.6772, 8.2225, 24.2046, 72.145, 98.0494, 106.7267, 104.568)
datos_actividad <- data.frame(concentracion, actividad)
head(datos_actividad)

ggplot(data = datos_actividad, aes(x = concentracion, y = actividad)) +
  geom_point() + 
  labs(title = "Actividad vs concentración") +
  theme_bw()

Con frecuencia, las concentraciones empleadas en los experimentos de screening son diluciones 1:10, por lo que la representación gráfica mejora si se emplea el logaritmo.

datos_actividad$log_conc <- log10(concentracion)
ggplot(data = datos_actividad, aes(x = log_conc, y = actividad)) +
  geom_point() + 
  labs(title = "Actividad vs logaritmo concentración") +
  theme_bw()

Con la función drm() se obtiene el ajuste de la curva, entre sus argumentos destacan:

• formula: descripción de las variables que forman el modelo en forma de ‘response ~ dose’.

• curveid: vector numérico o factor que actúa como identificador para diferenciar entre varias curvas (en caso de que las haya).

• data: dataframe que contiene los datos.

• fc: tipo de función empleada para crear el modelo (curva). La función drm() permite generar curvas empleando varios tipos de modelos, LL.4 y LL.5 para regresión logística de 4 y 5 parámetros respectivamente, y W1.4 para Weibull. Cada una de estas funciones recibe como argumento una lista con el valor de los parámetros. Si se desea que el modelo encuentre el valor óptimo de los parámetros se les da el valor NA.

• na.action: función que trate los valores ausentes, por defecto se emplea na.omit.

library(drc)

# Ajuste de la curva dosis respuesta con un modelo logístico de 4 parámetros en
# el que no se fija ninguno de ellos.
curve_fit <- drm(formula = actividad ~ concentracion, data = datos_actividad, 
                 na.action = na.omit,
                 fct = LL.4(fixed = c(NA,NA,NA,NA),
                            names = c("Hill","Bottom","Top","IC50")))

El summary del objeto devuelto por drc muestra el valor estimado de cada uno de los parámetros.

summary(curve_fit)

## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##                     Estimate Std. Error t-value   p-value    
## Hill:(Intercept)    -2.12297    0.36954 -5.7449 0.0045501 ** 
## Bottom:(Intercept)  -0.21912    3.30744 -0.0662 0.9503582    
## Top:(Intercept)    106.96927    3.82516 27.9646 9.728e-06 ***
## IC50:(Intercept)     1.76658    0.18704  9.4447 0.0007008 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  5.395822 (4 degrees of freedom)

Una vez obtenido el ajuste, se puede predecir la actividad esperada para nuevas concentraciones no probadas en el laboratorio, junto con su intervalo de confianza.

predict(object = curve_fit,
        newdata = data.frame(concentracion = 5),
        interval = "confidence")

## Prediction      Lower      Upper 
##   96.36089   86.94579  105.77599

Para extraer el valor de la IC50 relativa (concentración a la que el compuesto alcanza el 50% de su actividad máxima), o absoluta (concentración a la que el fármaco consigue una actividad del 50%) se emplea la función ED().

ED(object = curve_fit, respLev = 50, type = "relative")

## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50  1.76658    0.18704

Aunque se pude aplicar la función plot() a un objeto drc, suele ser preferible recurrir a la librería ggplot2.

# Se predice un grid de concentraciones entre la concentración mínima y máxima.
# Con estos valores se representa la curva del modelo. Cuantos más puntos se
# interpolen, mejor resolución tiene la curva.

grid_concentraciones <- seq(from = min(datos_actividad$concentracion),
                            to = max(datos_actividad$concentracion),
                            length.out = 1000)
grid_concentraciones <- data.frame(grid_concentraciones)

predicciones <- predict(object = curve_fit,
                        newdata = grid_concentraciones,
                        interval = "confidence")

# Se añade al dataframe predicciones el valor de las concentraciones.
predicciones <- cbind(grid_concentraciones, predicciones)
colnames(predicciones)[colnames(predicciones) == "grid_concentraciones"] <-
                                                                 "concentracion"
head(predicciones)

# Combinando las observaciones iniciales con los dachos predichos, se representan
# los puntos y la curva.

ggplot(data = datos_actividad,
       aes(x = log10(concentracion), y = actividad)) +
    geom_point(size = 3, shape = 1) +
    # Se añade la curva
    geom_ribbon(data = predicciones,
                aes(x = log10(concentracion), y=Prediction, ymin=Lower, ymax=Upper),
                alpha = 0.2) +
    geom_line(data = predicciones,
              aes(x=log10(concentracion), y= Prediction),
              colour = "red", size=0.8) +
    coord_cartesian(ylim = c(-20, 120)) +
    labs(title = "Curva dosis respuesta") +
    theme_bw()

Si existen varios valores de actividad para una misma concentración, por ejemplo, porque el experimento tiene varias repeticiones, la función drc las tiene en cuenta automáticamente.

concentracion <- c(0.0108, 0.0488, 0.2195, 0.9877, 2.643, 4.4444, 20.0, 50.0,
                   0.0108, 0.0488, 0.2195, 0.9877, 2.643, 4.4444, 20.0, 50.0)
actividad     <- c(-1.6941, -5.6772, 8.2225, 24.2046, 72.145, 98.0494, 106.7267,
                   104.568, 4.771, 3.834, 13.686, 36.412, 70.215, 104.320,
                   109.035, 108.163)
datos_actividad <- data.frame(concentracion, actividad)
datos_actividad$log_conc <- log10(concentracion)

curve_fit <- drm(formula = actividad ~ concentracion, data = datos_actividad, 
                 na.action = na.omit,
                 fct = LL.4(fixed = c(NA,NA,NA,NA),
                            names = c("Hill","Bottom","Top","IC50")))

grid_concentraciones <- seq(from = min(datos_actividad$concentracion),
                            to = max(datos_actividad$concentracion),
                            length.out = 1000)
grid_concentraciones <- data.frame(grid_concentraciones)

predicciones <- predict(object = curve_fit,
                        newdata = grid_concentraciones,
                        interval = "confidence")
predicciones <- cbind(grid_concentraciones, predicciones)
colnames(predicciones)[colnames(predicciones) == "grid_concentraciones"] <-
                                                                 "concentracion"

ggplot(data = datos_actividad,
       aes(x = log10(concentracion), y = actividad)) +
    geom_point(size = 3, shape = 1) +
    geom_ribbon(data = predicciones,
                aes(x = log10(concentracion), y=Prediction, ymin=Lower, ymax=Upper),
                alpha = 0.2) +
    geom_line(data = predicciones,
              aes(x=log10(concentracion), y= Prediction),
              colour = "red", size=0.8) +
    coord_cartesian(ylim = c(-20, 120)) +
    labs(title = "Curva dosis respuesta") +
    theme_bw()

En la práctica, los resultados obtenidos en los experimentos suelen contener bastante ruido, lo que hace que los datos no se ajusten bien a una curva sigmoide, produciendo valores para algunos de los parámetros que no tienen sentido desde el punto de vista biológico, por ejemplo, que el top de la curva sea muy superior al 100% o que el bottom sea muy inferior a 0%. En estos casos se pueden fijar algunos de los parámetros de la curva.

# Ajuste de una curva dosis respuesta fijando el bottom a cero y el top a 100.
curve_fit <- drm(formula = actividad ~ concentracion, data = datos_actividad, 
                 na.action = na.omit,
                 fct = LL.4(fixed = c(NA,0,100,NA),
                            names = c("Hill","Bottom","Top","IC50")))

# Los parámetros fijados no aparecen en el summary
summary(curve_fit)

## 
## Model fitted: Log-logistic (ED50 as parameter) (2 parms)
## 
## Parameter estimates:
## 
##                  Estimate Std. Error t-value   p-value    
## Hill:(Intercept) -2.12637    0.33865  -6.279 2.025e-05 ***
## IC50:(Intercept)  1.47841    0.14288  10.347 6.112e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  7.846383 (14 degrees of freedom)

grid_concentraciones <- seq(from = min(datos_actividad$concentracion),
                            to = max(datos_actividad$concentracion),
                            length.out = 1000)
grid_concentraciones <- data.frame(grid_concentraciones)

predicciones <- predict(object = curve_fit,
                        newdata = grid_concentraciones,
                        interval = "confidence")
predicciones <- cbind(grid_concentraciones, predicciones)
colnames(predicciones)[colnames(predicciones) == "grid_concentraciones"] <-
                                                                 "concentracion"

ggplot(data = datos_actividad,
       aes(x = log10(concentracion), y = actividad)) +
    geom_point(size = 3, shape = 1) +
    geom_ribbon(data = predicciones,
                aes(x = log10(concentracion), y=Prediction, ymin=Lower, ymax=Upper),
                alpha = 0.2) +
    geom_line(data = predicciones,
              aes(x=log10(concentracion), y= Prediction),
              colour = "red", size=0.8) +
    coord_cartesian(ylim = c(-20, 120)) +
    labs(title = "Curva dosis respuesta") +
    theme_bw()

# Se hace un detach del paquete drc y MASS para que no haya conflicto de funciones.
detach("package:drc", unload=TRUE)

Actividad de fármacos

El proyecto Genomics of Drug Sensitivity in Cancer se creó como una colaboración entre los centros de investigación Wellcome Sanger Institute y el Center for Molecular Therapeutics, Massachusetts General Hospital Cancer Center. El objetivo de este proyecto fue lanzar un plataforma a gran escala que permitiera analizar la actividad antitumoral de moléculas en una cantidad de líneas celulares lo suficientemente grande para que representaran a la población de pacientes de cáncer. Como resultado, se analizó un total de 266 moléculas en más de 1000 líneas celulares cancerígenas. En el momento de redactar este documento, la versión disponible es (v17.3).

• PANCANCER_IC.csv (28.4Mb): archivo coma delimitado con información detallada la actividad de cada uno de los fármacos en cada una de las líneas celulares. Las variables empleadas en el desarrollo de este documento son:

  • Drug name: identificador único (nombre) del fármaco.

  • Cell line name: identificador único (nombre) de la línea celular en la que se ha registrado la actividad.

  • IC50: valor de actividad.

Si bien el objetivo es que los investigadores empleen sus propios datos de actividad, entre las moléculas analizadas en el proyecto Genomics of Drug Sensitivity in Cancer, se incluyen fármacos que llevan años en el mercado y cuyas asociaciones con los patrones moleculares están bien descritos. Esto supone una fuente de información muy útil para validar nuevas estrategias de análisis.

En este documento, se emplean como ejemplo los datos del fármaco AZ628, un inhibidor Raf kinase desarrollado por AstraZeneca. Su eficacia ha sido probada para la inhibición de la proliferación tumoral especialmente en líneas con mutaciones BRAF V600E. Produce arresto del ciclo celular, lo que induce la apoptosis y, finalmente, la muerte celular.

Mutaciones y expresión genética

Para la identificación de biomarcadores, es necesario disponer, además de la información sobre la actividad, información sobre las mutaciones y los niveles de expresión genética de cada una de las líneas cancerígenas en las que se ha evaluado el fármaco.

Son varios los proyectos europeos y americanos que han puesto a disposición pública este tipo de información. COSMIC (Catalogue Of Somatic Mutations In Cancer) es una base de datos que recopila información molecular sobre 1015 líneas celulares y que se ha convertido en el mayor referente en cuanto a información de esta naturaleza. COSMIC ofrece acceso a toda su información sin costes de licencia a toda la comunidad de investigadores siempre y cuando no tengan fines comerciales. En el momento de redactar este documento, la versión disponible es (release v85, 8th May 2018). Los dos data sets empleados en este proyecto son:

• Complete mutation data (408.4Mb): archivo tab delimitado con información detallada sobre las mutaciones de cada una de las líneas celulares incluidas en el proyecto Cosmic Cell Lines. De entre las múltiples columnas que contiene el data set, para los análisis que se realizan en el programa aquí presentado se emplean:

  • Gene name: identificador único (nombre) del gen.

  • Sample name: identificador único de cada una de las líneas celulares.

  • Mutation Description: tipo de mutación a nivel de aminoácidos (substitution, deletion, insertion, complex, fusion etc.)

  • Primary site: tejido del que se ha extraído el tumor.

• Gene expression (650.1Mb): archivo tab delimitado con los niveles de expresión de todos los genes en cada una de las líneas celulares incluidas en el proyecto Cosmic Cell Lines. Los niveles de expresión se han obtenido con la tecnología * Affymetrix Human Genome U219 Array*.

  • Sample name: identificador único de cada una de las líneas celulares.

  • Gene Name: identificador único (nombre) del gen.

  • Z-score: nivel de expresión normalizado.

Pathways

ConsensusPathDB-human es una base de datos creada por el Max Planck Institute for Molecular Genetics que integra información sobre los genes que codifican las proteínas que participan en cada una de las rutas de señalización (pathways). Esta información es muy importante cuando se desea conocer si las alteraciones genéticas presentes un determinado cáncer están localizadas en una ruta de señalización concreta. EL contenido de la base de datos puede descargarse libremente en un archivo texto llamado CPDB_pathways_genes.tab. En su interior se encuentra la siguiente información:

• Pathway: nombre del pathway.

• Source: fuente de original de la que se ha obtenido la información.

• hgnc_symbol_ids: nombre de los genes que participan en el pathway.
  
  

Carga de datos

Se cargan los 3 sets de datos (mutaciones, expresión y pathways) empleados a lo largo del documento. Hay que tener en cuenta que todos ellos suman aproximadamente 1 GB, por lo que pueden saturar la memoria RAM. De ser así, es aconsejable cargarlos individualmente cuando sean necesarios y eliminarlos cuando no.

library(tidyverse)

# EXPRESIÓN GENÉTICA
# ==============================================================================
# Se cargan únicamente las columnas de interés que sea más rápido.

datos_expresion <- read_delim(file = "./datos/CosmicCLP_CompleteGeneExpression.txt",
                              col_names = TRUE,
                              col_types = cols_only(SAMPLE_NAME = "c",
                                                    GENE_NAME   = "c",
                                                    Z_SCORE     = "n"),
                              delim = "\t")
datos_expresion <- datos_expresion %>% rename(sample_name = SAMPLE_NAME,
                                              gene_name   = GENE_NAME,
                                              z_score     = Z_SCORE)
head(datos_expresion)

Al tratarse de valores normalizados, aquellos valores mayores que 2 se consideran atípicamente alto (gen con alta expresión) y aquellos con un valor menor que -2 se consideran atípicamente bajos (genes con baja expresión).

# MUTACIONES
# ==============================================================================
# Se cargan únicamente las columnas de interés que sea más rápido.

datos_mutaciones <- read_delim(file = "./datos/CosmicCLP_MutantExport.txt",
                               col_names = TRUE,
                               col_types = cols_only(
                                                `Gene name`    = "c",
                                                `Sample name`  = "c",
                                                `Primary site` = "c",
                                                `Mutation Description` = "c"),
                               delim = "\t")

datos_mutaciones <- datos_mutaciones %>% rename(
                                           gene_name    = `Gene name`,
                                           sample_name  = `Sample name`,
                                           primary_site = `Primary site`,
                                           mutation     = `Mutation Description`
                                          )
# Se filtran los transcritos alternativos: sus nombres siguen el siguiente patrón.
# Ejemplo: ASTN1 -> ASTN1_ENST00000367654, todos tienen el separador "_"o ".".
library(stringr)
datos_mutaciones <- datos_mutaciones %>%
                    filter(!(str_detect(string = gene_name, pattern = "_"))) %>%
                    filter(!(str_detect(string = gene_name, pattern = "[.]")))
head(datos_mutaciones)

# PATHWAYS
# ==============================================================================
datos_pathways <- read_delim(file = "./datos/CPDB_pathways_genes.txt",
                             col_names = TRUE,
                             col_types = cols_only(
                                                pathway = "c",
                                                source  = "c",
                                                hgnc_symbol_ids = "c"
                                                ),
                               delim = "\t")

head(datos_pathways)

# ACTIVIDAD DE AZ628
# ==============================================================================
datos_actividad <- read_csv(file = "./datos/activity_AZ628_panCancer.csv")
head(datos_actividad)

Líneas celulares

Número total de líneas celulares para las que se dispone de información genética:

datos_mutaciones %>% pull(sample_name) %>% unique() %>% length()

## [1] 1020

Número total de líneas celulares por histopatología (tejido de origen):

datos_mutaciones %>% select(sample_name, primary_site) %>% unique() %>%
                     group_by(primary_site) %>% count() %>%
                     ggplot(aes(x = reorder(primary_site, n), y = n)) +
                      geom_col(fill = "gray", color = "black") + 
                      coord_flip() +
                      theme_bw() + 
                      labs(title = "Número de líneas por histopatología",
                           x = "Histopatología",
                           y = "Número líneas celulares")

Número de mutaciones

mutaciones_por_linea <- datos_mutaciones %>%
                        select(sample_name, gene_name) %>%
                        group_by(sample_name) %>%
                        count() %>%
                        arrange(desc(n))
head(mutaciones_por_linea)

ggplot(data = mutaciones_por_linea, aes(x = n)) +
  geom_density(fill = "gray") + 
  labs(title = "Número de mutaciones por línea celular") +
  theme_bw()

summary(mutaciones_por_linea$n)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0   285.0   424.5   709.5   600.5 20675.0

rm(mutaciones_por_linea)

Puede observarse que, la distribución del número de mutaciones por línea celular, es muy asimétrica. La gran mayoría de líneas acumulan en torno a 600-700 mutaciones, con algunas excepciones en las que se dispara esta cantidad.

Los valores anteriores pueden resultar sospechosos puesto que, algunas líneas, tienen un número de mutaciones superior al número de genes. Esto es posible ya que, un mismo gen, puede contener múltiples mutaciones. Véase ahora el número de genes mutados por línea celular.

genes_mut_por_linea <- datos_mutaciones %>%
                        select(sample_name, gene_name) %>%
                        unique() %>% 
                        group_by(sample_name) %>%
                        count() %>%
                        arrange(desc(n))
head(genes_mut_por_linea)

ggplot(data = genes_mut_por_linea, aes(x = n)) +
  geom_density(fill = "gray") + 
  labs(title = "Número de genes mutados por línea celular") +
  theme_bw()

summary(genes_mut_por_linea$n)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0   274.5   405.5   612.1   565.0 10295.0

rm(genes_mut_por_linea)

Número medio de mutaciones por tejido (un gen mismo gen puede tener múltiples mutaciones):

mutaciones_por_tejido <- datos_mutaciones %>%
                         group_by(primary_site, sample_name) %>%
                         count() %>%
                         ungroup() %>%
                         group_by(primary_site) %>%
                         summarise(media = mean(n)) %>%
                         arrange(desc(media))

ggplot(data = mutaciones_por_tejido,
       aes(x = reorder(primary_site, media), y = media)) +
  geom_col(fill = "gray", color = "black") + 
  coord_flip() +
  theme_bw() + 
  labs(title = "Número medio de mutaciones por tejido",
       x = "Tejido")

summary(mutaciones_por_tejido$media)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   311.7   386.1   479.3   664.1   629.7  2697.9

rm(mutaciones_por_tejido)

datos_mutaciones %>%
  group_by(primary_site, sample_name) %>%
  count() %>%
  group_by(primary_site) %>%
  mutate(mediana = median(n)) %>% 
  ungroup() %>% 
  ggplot(aes(x = reorder(primary_site, mediana),  y = n, color = primary_site)) + 
    coord_flip() +
    geom_boxplot() +
    labs(title = "Número de mutaciones por tejido",
         x = "Tejido") + 
    theme_bw() +
    theme(legend.position = "none")