• Fenotipo: SM300-Efficiency
  • Modelo: Random Forest
  • Significancia: No
  • Datos: 54 SNP’s
  • Finalidad: Realiza un modelo por cada gen (54 en total)

Preprocesamiento de datos

Librerias

library(readxl)
library(randomForest)
library(caret)
library(dplyr)

Lectura de los datos

fenotipo <- "SM300_Efficiency"
fenotipos_path <- "C:/Users/ca_re/OneDrive/Desktop/Memoria/Desarrollo/fenotipos.xlsx"
lectura_datos <- read_excel(fenotipos_path, 
                            col_types = c(rep("text", 55), rep("numeric", 12)),
                            na = "NA")

Preparación de los datos

Se obtiene el valor para fenotipo SM300 Efficiency y valor para cada cepa en los 54 SNP’s. Se eliminan valores que no poseen información. Además se alamcenan los nombres del gen de los 54 SNP’s

datos <- setNames(data.frame(lectura_datos[[fenotipo]],lectura_datos[, 2:55]), 
                  c(fenotipo, names(lectura_datos[, 2:55]))) #Se obtiene valor fenotipo y bases de los SNP's (genes)

datos <- na.omit(datos) # Se eliminan filas que no contengan datos
datos <- datos %>% mutate_at(vars(2:55), as.factor) # Se indican las columnas que son factor
nombres_genes <- names(datos[,2:55]) # Se almacenan los nombres delos genes
datos

Funciones Modelo de regresión

Converción de datos

Debido a que los valores de cada uno de los SNP’s son bases nitrogenadas se deben transformar los datos a variables dummy para poder ser utilizadas en el modelo de regresión.

# Descripción: Convierte las columnas categoricas de un dataframe a valores a dummy
# Argumentos: datos (df a convertir), nombre del fenotipo 
# Valor retorno: dataFrame (columnas del dataFrame convertidas a variables dummy)
convertir_dummy <- function(datos, nombre_fenotipo, nombre_gen){
  datos_dummy <- as.data.frame(model.matrix(~ . - 1, data = datos))
  return(datos_dummy)
}

División datos

Al ser un conjunto de datos reducidos se utiliza dividen los datos en 70% entrenamiento y un 30% de prueba.

# Descripción: Toma un conjunto de datos y los divide en datos de entrenamiento y prueba
# Argumentos: datos(conjunto de datos a separar)
# Valor retorno: lista(con el conjunto de entrenamieto y prueba)
separacion_datos <- function(datos){
  set.seed(123)
  train_index <- sample(1:nrow(datos), 0.7 * nrow(datos))
  train <- datos[train_index, ]
  test <- datos[-train_index, ]
  return(list(train = train, test = test))
}

Entrenamiento

Se realiza un entrenamiento de un modelo Random Forest, tomando los datos de entrenamiento (70% de los datos totales) con una validación cruzada simple debido a la cantidad reducida de datos, evitando el sobreajuste o subajuste de los datos, a través de la creación de varias particiones de los datos, entrenando y probando el modelo con cada una de ellas, para esta ocasión se utilizan 10 particiones.

# Descripción: Entrena un modelo Random Forest sobre un conjunto de datos de entrenamiento
# Argumentos: train(conjunto de datos de entrenamiento), fenotipo (nombre del fenotipo, variable independiente)
# Valor retorno: modelo entrenado
entrenamiento <- function(train, fenotipo){
  model <- train(as.formula(paste(fenotipo, "~ .")), 
                 data = train, 
                 method = "rf", 
                 trControl = trainControl(method = "cv", number = 10),
                 tuneLength = 2)
  return(model)
}

Predicciones

# Descripción: Se predicen valores con el modelo anteriormente entrenado
# Argumentos: model (modelo entrenado), test (datos de prueba)
# Valor retorno: predicciones
predicciones <- function(model, test){
  predictions <- predict(model, newdata = test)
  return(predictions)
}

