Predicción de la supervivencia de Adenocarcinoma de pulmón mediante la integración de perfiles de mutación y anotación clínica

Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile

1. Introducción

El cáncer es una enfermedad genética compleja que surge de los efectos combinados de múltiples variaciones genéticas y epigenéticas (Zhang et al. 2016). Estos cambios pueden influir en la vulnerabilidad celular, provocando la transformación de una célula normal a una cancerosa anormal (Dimitrakopoulos and Beerenwinkel 2017). Estudios recientes han determinado que los cánceres son entidades altamente mutadas, donde los genes afectados están implicados principalmente en la proliferación, regeneración y apoptosis celular, impulsando el crecimiento maligno dentro de los tumores (Kou et al. 2016).
Con el advenimiento de Next Generation Sequencing (NGS) y el desarrollo de la Genómica se ha revolucionado el entendimiento de la biología del cáncer, mejorando el diagnóstico y las terapias, dando paso a la era de la medicina de precisión y la genómica del cáncer (Garraway and Lander 2013).
La medicina de precisión es definida como un enfoque clínico que busca seleccionar el método terapéutico más adecuado para cada paciente considerando diferentes parámetros clínicos y biomarcadores en amplios campos clínicos, como por ejemplo enfermedades metabólicas y cardiovasculares. Ahora, debido a la incipiente investigación de la genómica del cáncer, la oncología clínica será uno de los campos más beneficiados por este nuevo enfoque de investigación que integra el estudio del genoma (Kou et al. 2016). Específicamente, la genómica del cáncer involucra estudios sistemáticos del genoma humano para la identificación de alteraciones genéticas recurrentes en tipos específicos de cáncer. Estas alteraciones genéticas consisten en Single Nucleotide Variants (SNVs), pequeñas inserciones y deleciones (Indels), fusión de genes, Copy-Number Variations (CNVs) y grandes re-arreglos cromosomales, también llamados Structural Variants (Cheng, Zhao, and Zhao 2016).
Actualmente, los investigadores se han enfocado en la identificación de mutaciones “Drivers”, en genes asociados fuertemente al cáncer, que confieren una ventaja proliferativa selectiva a la célula con cáncer con respecto a células normales. Este tipo de mutación estaría causalmente implicada en la oncogénesis, siendo los principales blancos para el desarrollo de nuevas terapias (Nik-Zainal 2014).
Mediante estos estudios genómicos se ha podido determinar que el cáncer no es sólo complejo, sino que también es altamente heterogéneo, debido a que los mecanismos genéticos pueden variar entre pacientes del mismo tipo patológico (Yang et al. 2014). De esta forma, se ha podido identificar variantes genéticas predisponentes y subtipos tumorales, basándose en la firma molecular heterogénea de cada paciente (Riazalhosseini and Lathrop 2016).

2. Problemática

La discriminación de mutaciones “Drivers” de las mutaciones aleatorias o “Passenger” que no desempeñan un papel significativo en el desarrollo del cáncer es actualmente un desafío. Así mismo, éstas y otras mutaciones presentes en un paciente podrían estar asociadas con el diagnóstico y/o pronóstico del cáncer, confiriéndoles un valor predictivo independientemente de la histología del tumor.

3. Hipótesis

Los patrones de mutaciones en los genes Drivers del cáncer y su relación con la anotación de casos clínicos pueden ser utilizados para diseñar métodos predictivos de supervivencia de pacientes con cáncer.

4. Avances y desafíos en la genómica del cáncer

