GAM models for Neuroscience - Code 101

GAMS FOR NEURODEGENERATIVE DISEASES

Vale, esta es un summary de la primera parte de la tesis, o al menos de la idea original que teníamos pensada y que ahora está ligeramente aparcada. En muy resumidas cuentas la idea es es crear la database inicial a partir de los PET y MRI ya normalizados con Matlab y SPM. Una vez esta es un data.frame estable con todos nuestros datos (como también será para los SCC) podremos hacer modelos de regresión no paramétricos en que plotear funciones que unan puntos con la misma intensidad de actividad cerebral estimada.

El objetivo inicial es obtener imágenes (gam.vis y gam.check) similares a las que se exponen debajo. En el código que presento aquí se resume el proceso para obtenerlas, aunque todavía queda una etapa extra por llevar a cabo que sería la de comparar estos resultados con los nativos obtenidos con SPM, un software implementado con Matlab que es el estándar de mercado en el mundillo del diagnóstico por imagen médica (por lo menos para el caso de la neuroimagen). Mientras no pueda demostrar fehacientemente que mi metodología es mejor que la de SPM, solo serán imágenes muy bonitas.

#install.packages(c("cowplot","magick"))
library(magick);library(cowplot);library(ggplot2)

p1 <- ggdraw() + draw_image("GAMexample1.png")
p2 <- ggdraw() + draw_image("GAMexample2.png")
p3 <- ggdraw() + draw_image("GAMexample3.png")
p4 <- ggdraw() + draw_image("GAMexample4.png")

plot_grid(
  p1, p2,p3, p4,
  ncol = 2)

Claro que este es el mismo problema que son los SCC’s, la comparabilidad a nivel matemático. Claro que en SCCs voy un paso por delante ya que sí que puedo obtener imágenes que, a nivel visual, sean prácticamente idénticas en formato y por lo menos demostrar visualmente que los resultados son similares ya es un resultado positivo. Quizás no sea el resultado final deseado pero al menos estaría demostrando que hay otra vía de obtener resultados similares, más simple, y desde un approach de datos funcionales.

Pero ahora lo que nos atiene es obtener dichos modelos sobre diferentes cortes y ver su utilidad a nivel visual y qué nos indican. Cómo compararlas entre ellas y cómo compararlas con SPM es otra historia bien distinta.

PREÁMBULO:

#install.packages(c("mgcv","gamair","oro.nifti","memsic"))
library(mgcv);library(gamair);library(oro.nifti);library(memisc)

VISUALIZACIONES BÁSICAS:

Para la visualización de imágenes PET o MRI yo no recomendaría en absoluto utilizar R. Existen una serie de paquetes para trabajar con imágenes pero está claro que el objetivo de estos es trabajar con los datos subyacentes a la imagen, no tanto la imagen en si. Para trabajar con imágenes a un nivel visual recomendaría más bien el uso de SPM (clásico) o sino de MRIcro o de Mango. Cualquiera de ellas se parece más a un editor de imágenes estándar, algo así como un photoshop de cerebros. De todos modos, a veces puede ser necesario comprobar sin más que todo está correcto, extraer algún valor puntual de un pixel, hacer una visualización ortogonal…

Aunque el número de funciones disponibles es limitado, aquí van algunas de las más útiles:

LECTURA DE IMAGENES Y DATOS GENERALES:

img_003_S_1059 <- readNIfTI("003_S_1059", verbose = FALSE, warn = -1, reorient = TRUE,
                 call = NULL, read_data = TRUE)
img_003_S_1059

## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 16 (FLOAT32)
##   Bits per Pixel  : 32
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 2 (Aligned_Anat)
##   Sform Code      : 2 (Aligned_Anat)
##   Dimension       : 79 x 95 x 79
##   Pixel Dimension : 2 x 2 x 2
##   Voxel Units     : mm
##   Time Units      : sec

VISUALIZACIONES MÁS ÚTILES

Una lista de las diferentes visualizaciones que podemos obtener aunque, como bien he dicho, son mucho más completas e interactivas con otros programas.

Axial:

image(img_003_S_1059) #axial por defecto

Coronal:

image(img_003_S_1059,plane="coronal") #coronal

Sagital:

image(img_003_S_1059,plane="sagittal") #sagital

Plano ortográfico:

orthographic(img_003_S_1059,text="Plano Ortografico PET") #plano ortografico

## Warning in min(x, na.rm = na.rm): ningún argumento finito para min; retornando
## Inf

## Warning in max(x, na.rm = na.rm): ningun argumento finito para max; retornando -
## Inf

DATOS ASOCIADOS:

Estas imágenes, por supuesto, tienen una serie de datos asociados (sino ya me dirás qué clase de data science iba yo a hacer aquí). En general estos datos solo me valen de algo en bulk para sacar un data.frame útil. No obstante, en potencia puedo pedir el valor de un dato concreto en una coordenada concreta y esto puede ser útil para comprobar en ciertos escalones del desarrollo del código que todo vaya correctamente. Por si a caso, aquí queda:

Datos asociados a la imagen:

img_data=img_data(img_003_S_1059)
qplot(as.vector(img_data), 
      geom="histogram",
      main = "Histogram for PET values", 
      xlab = "PET levels",  
      fill=I("blue"), 
      col=I("red"), 
      alpha=I(.2))

negatives=img_data[img_data<0]
length(negatives) # de hecho hay muchos valores inferiores a 0 por lo que lo más rápido será, en la database, igualarlos a 0 (MRIcro los identifica como 0's d ehecho)

## [1] 7590

Aquí por ejemplo ya vemos en el histograma que hay muchos valores cercanos al cero o incluso negativos por lo que más adelante habrá que tenerlo en cuenta. Esto es probablemente derivado de un proceso de normalización que ha reducido valores nulos por debajo de 0.

Coordenadas concretas en la imagen:

#Diferentes coordenadas de la imagen con valores negativos (comprobacion)
img_data[1,2,3]

## [1] 0

img_003_S_1059[45,45,16]

## [1] 0.5555722

img_003_S_1059[10,70,40];img_data[10,70,40] # valdría cualquiera de los dos comandos

## [1] 0.3873922

## [1] 0.3873922

img_data[66,34,24]

## [1] 0.8744069