Evalaución

  • Coeficiente de determinación (R2): medida estádistica que refelja la relación entre las variables y qué tan bien predice los valores de la variable independiente. Se aproxima el resultado con 4 dígitos. Su ecuación es: \[R^2 = \frac{\sum_{t=1}^{T}(y_t - \hat{y_t})^2}{\sum_{t=1}^{T}(y_t - \bar{y})^2}\]
  • Coeficiente de Pearson: utilizado para cuantificar la relación lineal entre dos variables cuantitativas continuas.
  • Error cuadrático medio (MSE): utilizado para evaluar la calidad de un modelo de regresión. Indicando cuanto se desvia el modelo generado de los datos observados. Su ecuación es: \[MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y_i})^2\]
# Descripción: Se calcula r2, coeficiente de pearson y mse
# Argumentos: predictions (predicciones), test (df de prueba), fenotipo (nombre del fenotipo)
# Valor retorno: lista (con valores test y train)
evaluacion_coeficientes <- function(predictions, test, fenotipo){
  # Coeficiente de determinación
  n = length(predictions)
  x_mean = mean(predictions)
  y_mean = mean(test[[fenotipo]])
  ssx = sum((predictions - x_mean)^2)
  ssy = sum((test[[fenotipo]] - y_mean)^2)
  sxy = sum((predictions - x_mean)*(test[[fenotipo]] - y_mean))
  r2 = round((sxy^2) / (ssx*ssy),5)

  # Coeficiente de Pearson
  pearson <- round(cor(test[[fenotipo]],predictions)^2,5)

  # Error cuadrático medio (MSE)
  mse <- round(mean((test[[fenotipo]] - predictions)^2),5)
  
  # Se organizan los resultados
  #resultados <- data.frame(Gen = nombres_genes[i], R2 = r2, Pearson = pearson, MSE = mse)
  resultados <- data.frame(Gen = nombres_genes[i], R2 = r2, MSE = mse)
  return(resultados)
}

Ejecución

Se realizan 54 modelos (uno para cada SNP) a tarvés de las funciones nombradas anteriormente. Finalmente, se plasman los resultados en un dataframe.

  resultados <- data.frame(matrix(ncol=5,nrow=0)) # Dataframe donde se guardaran los resultados
  for(i in 1:length(nombres_genes)){
    # Se obtiene la información del gen 
    datos_gen <- data.frame(matrix(ncol=2,nrow=344))
    datos_gen <- setNames(data.frame(datos[[fenotipo]],datos[,i+1]),
                          c(fenotipo, colnames(datos)[i+1]))
    
    # Se transforma la columna del SNP a variable dummy para ser trabajada en el modelo
    datos_transformados <- convertir_dummy(datos_gen,fenotipo, colnames(datos)[i+1])
    
    # Se separan los datos en entrenamiento y prueba (70/30)
    datos_separados <- separacion_datos(datos_transformados)
    train <- datos_separados$train
    test <- datos_separados$test
    
    # Se realiza un modelo Random forest (con validación cruzada simple)
    model <- entrenamiento(train,fenotipo)
    
    # Se utilizan los datos de prueba
    predictions <- predicciones(model,test)
    
    # Se realiza la evaluación de las predicciones del modelo
    evaluacion <- evaluacion_coeficientes(predictions, test, fenotipo)
    
    # Se almacenan los resultados obtenidos del gen
    resultados <- rbind(resultados, evaluacion)
  }
resultados

Resultados

#Mayor R²
indice_max <- which.max(resultados$R2)
cat("Mayor R²:",resultados$R2[indice_max], ", MSE:",resultados$MSE[indice_max],", Nombre gen: ",resultados$Gen[indice_max])
## Mayor R²: 0.22245 , MSE: 0.01455 , Nombre gen:  YLR425W
#Menor MSE
indice_min <- which.min(resultados$MSE)
cat("\nMenor MSE:",resultados$MSE[indice_min], ", R²:",resultados$R2[indice_min],", Nombre gen: ",resultados$Gen[indice_min])
## 
## Menor MSE: 0.01455 , R²: 0.22245 , Nombre gen:  YLR425W