Actualmente, se han realizado varios estudios que buscan caracterizar el escenario genómico para los principales cánceres, como por ejemplo cáncer de pulmón, mama, próstata, estómago entre otros, los cuales tienen una alta incidencia y mortalidad. Esta caracterización genómica consiste principalmente en la identificación de genes driver y pathways, los cuales poseen un rol clave en la generación y desarrollo del cáncer. En el trabajo de (Li et al. 2016) se identificaron los principales genes drivers en Adenocarcinoma de Pulmón (LUAD, Lung Adenocarcinoma) en población asiática. Otro estudio en LUAD se propuso caracterizar el escenario genómico de esta enfermedad mediante la búsqueda de señales de mutación asociadas con la progresión del tumor, logrando identificar genes drivers por medio de la integración de datos epigenómicos, genómicos, la evolución clonal e información de las características clínicas.
En otro reciente estudio, se describió la prevalencia y ubicación de mutaciones somáticas y re-ordenamientos cromosómicos en muestras de tumor en pacientes con cáncer de próstata. Además, se identificaron mutaciones propias de pacientes afroamericanos, reportando una nueva fusión de genes presente en el 17% de los pacientes del estudio (Lindquist et al. 2016).
En la investigación de (Heo et al. 2017) se identificaron mutaciones somáticas en pacientes coreanos con Leucemia Mieloide Aguda (LMA, Acute Myeloid Leukemia), concluyendo que los subtipos morfológicos de la LMA pueden ser reflejados por patrones específicos de alteraciones genómicas. Además, los autores exponen la necesidad de estudios de este tipo, los cuales son muy útiles para el diagnóstico temprano de estas enfermedades complejas.
En el trabajo de (Pereira et al. 2016), no sólo identifican los 40 genes con mutaciones drivers en 2.433 muestras de tumor en cáncer primario de mama, si no que buscaron patrones de asociación entre eventos mutacionales somáticos evaluando el nivel de relación con la sobrevida y patrones clínicos.
Si bien la identificación de genes drivers es un importante avance para lograr comprender la biología del cáncer, determinar patrones de genes mutados que permitan estratificar los pacientes a nivel genómico es el siguiente paso y nuevo desafío.
Recientes estudios han planteado que la clasificación del cáncer en subtipos clínica y biológicamente significativos, puede aumentar la eficiencia de la predicción del diagnóstico y sobrevida, así como también, mejorar la estrategia terapéutica. Sin embargo, este enfoque de investigación es muy reciente, y sólo se ha explorado el genoma de los principales cánceres, dejando una amplia lista de otros tipos de cánceres para su caracterización utilizando su información genómica. Por tal motivo, el presente trabajo tiene por objetivo crear perfiles genómicos de muestras de Adenocarcinoma de Pulmón de The Cancer Genome Atlas (TCGA) para la predicción de sobrevida mediante la integración de anotación clínica y patrones mutacionales.

5. Descripción de los datos.

Los datos con los que trabajaremos corresponden a 230 registros de pacientes de Adenocarcinoma de Pulmón extraídos de Cbioportal for Cancer Genomic. Por cada paciente se tiene la cantidad de mutaciones en cada uno de los genes de su secuencia de ADN: ABL1, AKT1, AKT3, ARAF, AXL, BRAF, CCND1, CDK4, CDK6, etc. Además por cada paciente se registran una serie de datos clínicos: género, edad, tipo de cáncer, estadío del cáncer, estatus como fumador, cuantos meses sobrevivió luego del diagnóstico, entre otros.

#Cargando datos
library(data.table)
setwd("../data")
mutFreq <- fread("luad2014_mut_freq_clinical_data.txt")
caddSco <- fread("luad2014_cadd_score_clinical_data.txt")

A continuación se muestran algunas de las principales características clínicas de los datos utilizados.

library(ggplot2)
require(cowplot)
blank_theme <- theme_minimal() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.border = element_blank(),
    panel.grid=element_blank(),
    axis.ticks = element_blank(),
    axis.text.x=element_blank(),
    plot.title=element_text(size=12, face="bold")
  )
#Gender
data <- setNames(nm = c("Gender", "Freq"), as.data.frame(table(mutFreq$GENDER, useNA="ifany")))
genderPieChart <- ggplot(data, aes(x="", y=Freq, fill=Gender)) +
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start=0) +
  blank_theme +
  geom_text(aes(label = Freq), position = position_stack(vjust = 0.5), size=5) +
  ggtitle("Distribución de Género")
#Overall_Survival_Status
y<-table(mutFreq$OS_STATUS, useNA="ifany")
data <- setNames(nm = c("Overall_Survival_Status", "Freq"), as.data.frame(y))
overallSurvPieChart <- ggplot(data, aes(x="", y=Freq, fill=Overall_Survival_Status)) +
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start=0) +
  blank_theme +
  geom_text(aes(label = Freq), position = position_stack(vjust = 0.5), size=5) + 
  ggtitle("Estado de Supervivencia")
#AGE
ageHistChart <- ggplot(mutFreq, aes(x=AGE)) +
  geom_histogram(binwidth = 5, na.rm = TRUE) +
  geom_vline(xintercept = mean(mutFreq$AGE, na.rm = TRUE), col="blue") +
  ggtitle("Distribución de Edad")
#Tumor_stage_2009
data <- setNames(nm = c("Tumor_stage_2009", "Freq"), as.data.frame(table(mutFreq$TUMOR_STAGE_2009, useNA="ifany")))
tumorStagePieChart <- ggplot(data, aes(x="", y=Freq, fill=Tumor_stage_2009)) +
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start=0) +
  blank_theme +
  geom_text(aes(label = Freq), position = position_stack(vjust = 0.5), size=3) + 
  ggtitle("Estadío del tumor")
