Vamos a mirar ADNI, a ver que sale. En principio habria que ajustar tanot por edad como por volumen craneal. A primera vista pareceria una buena opcion, pero una vez que ajusto por edad y volumen total,

library(ADNIMERGE)
pop <- adnimerge[adnimerge$VISCODE=="bl" & (adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
m1 <- lm(pop$Ventricles ~ pop$ICV)
pop$adj_Ventricles = pop$Ventricles - m1$coefficients[[2]]*(pop$ICV- mean(pop$ICV,na.rm=TRUE))
m2 <- lm(pop$adj_Ventricles ~ pop$AGE)
pop$vc_Ventricles = pop$adj_Ventricles - m2$coefficients[[2]]*(pop$AGE-mean(pop$AGE,na.rm=TRUE))

Y queda entonces,

uncorrected	corrected

Los valores estan demasiado juntos y la correccion los agrupa mas.

al escoger el hipocampo estoy presuponiendo que hay AD, pero vamos a ver que queda,

uncorrected	corrected

Pues no esta mal.

   UCSFVOLs <- c("FLDSTRENG", "FSVERSION", "IMAGEUID", "Ventricles", 
        "Hippocampus", "WholeBrain", "Entorhinal", "Fusiform", 
        "MidTemp", "ICV")

        ucsfvol$ICV <- ucsfvol$ST10CV
        ucsfvol$Hippocampus <- (ucsfvol$ST29SV + ucsfvol$ST88SV)
        ucsfvol$Entorhinal <- (ucsfvol$ST24CV + ucsfvol$ST83CV)
        ucsfvol$Fusiform <- (ucsfvol$ST26CV + ucsfvol$ST85CV)
        ucsfvol$MidTemp <- (ucsfvol$ST40CV + ucsfvol$ST99CV)
        ucsfvol$Ventricles <- (ucsfvol$ST30SV + ucsfvol$ST37SV + 
            ucsfvol$ST89SV + ucsfvol$ST96SV)
        ucsfvol$WholeBrain <- apply(ucsfvol[, c("ST128SV", "ST17SV", 
            "ST18SV", "ST61SV", "ST16SV", "ST53SV", "ST42SV", 
            "ST29SV", "ST12SV", "ST11SV", "ST65SV", "ST76SV", 
            "ST77SV", "ST120SV", "ST75SV", "ST112SV", "ST101SV", 
            "ST88SV", "ST71SV", "ST70SV", "ST124SV", "ST147SV", 
            "ST148SV", "ST150SV", "ST151SV")], 1, sum)

Usando el diccionario en ADNI podemos hacer algo como,

Fusiform =  rh.fusiform.GrayVol + lh.fusiform.GrayVol
MidTemp = rh.middletemporal.GrayVol + lh.middletemporal.GrayVol

y entonces,

pop$ND = ifelse(pop$DX == "Dementia", 1, 0)
xt <- pop[, c("Hippocampus", "Entorhinal", "Ventricles", "Fusiform", "MidTemp","AGE","ICV","DX")]

Esta es la cosa. Tengo un monton de informacion que podria clasificarme a los sujetos entre Neurodegeneracion SI y Neurodegeneracion NO. Y si tres variables juntas me pueden dar la info mas o menos adecuadad lo suyo es que meta todo lo relevante e intente usarlo al mismo tiempo. Lo primero deberia ser construir el vector de training a partir de los sujetos de ADNI.

> library("e1071")
> library("caret")
> library("caTools")
> library("ADNIMERGE")
> pop <- adnimerge[adnimerge$VISCODE=="bl" & (adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
> pop$ND = ifelse(pop$DX == "Dementia", 1, 0)
> xt <- pop[, c("Hippocampus", "Entorhinal", "Ventricles","AGE","ICV","ND")]

Este es el vector de training pero voy a probar el algoritmo tambien asi que lo que hare es dividirlo y ver como funciona. Voy a tomar el 70% para entrenar el algoritmo y el restante 30_% para ver los resultados.

> split <- sample.split(xt, SplitRatio = 0.7)
> train_cl <- subset(xt, split == "TRUE")
> test_cl <- subset(xt, split == "FALSE")
> set.seed(120)
> classifier_cl <- naiveBayes(ND ~ ., data = train_cl)
> y_pred <- predict(classifier_cl, newdata = test_cl)
> cm <- table(test_cl$ND, y_pred)
> confusionMatrix(cm)
Confusion Matrix and Statistics

   y_pred
      0   1
  0 248  29
  1  46  84
                                          
               Accuracy : 0.8157          
                 95% CI : (0.7746, 0.8522)
    No Information Rate : 0.7224          
    P-Value [Acc > NIR] : 7.984e-06       
                                          
                  Kappa : 0.5609          
                                          
 Mcnemar's Test P-Value : 0.06467         
                                          
            Sensitivity : 0.8435          
            Specificity : 0.7434          
         Pos Pred Value : 0.8953          
         Neg Pred Value : 0.6462          
             Prevalence : 0.7224          
         Detection Rate : 0.6093          
   Detection Prevalence : 0.6806          
      Balanced Accuracy : 0.7935          
                                          
       'Positive' Class : 0 

Y no es perfecto pero esta bien!

Asi que para evaluar los sujetos habria que construir una matriz de datos identica a la de ADNI. Esto no es tan dificil, asi que vamos,

[osotolongo@brick03 bioface]$ xnat_pullfs.pl -s aseg -x bioface19 -o bf_base_aseg.csv
[osotolongo@brick03 bioface]$ xnat_pullfs.pl -s aparc -x bioface19 -o bf_base_aparc.csv
[osotolongo@brick03 bioface]$ join -t, bf_base_aseg.csv bf_base_aparc.csv > bf_base.csv
[osotolongo@brick03 bioface]$ awk -F"," '{print $3","$5}' bioface_np.csv | sed 's/Subjecte/Subject_ID/; s/Edad/AGE/' > bf_age.csv
[osotolongo@brick03 bioface]$ join -t, bf_base.csv bf_age.csv > bf_data.csv

Lo demas voy a hacerlo en R,

> base <- read.csv("bf_data.csv")
> base$Hippocampus = base$Left.Hippocampus + base$Right.Hippocampus
> base$Entorhinal = base$lh.entorhinal.GrayVol + base$rh.entorhinal.GrayVol
> base$Ventricles <- base$Left.Inf.Lat.Vent + base$Right.Inf.Lat.Vent + base$Left.Lateral.Ventricle + base$Right.Lateral.Ventricle
> base$ICV = base$eTIV

Ahora el procedimiento,

> classifier_cl <- naiveBayes(ND ~ ., data = xt)
> base$ND <- predict(classifier_cl, newdata = base)

A ver como queda,

> plot(base$AGE, base$Hippocampus, main = "Hippocampus Volume versus Age", xlab="Age", ylab="HV", pch=19, col=ifelse(base$ND==1,"red","blue"))

Nice!

A partir de la base de datos de ADNI, tomo los sujetos diagnosticados como Dementia y CN en el baseline, así como su edad y ciertos valores calculados a partir de la segmentacion. Tomo el diagnostico de Dementia como positivo en neurodegeneracion y negativo en caso contrario.

library("e1071")
library("caret")
library("caTools")
library("ADNIMERGE")
pop <- adnimerge[adnimerge$VISCODE=="bl" & (adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
pop$ND = as.factor(ifelse(pop$DX == "Dementia", 1, 0))
xt <- pop[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "ND")]

Con estos datos, entreno el clasificador,

classifier_cl <- naiveBayes(ND ~ ., data = xt)

Ahora, tomo un proyecto y reuno los datos de segmentacion y demograficos. No es neceario todo pero la edad y la segmentacion basica si. Esto puede cambiar de proyeto en proyecto pero es mas o menos asi, (lo repito por tenerlo todo junto)

xnat_pullfs.pl -s aseg -x bioface19 -o bf_base_aseg.csv
xnat_pullfs.pl -s aparc -x bioface19 -o bf_base_aparc.csv
join -t, bf_base_aseg.csv bf_base_aparc.csv > bf_base.csv
awk -F"," '{print $3","$5}' bioface_np.csv | sed 's/Subjecte/Subject_ID/; s/Edad/AGE/' > bf_age.csv
join -t, bf_base.csv bf_age.csv > bf_data.csv

Cargo estos datos y construyo las variables identicas a la de los datos de ADNIMERGE.

setwd("/home/data/bioface")
base <- read.csv("bf_data.csv")
base$Hippocampus = base$Left.Hippocampus + base$Right.Hippocampus
base$Entorhinal = base$lh.entorhinal.GrayVol + base$rh.entorhinal.GrayVol
base$Ventricles <- base$Left.Inf.Lat.Vent + base$Right.Inf.Lat.Vent + base$Left.Lateral.Ventricle + base$Right.Lateral.Ventricle
base$MidTemp = base$lh.middletemporal.GrayVol + base$rh.middletemporal.GrayVol
base$ICV = base$eTIV

Ahora puedo calcular el valor de (N),

base$ND <- predict(classifier_cl, newdata = base)

o la probabilidad de que exista neurodegeneración

base$post <- predict(classifier_cl, newdata = base, type="raw")

Clasificación	Probabilidades de N+

Y esto voy a hacerlo para FACEHBI también para compararlo,

setwd("/home/data/facehbi")
read.csv("face_data.csv") -> face
face$Hippocampus = face$Left.Hippocampus + face$Right.Hippocampus
face$Entorhinal = face$lh.entorhinal.GrayVol + face$rh.entorhinal.GrayVol
face$Ventricles <- face$Left.Inf.Lat.Vent + face$Right.Inf.Lat.Vent + face$Left.Lateral.Ventricle + face$Right.Lateral.Ventricle
face$MidTemp = face$lh.middletemporal.GrayVol + face$rh.middletemporal.GrayVol
face$ICV = face$eTIV
face$ND <- predict(classifier_cl, newdata = face)
face$post <- predict(classifier_cl, newdata = face, type="raw")

Clasificación	Probabilidades de N+

Vamos a intentar con métodos no lineales. Ojo, estos métodos no tienen interpretación probabilistica. Son construcciones matemáticas que funcionan.

BIOFACE

rf_cl <- randomForest(ND ~ ., data = xt, na.action = na.omit)
base$ND <- predict(rf_cl, newdata = base)
base$post <- predict(rf_cl, newdata = base, type="prob")
plot(base$AGE, base$Hippocampus, main = "Hippocampus Volume versus Age", xlab="Age", ylab="HV", pch=19, col=ifelse(base$ND==1,"red","green"))
hist(base$post[,2], col="red", density = 50, probability=TRUE, breaks = 10, main="Probabilities profile", xlab = "Probability of N+")

Clasificación	Probabilidades de N+

FACEHBI

face$ND <- predict(rf_cl, newdata = face)
face$post <- predict(rf_cl, newdata = face, type = "prob")
plot(face$AGE, face$Hippocampus, main = "Hippocampus Volume versus Age", xlab="Age", ylab="HV", pch=19, col=ifelse(face$ND==1,"red","green"))
hist(face$post[,2], col="red", density = 50, probability=TRUE, breaks = 10, main="Probabilities profile", xlab = "Probability of N+")

Clasificación	Probabilidades de N+

BIOFACE

svm_cl <- svm(ND ~ ., data = xt, na.action = na.omit)
base$ND <- predict(svm_cl, newdata = base)
plot(base$AGE, base$Hippocampus, main = "Hippocampus Volume versus Age", xlab="Age", ylab="HV", pch=19, col=ifelse(base$ND==1,"red","green"))

FACEHBI

face$ND <- predict(svm_cl, newdata = face)
plot(face$AGE, face$Hippocampus, main = "Hippocampus Volume versus Age", xlab="Age", ylab="HV", pch=19, col=ifelse(face$ND==1,"red","green"))

¿Son las afecciones vasculares parte de la neurodegeneracion? ¿Debe ser considerado el aumento en el volumen de WMH una parte fundamental de la evaluacion de N?

Pues si queremos hacerlo de esta manera, hay basicamente dos maneras. Utilizando las WMH de Freesurfer (T1w) o las calculadas por un soft independiente (T2w). Si tenemos un FLAIR 3D podemos hacerlo de ambas maneras pero en caso contrario.

Las WMH vienen ya calculadas en el procedimiento standard asi que esto deberia ser painless. Pero en la DB de ADNIMERGE, no vienen integradas en la tabla de resumen que estamos usando, asi que hay que incluirlas. Pero revisando el procedimiento de integracion (adnimerger) veo que hay que arreglar un poco las VISCODE.

ucsffsx -> ucsffsx_t
ucsffsx_t$VISCODE <- ifelse(ucsffsx_t$VISCODE == "sc", "bl", ucsffsx_t$VISCODE)
a1 <- merge(pop, ucsffsx_t, by=c("RID", "VISCODE"))
a1t <- a1[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "ST128SV", "ND")]
...
...
axt <- rbind(a1t, a2t, a3t)
xt <- rename(axt, "WMH" = "ST128SV")

y esto hay que hacerlo para cada una de las tablas y despues concatenar los resutados. Ahora habria tambien que introducir las WMH en la DB del test,

base$WMH = base$WM.hypointensities

que se dice un poco mas rapido de lo que se hace, tambien es cierto,

library("e1071")
library("caret")
library("caTools")
library("tidyverse")
library("ADNIMERGE")
 
input_file="input_data.csv"
output_file="classifier_output_w.csv"
output_fig="classifier_output_hippocampus_w.ps"
pop <- adnimerge[(adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
pop$ND = as.factor(ifelse(pop$DX == "Dementia", 1, 0))
ucsffsx -> ucsffsx_t
ucsffsx_t$VISCODE <- ifelse(ucsffsx_t$VISCODE == "sc", "bl", ucsffsx_t$VISCODE)
a1 <- merge(pop, ucsffsx_t, by=c("RID", "VISCODE"))
ucsffsx51 -> ucsffsx_t
ucsffsx_t$VISCODE <- ifelse(ucsffsx_t$VISCODE == "sc", "bl", ucsffsx_t$VISCODE)
a2 <- merge(pop, ucsffsx_t, by=c("RID", "VISCODE"))
ucsffsx6 -> ucsffsx_t
ucsffsx_t$VISCODE <- ifelse(ucsffsx_t$VISCODE == "sc", "bl", ucsffsx_t$VISCODE)
a3 <- merge(pop, ucsffsx_t, by=c("RID", "VISCODE"))
a1t <- a1[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "ST128SV", "ND")]
a2t <- a2[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "ST128SV", "ND")]
a3t <- a3[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "ST128SV", "ND")]
axt <- rbind(a1t, a2t, a3t)
xt <- rename(axt, "WMH" = "ST128SV")
classifier_cl <- naiveBayes(ND ~ ., data = xt)
base <- read.csv(input_file)
base$Hippocampus = base$Left.Hippocampus + base$Right.Hippocampus
base$Entorhinal = base$lh.entorhinal.GrayVol + base$rh.entorhinal.GrayVol
base$Ventricles <- base$Left.Inf.Lat.Vent + base$Right.Inf.Lat.Vent + base$Left.Lateral.Ventricle + base$Right.Lateral.Ventricle
base$MidTemp = base$lh.middletemporal.GrayVol + base$rh.middletemporal.GrayVol
base$ICV = base$eTIV
base$WMH = base$WM.hypointensities
base$ND <- predict(classifier_cl, newdata = base)
base2e = base[, c("Subject_ID", "ND")]
write.csv(base2e, file=output_file, row.names=FALSE, quote=FALSE)
a <- lm(base$Hippocampus ~ base$ICV)
base$aHV = base$Hippocampus - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aHV, main = "Hippocampus volume versus Age", xlab="Age", ylab="adjusted HV", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()
a <- lm(base$MidTemp ~ base$ICV)
base$aMidTemp = base$MidTemp - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
output_fig="classifier_output_middletemporal_w.ps"
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aMidTemp, main = "Middle temporal cortex volume versus Age", xlab="Age", ylab="adjusted MidTemp Volume", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()
a <- lm(base$Entorhinal ~ base$ICV)
base$aEntorhinal = base$Entorhinal - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
output_fig="classifier_output_entorhinal_w.ps"
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aEntorhinal, main = "Entorhinal cortex volume versus Age", xlab="Age", ylab="adjusted Entorhinal Volume", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()

OK, ya puestos también podemos hacerlo con las WMH calculadas desde T2W. Yo uso PGS, para las WMH de T2w, asi que tengo que pegarlas a las tablas,

(Ejemplo con facehbi)

[osotolongo@brick03 facehbi]$ sed 's/;/,/' facehbi_mri.csv > codes.csv
[osotolongo@brick03 facehbi]$ sed -i '1iSubject,Subject_ID' codes.csv
[osotolongo@brick03 facehbi]$ join -t, codes.csv facehbi_wmh_metrics.csv | awk -F',' '{print $2","$3}' > facehbi_wmh.csv
[osotolongo@brick03 facehbi]$ ./xnat_get_age.pl -x facehbi
[osotolongo@brick03 facehbi]$ xnat_pullfs.pl -s aseg -p facehbi -o base_aseg.csv
[osotolongo@brick03 facehbi]$ xnat_pullfs.pl -s aparc -p facehbi -o base_aparc.csv
[osotolongo@brick03 facehbi]$ join -t, base_aseg.csv base_aparc.csv > base_full.csv
[osotolongo@brick03 facehbi]$ (head -n 1 facehbi_age_data.csv && tail -n +2 facehbi_age_data.csv | sort -t, -n)> facehbi_age_data_sorted.csv
[osotolongo@brick03 facehbi]$ join -t, facehbi_age_data_sorted.csv facehbi_wmh.csv > facehbi_age_wmh.csv
[osotolongo@brick03 facehbi]$ join -t, base_full.csv facehbi_age_wmh.csv > input_data.csv

Claro que aqui hay que usar la tabla de UCD WMH,

pop <- adnimerge[(adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
ucd_wmh -> ucd_t
ucd_t$VISCODE <- ifelse(ucd_t$VISCODE == "scmri", "bl", ucd_t$VISCODE)
ucd_t$VISCODE <- ifelse(ucd_t$VISCODE == "sc", "bl", ucd_t$VISCODE)
pop <- merge(pop, ucd_t, by=c("RID", "VISCODE"))

y despues modificar un poco las de nuestro proyecto para que cuadren las unidades,

base$TOTAL_WMH = base$WMH*0.001

OK, ya puestos también podemos hacerlo con las WMH calculadas desde T2W

library("e1071")
library("caret")
library("caTools")
library("ADNIMERGE")
input_file="input_dataw.csv"
output_file="classifier_output_w2.csv"
output_fig="classifier_output_hippocampus_w2.ps"
pop <- adnimerge[(adnimerge$DX=="Dementia" | adnimerge$DX=="CN"),];
ucd_wmh -> ucd_t
ucd_t$VISCODE <- ifelse(ucd_t$VISCODE == "scmri", "bl", ucd_t$VISCODE)
ucd_t$VISCODE <- ifelse(ucd_t$VISCODE == "sc", "bl", ucd_t$VISCODE)
pop <- merge(pop, ucd_t, by=c("RID", "VISCODE"))
pop$ND = as.factor(ifelse(pop$DX == "Dementia", 1, 0))
xt <- pop[, c("Hippocampus", "Entorhinal", "Ventricles", "MidTemp", "AGE", "ICV", "TOTAL_WMH", "ND")]
classifier_cl <- naiveBayes(ND ~ ., data = xt)
base <- read.csv(input_file)
base$Hippocampus = base$Left.Hippocampus + base$Right.Hippocampus
base$Entorhinal = base$lh.entorhinal.GrayVol + base$rh.entorhinal.GrayVol
base$Ventricles <- base$Left.Inf.Lat.Vent + base$Right.Inf.Lat.Vent + base$Left.Lateral.Ventricle + base$Right.Lateral.Ventricle
base$MidTemp = base$lh.middletemporal.GrayVol + base$rh.middletemporal.GrayVol
base$ICV = base$eTIV
base$TOTAL_WMH = base$WMH*0.001
base$ND <- predict(classifier_cl, newdata = base)
base$posterior <- predict(classifier_cl, newdata = base, type="raw")
base2e = base[, c("Subject_ID", "ND", "posterior")]
write.csv(base2e, file=output_file, row.names=FALSE, quote=FALSE)
a <- lm(base$Hippocampus ~ base$ICV)
base$aHV = base$Hippocampus - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aHV, main = "Hippocampus volume versus Age", xlab="Age", ylab="adjusted HV", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()
a <- lm(base$MidTemp ~ base$ICV)
base$aMidTemp = base$MidTemp - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
output_fig="classifier_output_middletemporal_w2.ps"
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aMidTemp, main = "Middle temporal cortex volume versus Age", xlab="Age", ylab="adjusted MidTemp", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()
a <- lm(base$Entorhinal ~ base$ICV)
base$aEntorhinal = base$Entorhinal - a$coefficients[[2]]*(base$ICV - mean(base$ICV, na.rm=TRUE))
output_fig="classifier_output_entorhinal_w2.ps"
postscript(output_fig, width=1024, height=600, bg="white")
plot(base$AGE, base$aEntorhinal, main = "Entorhinal cortex volume versus Age", xlab="Age", ylab="adjusted Entorhinal", pch=19, col=ifelse(base$ND==1,"red","green"))
dev.off()

(NO WAY, explicar como se hace que hay que meter el ucd_wmh de ADNI y sacar wl WMH con PGS)

Con BIOFACE:

no WMH	T1w WMH	T2w WMH

Con FACEHBI:

no WMH	T1w WMH	T2w WMH