plot_grid(plot_grid(genderPieChart, overallSurvPieChart, align='hv'), plot_grid(ageHistChart, tumorStagePieChart), nrow = 2)

#TOBACCO_SMOKING_HISTORY_INDICATOR
y<-table(mutFreq$TOBACCO_SMOKING_HISTORY_INDICATOR, useNA="ifany")
data <- setNames(nm = c("TOBACCO_SMOKING_HISTORY_INDICATOR", "Freq"), as.data.frame(y))
smokingPieChart <- ggplot(data, aes(x="", y=Freq, fill = TOBACCO_SMOKING_HISTORY_INDICATOR)) +
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start=0) +
  blank_theme +
  geom_text(aes(label = Freq), position = position_stack(vjust = 0.5), size=3) +
  ggtitle("Estatus de fumador")
#Overall Survival VS TOBACCO_SMOKING_HISTORY_INDICATOR
os <- ecdf(mutFreq$OS_MONTHS)
x <- mutFreq[, c("TOBACCO_SMOKING_HISTORY_INDICATOR", "OS_MONTHS")]
x <- na.omit(x)
library(plyr)
x <- arrange(x, TOBACCO_SMOKING_HISTORY_INDICATOR, OS_MONTHS)
x.ecdf <- ddply(x, .(TOBACCO_SMOKING_HISTORY_INDICATOR), transform, ecdf = ecdf(OS_MONTHS)(OS_MONTHS) )
SurvivalCumBySmokingStepChart <- ggplot() +
  geom_step(data = mutFreq, aes(OS_MONTHS, y = 100 * (1 - os(OS_MONTHS)), colour = "All"), size = 1.2) +
  geom_step(data = x.ecdf, aes(OS_MONTHS, y = 100 * (1 - ecdf), colour = TOBACCO_SMOKING_HISTORY_INDICATOR)) +
  labs(title = "Sobrevivencia", colour = "TOBACCO_SMOKING_HISTORY", x = "Meses Sobrevividos", y = "Sobrevivencia (%)")
plot_grid(smokingPieChart, SurvivalCumBySmokingStepChart, align='hv', nrow = 2)

#Overall Survival VS Gender
x <- mutFreq[, c("GENDER", "OS_MONTHS")]
x <- na.omit(x)
library(plyr)
x <- arrange(x, GENDER, OS_MONTHS)
x.ecdf <- ddply(x, .(GENDER), transform, ecdf = ecdf(OS_MONTHS)(OS_MONTHS) )
SurvivalCumByGenderStepChart <- ggplot() +
  geom_step(data = x.ecdf, aes(OS_MONTHS, y = 100 * (1 - ecdf), colour = GENDER), size = 0.8) +
  labs(title = "Sobrevivencia", colour = "GENDER", x = "Meses Sobrevividos", y = "Sobrevivencia (%)")
SurvivalByGenderAndTumorSubtypeBarChart <- ggplot(aes(y = OS_MONTHS, x = HISTOLOGICAL_SUBTYPE, fill = GENDER), data = na.omit(mutFreq[, c("OS_MONTHS", "HISTOLOGICAL_SUBTYPE", "GENDER")])) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8), axis.text.y = element_text(size=8), legend.text = element_text(size = 8), legend.title = element_text(size = 10)) +
  xlab("Subtipo de Tumor") +
  ylab("Meses Sobrevividos")
plot_grid(SurvivalCumByGenderStepChart, SurvivalByGenderAndTumorSubtypeBarChart, align='hv', nrow = 2)

SurvivalSmokingGenderBoxplot <- ggplot(aes(y = OS_MONTHS, x = TOBACCO_SMOKING_HISTORY_INDICATOR, fill = GENDER), data = mutFreq) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
  xlab("Tipo de fumador") +
  ylab("Overall survival")
#SurvivalSmokingOSStatusBoxplot <- ggplot(aes(y = OS_MONTHS, x = TOBACCO_SMOKING_HISTORY_INDICATOR, fill = OS_STATUS), data = mutFreq) +
#  geom_boxplot() +
#  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
#  xlab("Tipo de fumador") +
#  ylab("Overall survival")
#plot_grid(SurvivalSmokingGenderBoxplot, SurvivalSmokingOSStatusBoxplot, align='v')
#SurvivalTumorStageGenderBoxplot <- ggplot(aes(y = OS_MONTHS, x = TUMOR_STAGE_2009, fill = GENDER), data = mutFreq) +
#  geom_boxplot() +
#  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
#  xlab("Estadio de tumor") +
#  ylab("Overall survival")
#SurvivalTumorStageOSStatusBoxplot <- ggplot(aes(y = OS_MONTHS, x = TUMOR_STAGE_2009, fill = OS_STATUS), data = mutFreq) +
#  geom_boxplot() +
#  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
#  xlab("Estadio de tumor") +
#  ylab("Overall survival")
#plot_grid(SurvivalTumorStageGenderBoxplot, SurvivalTumorStageOSStatusBoxplot, align='v')
SurvivalTumorSubtypeGenderBoxplot <- ggplot(aes(y = AGE, x = HISTOLOGICAL_SUBTYPE, fill = GENDER), data = mutFreq) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
  xlab("tipo de tumor") +
  ylab("Overall survival")
#SurvivalTumorSubtypeOSStatusBoxplot <- ggplot(aes(y = OS_MONTHS, x = HISTOLOGICAL_SUBTYPE, fill = OS_STATUS), data = mutFreq) +
#  geom_boxplot() +
#  theme(axis.text.x = element_text(angle = 50, hjust = 1, size=8)) +
#  xlab("tipo de tumor") +
#  ylab("Overall survival")
plot_grid(SurvivalSmokingGenderBoxplot, SurvivalTumorSubtypeGenderBoxplot, align='v')

En la siguiente tabla se muestra, por cada gen la cantidad total de mutaciones presentes en la data, cantidad de pacientes que presentan mutaciones en este gen y el porciento que representa del total.

#Get total of mutations by gene
genes <- setNames(nm = c("Gene", "Mut"), stack(colSums(mutFreq[, 23:14780]))[2:1])
#Get total of CADD Score by gene
pat <- setNames(nm = c("Gene", "Score"), stack(colSums(caddSco[, 2:14759]))[2:1])
pat <- merge(pat, genes[, c("Gene", "Mut")], by="Gene")
#Add total of patients by mutated gene
pat <- merge(pat, setNames(nm = c("Gene", "Pat"), stack(colSums(caddSco[,  2:14759] / caddSco[,  2:14759], na.rm = TRUE))[2:1]), by="Gene")
pat<- pat[order(pat$Score, decreasing = TRUE), ]
#Add CADD Score of mutated gene by patient
pat$Freq <- apply(data.table(pat[, 4]), 1, function(x){return(x * 100 / 230)})
pat

5.1 Pre-procesamiento

La mayoría de estudios genómicos se han basado en medidas o indicadores que permiten interpretar variantes genómicas y relacionarlas con la generación y progresión del cáncer, por ejemplo, la frecuencia mutacional. En varios estudios esta medida ha sido clave para la identificación de genes drivers. Sin embargo, otras investigaciones proponen utilizar diferentes medidas que correlacionen la funcionalidad molecular como la patogenicidad de las mutaciones en los genes. Estas medidas son conocidas como Measures Deleteriousness.
En el trabajo de (Vural, Wang, and Guda 2016) se utilizó un método para integrar anotación funcional, conservación e información del modelo del gen en un único score llamado el Combined Annotation Dependent Depletion (CADD) score (Kircher et al. 2014), que es una Mesasures Deleteriousness. Este score se correlaciona con la diversidad alélica, anotaciones, patogenicidad de las variantes codificantes, los efectos reguladores medidos experimentalmente y variantes causales de enfermedades dentro del genoma. Este score toma valores en el rango +- infinito, donde un valor alto indica grandes efectos deletéreos de la mutación. En este trabajo se propone crear dos tipos de perfiles basados en la frecuencia mutacional y en una medida de patogenicidad de mutaciones con el objetivo de evaluarlas como medidas predictoras de sobrevida.

5.1.1 Perfil basado en Mutation score

Este tipo de perfil corresponde a una matriz de Mutation Score (n x p), donde n es el número de muestras y p el número de genes. Cada entrada (i, j) contendrá la suma de los CADD score de las mutaciones encontradas para la muestra i en el gen j, esta suma es conocida como Mutation Score. Los Mutation Score serán transformados a una escala de 0 a infinito positivo.

5.2.2 Perfil basado en Fecuencia mutacional

A este perfil corresponde otra matriz con las mismas características estructurales (n x p), donde cada celda (i, j) contendrá la frecuencia relativa de mutaciones para la muestra i en el gen j.

6. Metodología

6.1 Perfiles mutacionales

En cuanto a los datos mutacionales, se crearán perfiles mutacionales basados en el Mutation score, descrito previamente. Se Obtendrá una matriz donde las filas corresponden a las muestras y las columnas a los genes con su correspondiente score de malignidad. Esta matriz será normalizada entre valores de 0 y 1, utilizando el método de máximos y mínimos.

En la siguiente figura se muestra el flujo de trabajo que se realizará para los datos mutacionales.

Flujo de trabajo para datos datos mutacionales

6.2 Datos clínicos

En cuando a los datos clínicos, estos pasarán por un proceso de Binarización, un método para codificar variables categóricas, donde cada categoria para cada atributo será convertido en un nuevo atributo con valores binarios (0,1) para denotar la presencia y ausencia de esa categoria para cada muestra. También, los datos númericos, en este caso la edad, será normalizada según el método de máximos y mínimos. Posteriormente, se analizará la distribución de la sobrevida, ya que para los métodos de predicción que se utlizarán se requiere que la variable respuesta, en este caso sobrevida, posea una distribución Normal. Por lo tanto, se evaluará que distribución posee esta variable, se corroborará la necesidad de eliminar outliers, si existen, los cuales pueden tener un efecto negativo en el análisis, y también se verá si es necesario aplicar una transformación (logaritmo) para que la sobrevida tenga una distribución del tipo Normal. En la siguiente figura se muestra el flujo de trabajo que se realizará para los datos clínicos

Flujo de trabajo para datos datos mutacionales

luego, las matrices resultantes de los datos mutacionales y clínicos serán unidos por el Id de las muestras, creandose una matriz o perfil clínico-mutacional. Posteriormente, esta matriz será dividida en dos set de datos, el de entrenamiento, con el 75% de los datos y el de test con el 25% restante de forma aleatoria. Los algoritmos que se aplicarán serán entrenados con el set de entrenamiento y luego se aplicarán con el set test. En este trabajo se aplicarán algoritmos de predicción y de clasificación.

Flujo de trabajo para la predicción de la sobrevida

6.3 Algoritmos de predicción

6.3.1 Principal Component Regression (PCR)

El método Principal Component Regression (PCR), a diferencia de las técnicas de regresión tradicionales que trabajan directamente sobre los datos originales, es capaz de extraer los componentes principales (PCs), que son combinaciones lineales de los features originales, realizando luego una regresión lineal con los coeficientes de los PCs obtenidos. Por lo general, sólo un subconjunto de PCs es elegido para la regresión. La selección tradicional se basa en identificar los PCs con una alta varianza y que estén altamente asociados a la variable respuesta.Se utilizará la función prcomp para calcular los componentes principales (PCs), y la función lm para ajustar el modelo lineal.

6.3.2 Partial Least Squares Regression (PLSR)

Partial Least Squares Regression (PLSR) es una alternativa supervisada a PCR. Al igual que la PCR, PLSR es un método de reducción de dimensiones, que primero identifica un nuevo conjunto más pequeño de características que son combinaciones lineales de las características originales, luego se ajusta un modelo lineal a través de mínimos cuadrados a las nuevas características M. Sin embargo, a diferencia de PCR, PLSR hace uso de la variable de respuesta para identificar las nuevas características. Se usará la función plsr del paquete pls para realizar este análisis.

6.3.3 Regression Tree

Los árboles de predicción se utilizan para predecir una respuesta o clase \(Y\). para las entradas \(X_1, X_2, X_3, ..., X_n\). Si es una respuesta continua, se llama un árbol de regresión, si es categórico, se llama árbol de clasificación. En cada nodo del árbol, se comprueba el valor de una entrada \(X_i\) y dependiendo de la respuesta (binaria) se continúa hacia la sub-rama izquierda o hacia la derecha. Cuando se llega a una hoja, se encuentra la predicción (generalmente es una estadística simple del conjunto de datos que representa la hoja, como el valor más común de las clases disponibles). En este caso, se utlizará un árbol de predicción, ya que la sobrevida es númerica. Para esto se utilizará la función rpart del paquete rpart.

7. Resultados

En cuanto sobrevida, esta presentó una distribución exponencial, por lo tanto se aplicó una transformación logarítmica. También, se encontró dos outliers los cuales fueron removidos del set de datos.

Distribución original de la sobrevida

Distribución transformada y filtrada de la sobrevida

7.2 Principal Component analysis (PCA) con regresión lineal de los PCs significativos

Los datos fueron centralizados y escalados, como requerimiento de PCA. Luego, se calculó PCA usando la función prcomp(). Los primeros 10 componentes fueron los que explican una mayor proporción de la varianza de los datos. Luego, estos 10 PCs se ingresaron como predictores en un modelo de regresión lineal donde la variable respuesta es la sobrevida. El resultado de esto, entregó 5 PCs significativos, PC1, PC2, PC6, PC7 y PC10, los cuales fueron usados como predictores en un nuevo modelo de regresión lineal. Este último modelo fue usado para la predicción con el set set.

Flujo de trabajo para PCA con regresión lineal

Para los resultados de la predicción se calculo el Mean Squares Error, el cual fue de 0.8047338.

7.3 Partial Least Squares Regression (PLSR)

Para determinar el número de componentes que se utilizó en la regresión, se estima el Mean squared error of prediction (MSEP) usando el método de crosvalidación. A continuación, en la siguiente figura se muestra el MSEP para los componentes principales, donde en el 100PC hay un punto de inflección.

MSEP para los PCs

Posteriormente, se realizó la predicción usando el número de PCs con mejor MSEP y el set test. El MSE de la predicción con el set de datos fue de 6.183472.

7.4 Regression Tree

En la siguiente figura se muestra el modelo de regression tree creado con el set de entrenamiento. Se puede observar que la mayoría de los nodos corresponden a variables clínicas, y que los genes que allí se encuentran no son en su mayoría genes frecuentemente mutados. El valor de MSE de la predicción con el set test fue de 2.447874.

7.5 Algoritmos de Clasificación

La Sobrevida es una variable continua, pero para tratar el problema como un problema de clasificación se hizo necesario transformarla a una variable categórica. Esto se hizo creando 6 clases que cubrieran distintos rangos de tiempo en meses que representan la distribución de la Sobrevida. Se intentó en primer lugar que las clases estuvieran balanceadas, y luego que reflejaran rangos de tiempos representativos para el problema. A continuación se muestra una descripción de esta nueva variable categórica (OS_CLASS).

data <- caddSco
data$OS_CLASS <- as.factor(cut(data$OS_MONTHS, c(0, 1, 6, 12, 24, 36, 300), include.lowest = TRUE))
summary(data$OS_CLASS)

   [0,1]    (1,6]   (6,12]  (12,24]  (24,36] (36,300]     NA's 
      32       37       28       37       29       40       27 

qplot(x = OS_CLASS, data = data)

Para tratar el problema se usaron los siguientes algoritmos de clasificación:

• Dummy Classifier (Base)
  
• Logistic Regression (LR)
  
• Linear Discriminant Analysis (LDA)
  
• KNNeighbors Classifier (KNN)
  
• Decision Tree Classifier (CART)
  
• Gaussian Naive Bayes (NB)
  
• Support Vector Machine (SVM)

La experimentación se realizó usando solamente los genes “drivers”, es decir los 53 genes que se consideran más fuertemente asociados al desarrollo del cáncer. Para comprobar los resultados obtenidos, se realizó además un test estadístico para ver si los clasificadores eran significativamente diferentes, obteniéndose

# Load libraries
import pandas
from pandas.tools.plotting import scatter_matrix

import matplotlib.pyplot as plt

from sklearn import model_selection
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score

from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import SVC
from sklearn.dummy import DummyClassifier

from statsmodels.stats.weightstats import ttest_ind


import numpy as np
from sklearn.model_selection import train_test_split

## run_classifier recibe un clasificador, un dataset (X, y) 
## y opcionalmente la cantidad de resultados que se quiere obtener del clasificador
def run_classifier(clf, X, y, num_tests=10):
    scores = []
    
    for _ in range(num_tests):
        ### COMPLETAR ACÁ
        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)
        clf.fit(X_train, y_train)
        scores.append( clf.score(X_test, y_test) )  # X_test y y_test deben ser definidos previamente
    
    return np.array(scores)
    


## run_classifiers recibe un dataset (X, y)
def run_classifiers(X, y):
    # Spot Check Algorithms
    models = []
    models.append(('Base', DummyClassifier(strategy='stratified')))
    models.append(('LR', LogisticRegression()))
    #models.append(('LDA', LinearDiscriminantAnalysis()))
    models.append(('KNN', KNeighborsClassifier()))
    models.append(('CART', DecisionTreeClassifier()))
    models.append(('NB', GaussianNB()))
    models.append(('SVM', SVC()))

    result_list = []
    table_head = "NAMES" + "\t" + "Accuracy"
    for name, clf in models:
        table_head += "\t" + name
        accuracys = run_classifier(clf, X, y, 40)
        result_list.append((name, accuracys))

    print("+ indica diferencia significativa\n")
    print(table_head)
    for name1, results1 in result_list:
        table_row = name1 + "\t" + str(results1.mean())
        #print("Comparando %s - Accuracy: %.2f" % (name1, results1.mean()))
        for name2, results2 in result_list:
            table_row = table_row + "\t"
            if name1 == name2:
                continue

            _, p_value, __ = ttest_ind(results1, results2)  # t-test: test estadistico
        
            if p_value <= 0.05:
                sig = "+"
            else:
                sig = ""
            table_row = table_row + sig
            #print("%s:\t%.2f %s" % (name2, results2.mean(), sig))
        print(table_row)




# Load clinical dataset
dataset = pandas.read_csv("../data/clinical_data_bin_with_OS_Class_AGE_NORM.csv")


#Clinical Data
data = dataset.drop(['OS_MONTHS'], axis = 1).drop(['AGE_NORM'], axis = 1)
X = data.values[:, 1:(data.columns.size-1) ]
y = data.values[:, data.columns.size-1]

print("Resultados de Clasificadores con Datos Clínicos\n")
run_classifiers(X, y)
print()
print()



# Load CADD 53 driver genes data
dataset_gen_num_53 = pandas.read_csv("../data/gene_num_53.csv")
#dataset_gen_cadd = pandas.read_csv("../data/gene_cadd_norm.csv")


# CADD 53 Driver Genes Data
#print(dataset.loc[:, ['Row.names', 'OS_CLASS']])
data = pandas.merge(dataset_gen_num_53, dataset.loc[:, ['Row.names', 'OS_CLASS']], on = 'Row.names')
X = data.values[:, 1:(data.columns.size-1)]
y = data.values[:, data.columns.size-1]

print("Resultados de Clasificadores con Datos Genéticos (CADD)\n")
run_classifiers(X, y)
print()
print()


dataset_gen_num_53 = pandas.read_csv("../data/gene_num_53.csv")
#dataset_gen_cadd = pandas.read_csv("../data/gene_cadd_norm.csv")

data = pandas.merge(dataset_gen_num_53, dataset, on = 'Row.names').drop(['OS_MONTHS'], axis = 1).drop(['AGE_NORM'], axis = 1)
X = data.values[:, 1:(data.columns.size-1)]
y = data.values[:, data.columns.size-1]

print("Resultados de Clasificadores con Datos Clínicos y Genéticos\n")
run_classifiers(X, y)

Resultados de Clasificadores con Datos Cl攼㹤nicos

+ indica diferencia significativa

NAMES   Accuracy    Base    LR  KNN CART    NB  SVM
Base    0.166666666667      +       +   +   
LR  0.257352941176  +       +   +   +   +
KNN 0.185294117647      +               
CART    0.196078431373  +   +               +
NB  0.207843137255  +   +               +
SVM 0.166176470588      +       +   +   


Resultados de Clasificadores con Datos Gen攼㸹ticos (CADD)

+ indica diferencia significativa

NAMES   Accuracy    Base    LR  KNN CART    NB  SVM
Base    0.169607843137                      
LR  0.158823529412                      
KNN 0.16862745098                       
CART    0.164215686275                      
NB  0.150980392157                      
SVM 0.167156862745                      


Resultados de Clasificadores con Datos Cl攼㹤nicos y Gen攼㸹ticos

+ indica diferencia significativa

NAMES   Accuracy    Base    LR  KNN CART    NB  SVM
Base    0.169117647059      +       +       +
LR  0.247058823529  +       +   +   +   +
KNN 0.166176470588      +       +       
CART    0.21568627451   +   +   +       +   +
NB  0.157352941176      +       +       
SVM 0.150490196078  +   +       +       

Los resultados obtenidos no son muy buenos, destacándose los clasificadores Logistic Regression y Decision Tree como los de mejores resultados. Se hace necesario un análisis más detallado de las causas de estos resultados.

8. Conclusiones y trabajo futuro

• En los métodos de predicción, los perfiles clínico-mutacional tienen un menor MSE.
  
• En los métodos de clasificación, los perfiles clínico y clínico-mutacional mostraron un mejor accuracy, aunque no se obtuvieron muy buenos resultados en general.
  
• Usar otros set de datos con mayor cantidad de muestras.
  
• Detectar y explorar la estructura interna de los datos con técnicas de Minería de Datos (otros algoritmos de agrupamiento, etc.).
  
• Experimentación con otras técnicas de reducción de dimensionalidad, y otros algoritmos de regresión y de clasificación.
  
• Aplicar otras métricas para comparar los resultados de las predicciones.
  
• Identificar los componentes que contribuyen en la sobrevida.

9. Bibliografía

Cheng, Feixiong, Junfei Zhao, and Zhongming Zhao. 2016. “Advances in Computational Approaches for Prioritizing Driver Mutations and Significantly Mutated Genes in Cancer Genomes.” Briefings in Bioinformatics 17 (4): 642-56. doi:10.1093/bib/bbv068.
Dimitrakopoulos, Christos M., and Niko Beerenwinkel. 2017. “Computational Approaches for the Identification of Cancer Genes and Pathways.” Wiley Interdisciplinary Reviews. Systems Biology and Medicine 9 (1). doi:10.1002/wsbm.1364.
Garraway, Levi A., and Eric S. Lander. 2013. “Lessons from the Cancer Genome.” Cell 153 (1): 17-37. doi:10.1016/j.cell.2013.03.002.
Heo, Seong Gu, Youngil Koh, Jong Kwang Kim, Jongsun Jung, Hyung-Lae Kim, Sung-Soo Yoon, and Ji Wan Park. 2017. “Identification of Somatic Mutations Using Whole-Exome Sequencing in Korean Patients with Acute Myeloid Leukemia.” BMC Medical Genetics 18 (March). doi:10.1186/s12881-017-0382-y.
Kircher, Martin, Daniela M Witten, Preti Jain, Brian J O’Roak, Gregory M Cooper, and Jay Shendure. 2014. “A General Framework for Estimating the Relative Pathogenicity of Human Genetic Variants.” Nature Genetics 46 (3): 310-15. doi:10.1038/ng.2892.
Kou, Tadayuki, Masashi Kanai, Shigemi Matsumoto, Yasushi Okuno, and Manabu Muto. 2016. “The Possibility of Clinical Sequencing in the Management of Cancer.” Japanese Journal of Clinical Oncology 46 (5): 399-406. doi:10.1093/jjco/hyw018.
Li, Shiyong, Yoon-La Choi, Zhuolin Gong, Xiao Liu, Maruja Lira, Zhengyan Kan, Ensel Oh, et al. 2016. “Comprehensive Characterization of Oncogenic Drivers in Asian Lung Adenocarcinoma.” Journal of Thoracic Oncology 11 (12): 2129-40. doi:10.1016/j.jtho.2016.08.142.
Lindquist, Karla J., Pamela L. Paris, Thomas J. Hoffmann, Niall J. Cardin, Rémi Kazma, Joel A. Mefford, Jeffrey P. Simko, et al. 2016. “Mutational Landscape of Aggressive Prostate Tumors in African American Men.” Cancer Research 76 (7): 1860-68. doi:10.1158/0008-5472.CAN-15-1787.
Nik-Zainal, Serena. 2014. “Insights into Cancer Biology through next-Generation Sequencing.” Clinical Medicine 14 (Suppl 6): s71-77. doi:10.7861/clinmedicine.14-6-s71.
Pereira, Bernard, Suet-Feung Chin, Oscar M. Rueda, Hans-Kristian Moen Vollan, Elena Provenzano, Helen A. Bardwell, Michelle Pugh, et al. 2016. “The Somatic Mutation Profiles of 2,433 Breast Cancers Refine Their Genomic and Transcriptomic Landscapes.” Nature Communications 7 (May): 11479. doi:10.1038/ncomms11479.
Riazalhosseini, Yasser, and Mark Lathrop. 2016. “Precision Medicine from the Renal Cancer Genome.” Nature Reviews Nephrology 12. doi:10.1038/nrneph.2016.133.
Vural, Suleyman, Xiaosheng Wang, and Chittibabu Guda. 2016. “Classification of Breast Cancer Patients Using Somatic Mutation Profiles and Machine Learning Approaches.” BMC Systems Biology 10 (Suppl 3). doi:10.1186/s12918-016-0306-z.
Yang, William, Kenji Yoshigoe, Xiang Qin, Jun S Liu, Jack Y Yang, Andrzej Niemierko, Youping Deng, et al. 2014. “Identification of Genes and Pathways Involved in Kidney Renal Clear Cell Carcinoma.” BMC Bioinformatics 15 (Suppl 17): S2. doi:10.1186/1471-2105-15-S17-S2.
Zhang, Fan, Chunyan Ren, Kwun Kit Lau, Zihan Zheng, Geming Lu, Zhengzi Yi, Yongzhong Zhao, et al. 2016. “A Network Medicine Approach to Build a Comprehensive Atlas for the Prognosis of Human Cancer.” Briefings in Bioinformatics 17 (6): 1044-59. doi:10.1093/bib/bbw076.