Integrando ENIGMA Freesurfer QC

Para casi todos los proyectos de colaboracion, los procedimientos de ENIGMA se han tomado como el standard a seguir. Asi que vamos a intentar integrar estos procedimientos en el pipeline, sobre todo en el proyecto MRIFACE, que es el mas extenso y dificil de manipular. A partir de aqui deberia ser simple aplicar el procedimeinto al resto de proyectos.

Nota: El procedimiento es bastante básico y tiene dos evaluaciones simples: Good, Fail. Los minor errors pasaran como buenos, los major errors como malos, en general. No obstante, el objetivo es guardar todas las imagenes en XNAT como recurso del experimento, asi que si es necesario basta bajar las imagenes y volver a revisar el QC. Esto es muchisimo mas sencillo y rapido que volver a ejecutar todo el procedimiento para un puñado de MRIs cuyos resultados se quieran compartir a algun proyecto.

Arreglando scripts

Primero, hay que obtener los scripts y documentación de trabajo de los protocolos de enigma.

Bajando resultados de Freesurfer

El paso previo a ejecutar el QC es obtener los resultados del procesamiento y ordenarlos en un directorio. Para poder evaluar paulatinamente los sujetos, hay que permitir dos cosas,

elegir el directorio donde vamos a guardar los resultados
suministrar una lista de los sujetos a evaluar

Esto se hace tomando estos valores desde la linea de comandos (opcional),

while (@ARGV and $ARGV[0] =~ /^-/) {
        $_ = shift;
        last if /^--$/;
#...
        if (/^-i/) { $ifile = shift; chomp($ifile);}
        if (/^-o/) { $odir = shift; chomp($odir);}       
#...
}

Luego, de la posible lista de sujetos (o del proyecto completo), debo obtener los experimentos correspondientes. Cada experimento procesado debe tener el directorio completo de FS compactado y asociado como recurso. Este lo bajamos y lo descomprimimos en su sitio.

Todo esto lo hacemos reutilizando las bibiotecas del pipeline

Facilmente

get_fs_data.pl

#!/usr/bin/perl
 
# Copyright 2023 O. Sotolongo <asqwerty@gmail.com>
 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
 
# Main goal here is to get all the FS out put data from XNAT.
# Basically, for each target subject, get experiments, download the 
# gziped tar file with all the data and explode in a separated directory.
#
# And from there you can make your favorite FS QC procedure
#
 
use strict; use warnings;
use Cwd qw(cwd);
use NEURO4 qw(load_project print_help populate check_or_make getLoggingTime);
use XNATACE qw(xget_session xget_mri xlist_res xget_res_file xget_subjects);
use File::Temp qw(:mktemp tempdir);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
my $prj;
my $xprj;
my $ifile; my $efile;
my $wdir = cwd;
my $odir;
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
	$_ = shift;
	last if /^--$/;
	if (/^-i/) { $ifile = shift; chomp($ifile);}
	if (/^-o/) { $odir = shift; chomp($odir);}
	if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/reviewqc.hlp'; exit;}
	if (/^-x/) { $xprj = shift; chomp($xprj);} #nombre del proyecto en XNAT
	if (/^-p/) { $prj = shift; chomp($prj);} #nombre local del proyecto
}
 
# Mira, hay que meter el proyecto de XNAT con alguno de los dos switch
if ($prj and not $xprj) {
	my %pdata = load_project($prj);
	$xprj = $pdata{'XNAME'};
}
# O te vas a tomar por culo
die "Should supply XNAT project name or define it at local project config!\n" unless $xprj;
 
# y entonces, el file de input es el que tu digas o los tendras que mirar todos
# Estos son los experimentos a los que vamos a revisar o hacer el QC. 
# Puedes tomar la lista completa y reducirla a lo que quieres o simplemente 
# tomar la lista as is
my $tmp_dir = $ENV{TMPDIR};
my %xconf = xget_session();
my @pollos;
my %epollos;
$efile = $wdir.'/'.$xprj.'_experiment.list';
my $selection = '*/mri/{orig,orig_nu,aparc+aseg}.mgz */label/*.aparc.annot */surf/*.pial* */surf/*.orig */surf/*.white */mri/transforms/talairach.xfm */stats/*.aparc.stats */stats/aseg.stats';
#Añadir logging!!!!!
my $odebug = $wdir.'/.debug_'.$xprj.'_data_'.getLoggingTime().'.out';
open STDOUT, ">$odebug" or die "Can't redirect output";
open STDERR, ">&STDOUT" or die "Can't redirect errors";
 
if ($ifile and -f $ifile){
	$efile=$wdir.'/'.$xprj.'_experiment_custom.list';
	open IDF, "<$ifile";
	@pollos = <IDF>; chomp @pollos;
	close IDF;
	open TDF, ">$efile";
	foreach my $pollo (@pollos){
		$epollos{$pollo} = [ xget_mri($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $pollo) ];
		foreach my $pexp (@{$epollos{$pollo}}){
			print TDF "$pexp\n";
		}
	}
	close TDF;
}
unless ($odir and -d $odir) {
	# Si me das un directorio de input y este existe entonces he de avisar que existe y salir
	# Dejar que el luser arregle a su manera el problema 
	$odir=$wdir.'/'.$xprj.'_fsresults';
	if ( -d $odir ){
		# Si el default existe entonces aviso y salgo
		die "Output directory $odir already exists, move it or give me a new directory output";
	}
	mkdir $odir;
	#y ejecuto la preparacion
	if ($ifile and -f $ifile){
		# si me has dado una lista de los que quieres procesar
		# me limito a preparar esos
		my $tsdir = tempdir(TEMPLATE => $tmp_dir.'/tars_data.XXXXX', CLEANUP => 1);
		foreach my $pollo (@pollos){
			foreach my $pexp (@{$epollos{$pollo}}) {
				my $otfile = $tsdir.'/'.$pexp.'.tar.gz';
				my %fs_files = xlist_res($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $pexp, 'FS');
				foreach my $fsfile (sort keys %fs_files){
					if ($fsfile =~ /.*\.tar\.gz$/){
						xget_res_file($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $pexp, 'FS', $fsfile, $otfile);
					}
				}
				my $otdir = $odir.'/'.$pexp;
				mkdir $otdir;
				my $pre = 'tar xzf '.$otfile.' --strip-components=1 -C '.$otdir.' '.$selection;
				print "$pre\n";
				system($pre);
				unlink $otfile;
			}
		}
		unlink $tsdir;
	}else{	
		#pero si no hay lista los bajo todos
		my $tsdir = tempdir(TEMPLATE => $tmp_dir.'/tars_data.XXXXX', CLEANUP => 1);
		my %allchicks = xget_subjects($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj);
		foreach  my $pollo (sort keys %allchicks){
			$epollos{$pollo} = [ xget_mri($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $pollo) ];
			foreach my $pexp (@{$epollos{$pollo}}){
				my $otfile = $tsdir.'/'.$pexp.'.tar.gz';
				my %fs_files = xlist_res($xconf{'HOST'}, $xconf{'JSESSION'}, $pexp, 'FS');
				foreach my $fsfile (sort keys %fs_files){
                                	if ($fsfile =~ /.*\.tar\.gz$/){
                                        	xget_res_file($xconf{'HOST'}, $xconf{'JSESSION'}, $pexp, 'FS', $fsfile, $otfile);
                                	}
                        	}
				if (-e $otfile) { 
					my $otdir = $odir.'/'.$pexp;
					mkdir $otdir;
					my $pre = 'tar xzf '.$otfile.' --strip-components=1 -C '.$otdir.' '.$selection;
					print "$pre\n";
					system($pre);
					unlink $otfile;
				} else {
					print "No FS data for $pexp ($pollo)!\n";
				}
			}
		}
		unlink $tsdir;
	}
}else{
	die "This output directory already exists, choose another one!";
}

Las opciones del script get_fs_data.pl son,

-p prj : nombre local del proyecto (OBLIGATORIO)
-x xprj : nombre en XNAT del proyecto (Opcional si el proyecto local esta configurado correctamente)
-o output_dir : directorio de output (Opcional)
-i subjects.list : archivo con la lista de sujetos a bajar (Opcional)

tl;dr

[osotolongo@brick03 mri_face]$ ./get_fs_data.pl -p mriface

extract data

Basicamente es usar el script de ENIGMA para sacar SurfAVG y ThicAVG de los resultados individuales y ponerlos en una tabla, pero voy poner esto un poco mas amigable y voy a añadir un par de lineas al inicio del script para poder ejecutarlo desde cualquier sitio y me guarde los CSV en el directorio de trabajo.

#!/bin/bash
idir=$(realpath $1);
rdir=$(pwd);
 
[[ -z ${idir} ]] && { echo "Should supply input directory"; exit 1; }
[[ -d ${idir} ]] || { echo "Directory not found"; exit 404; }
cd ${idir};

El script queda mas o menos igual con alguna pequeña modificacion,

extract.sh

#!/bin/bash
idir=$(realpath $1);
rdir=$(pwd);
 
[[ -z ${idir} ]] && { echo "Should supply input directory"; exit 1; }
[[ -d ${idir} ]] || { echo "Directory not found"; exit 404; }
cd ${idir};
 
echo 'SubjID,L_bankssts_thickavg,L_caudalanteriorcingulate_thickavg,L_caudalmiddlefrontal_thickavg,L_cuneus_thickavg,L_entorhinal_thickavg,L_fusiform_thickavg,L_inferiorparietal_thickavg,L_inferiortemporal_thickavg,L_isthmuscingulate_thickavg,L_lateraloccipital_thickavg,L_lateralorbitofrontal_thickavg,L_lingual_thickavg,L_medialorbitofrontal_thickavg,L_middletemporal_thickavg,L_parahippocampal_thickavg,L_paracentral_thickavg,L_parsopercularis_thickavg,L_parsorbitalis_thickavg,L_parstriangularis_thickavg,L_pericalcarine_thickavg,L_postcentral_thickavg,L_posteriorcingulate_thickavg,L_precentral_thickavg,L_precuneus_thickavg,L_rostralanteriorcingulate_thickavg,L_rostralmiddlefrontal_thickavg,L_superiorfrontal_thickavg,L_superiorparietal_thickavg,L_superiortemporal_thickavg,L_supramarginal_thickavg,L_frontalpole_thickavg,L_temporalpole_thickavg,L_transversetemporal_thickavg,L_insula_thickavg,R_bankssts_thickavg,R_caudalanteriorcingulate_thickavg,R_caudalmiddlefrontal_thickavg,R_cuneus_thickavg,R_entorhinal_thickavg,R_fusiform_thickavg,R_inferiorparietal_thickavg,R_inferiortemporal_thickavg,R_isthmuscingulate_thickavg,R_lateraloccipital_thickavg,R_lateralorbitofrontal_thickavg,R_lingual_thickavg,R_medialorbitofrontal_thickavg,R_middletemporal_thickavg,R_parahippocampal_thickavg,R_paracentral_thickavg,R_parsopercularis_thickavg,R_parsorbitalis_thickavg,R_parstriangularis_thickavg,R_pericalcarine_thickavg,R_postcentral_thickavg,R_posteriorcingulate_thickavg,R_precentral_thickavg,R_precuneus_thickavg,R_rostralanteriorcingulate_thickavg,R_rostralmiddlefrontal_thickavg,R_superiorfrontal_thickavg,R_superiorparietal_thickavg,R_superiortemporal_thickavg,R_supramarginal_thickavg,R_frontalpole_thickavg,R_temporalpole_thickavg,R_transversetemporal_thickavg,R_insula_thickavg,LThickness,RThickness,LSurfArea,RSurfArea,ICV' > CorticalMeasuresENIGMA_ThickAvg.csv
echo 'SubjID,L_bankssts_surfavg,L_caudalanteriorcingulate_surfavg,L_caudalmiddlefrontal_surfavg,L_cuneus_surfavg,L_entorhinal_surfavg,L_fusiform_surfavg,L_inferiorparietal_surfavg,L_inferiortemporal_surfavg,L_isthmuscingulate_surfavg,L_lateraloccipital_surfavg,L_lateralorbitofrontal_surfavg,L_lingual_surfavg,L_medialorbitofrontal_surfavg,L_middletemporal_surfavg,L_parahippocampal_surfavg,L_paracentral_surfavg,L_parsopercularis_surfavg,L_parsorbitalis_surfavg,L_parstriangularis_surfavg,L_pericalcarine_surfavg,L_postcentral_surfavg,L_posteriorcingulate_surfavg,L_precentral_surfavg,L_precuneus_surfavg,L_rostralanteriorcingulate_surfavg,L_rostralmiddlefrontal_surfavg,L_superiorfrontal_surfavg,L_superiorparietal_surfavg,L_superiortemporal_surfavg,L_supramarginal_surfavg,L_frontalpole_surfavg,L_temporalpole_surfavg,L_transversetemporal_surfavg,L_insula_surfavg,R_bankssts_surfavg,R_caudalanteriorcingulate_surfavg,R_caudalmiddlefrontal_surfavg,R_cuneus_surfavg,R_entorhinal_surfavg,R_fusiform_surfavg,R_inferiorparietal_surfavg,R_inferiortemporal_surfavg,R_isthmuscingulate_surfavg,R_lateraloccipital_surfavg,R_lateralorbitofrontal_surfavg,R_lingual_surfavg,R_medialorbitofrontal_surfavg,R_middletemporal_surfavg,R_parahippocampal_surfavg,R_paracentral_surfavg,R_parsopercularis_surfavg,R_parsorbitalis_surfavg,R_parstriangularis_surfavg,R_pericalcarine_surfavg,R_postcentral_surfavg,R_posteriorcingulate_surfavg,R_precentral_surfavg,R_precuneus_surfavg,R_rostralanteriorcingulate_surfavg,R_rostralmiddlefrontal_surfavg,R_superiorfrontal_surfavg,R_superiorparietal_surfavg,R_superiortemporal_surfavg,R_supramarginal_surfavg,R_frontalpole_surfavg,R_temporalpole_surfavg,R_transversetemporal_surfavg,R_insula_surfavg,LThickness,RThickness,LSurfArea,RSurfArea,ICV' > CorticalMeasuresENIGMA_SurfAvg.csv
 
for subj_id in $(ls ${idir}); do #may need to change this so that it selects subjects with FS output
	if [[ -d ${subj_id} ]]; then 
		printf "%s,"  "${subj_id}" >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%s,"  "${subj_id}" >> CorticalMeasuresENIGMA_SurfAvg.csv
 
		for side in lh.aparc.stats rh.aparc.stats; do
 
			for x in bankssts caudalanteriorcingulate caudalmiddlefrontal cuneus entorhinal fusiform inferiorparietal inferiortemporal isthmuscingulate lateraloccipital lateralorbitofrontal lingual medialorbitofrontal middletemporal parahippocampal paracentral parsopercularis parsorbitalis parstriangularis pericalcarine postcentral posteriorcingulate precentral precuneus rostralanteriorcingulate rostralmiddlefrontal superiorfrontal superiorparietal superiortemporal supramarginal frontalpole temporalpole transversetemporal insula; do
 
				printf "%g," `grep -w ${x} ${subj_id}/stats/${side} | awk '{print $5}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
				printf "%g," `grep -w ${x} ${subj_id}/stats/${side} | awk '{print $3}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
 
			done
		done
 
		printf "%g," `cat ${subj_id}/stats/lh.aparc.stats | grep MeanThickness | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%g," `cat ${subj_id}/stats/rh.aparc.stats | grep MeanThickness | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%g," `cat ${subj_id}/stats/lh.aparc.stats | grep MeanThickness | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
		printf "%g," `cat ${subj_id}/stats/rh.aparc.stats | grep MeanThickness | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
		printf "%g," `cat ${subj_id}/stats/lh.aparc.stats | grep WhiteSurfArea | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%g," `cat ${subj_id}/stats/rh.aparc.stats | grep WhiteSurfArea | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%g," `cat ${subj_id}/stats/lh.aparc.stats | grep WhiteSurfArea | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
		printf "%g," `cat ${subj_id}/stats/rh.aparc.stats | grep WhiteSurfArea | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
		printf "%g" `cat ${subj_id}/stats/aseg.stats | grep IntraCranialVol | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_ThickAvg.csv
		printf "%g" `cat ${subj_id}/stats/aseg.stats | grep IntraCranialVol | awk -F, '{print $4}'` >> CorticalMeasuresENIGMA_SurfAvg.csv
 
		echo "" >> CorticalMeasuresENIGMA_ThickAvg.csv
		echo "" >> CorticalMeasuresENIGMA_SurfAvg.csv
	fi
done
 
mv CorticalMeasuresENIGMA_ThickAvg.csv ${rdir}/
mv CorticalMeasuresENIGMA_SurfAvg.csv ${rdir}/

y entonces podemos ejecutar desde cualquier parte e indicando el directorio de trabajo,

[osotolongo@brick03 mri_face]$ qc_scripts/extract.sh mriface_fsresults

y obtenemos,

[osotolongo@brick03 mri_face]$ ls CorticalMeasuresENIGMA*.csv
CorticalMeasuresENIGMA_SurfAvg.csv
CorticalMeasuresENIGMA_ThickAvg.csv

que seran nuestros archivos de trabajo para detectar los sujetos problematicos

outliers

Pues entonces vamos a buscar los outliers dentro de todos estos datos. Lo que hacemos es copir el script outliers.R en el directorio de trabajo,

[osotolongo@brick03 mri_face]$ cp qc_scripts/outliers.R ./

y ejecutar,

[osotolongo@brick03 mri_face]$ R --no-save --slave < outliers.R > outliers.log

que nos da un archivo bastante largo con los vaores que no encajan dentro de una distribucion normal,

[osotolongo@brick03 mri_face]$ head outliers.log
Please check the following subjects closely using the QC methods provided in Step 2 of the cortical protocols.
Subject XNAT05_E00042 has THICKNESS values too HIGH for the following structures: L_cuneus_thickavg
 Subject XNAT05_E00042 has THICKNESS values too HIGH for the following structures: R_cuneus_thickavg
 Subject XNAT05_E00042 has THICKNESS values too HIGH for the following structures: R_postcentral_thickavg
Subject XNAT05_E00047 has THICKNESS values too LOW for the following structures: L_caudalanteriorcingulate_thickavg
 Subject XNAT05_E00047 has THICKNESS values too LOW for the following structures: L_posteriorcingulate_thickavg
 Subject XNAT05_E00047 has THICKNESS values too LOW for the following structures: L_rostralanteriorcingulate_thickavg
 Subject XNAT05_E00047 has THICKNESS values too LOW for the following structures: R_posteriorcingulate_thickavg
Subject XNAT05_E00050 has THICKNESS values too HIGH for the following structures: L_superiorparietal_thickavg
Subject XNAT05_E00051 has THICKNESS values too LOW for the following structures: R_insula_thickavg

Ahora vamos a ver quien esta aqui,

[osotolongo@brick03 mri_face]$ grep "THICKNESS\|\|AREA" outliers.log | awk '{print $2}' | sort | uniq > outliers.list
[osotolongo@brick03 mri_face]$ wc -l outliers.list 
289 outliers.list

Y estos son los sujetos que tenemos que mirar. Aquellos que estan fuera de esta lista tienen todas las medidas dentro de la norma y no deberian tener problemas con la segmentacion.

So, lo primero que voy a hacer es marcar estos como OK poniendo este archivo como un resource en XNAT,

[osotolongo@brick03 mri_face]$ cat rating.json 
{"ResultSet":{"Result":[{"rating":"good","notes":"None"}]}

Asi que hacemos un code refactorization,

Este si que es sencillo

check_good_qc.pl

#!/usr/bin/perl
 
# Copyright 2023 O. Sotolongo <asqwerty@gmail.com>
 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
 
# Este es un script de apoyo para el QC de FS. Si ya tengo los outliers,
# puedo marcar los que no lo son como buenos en XNAT
 
use strict; use warnings;
use NEURO4 qw(load_project);
use XNATACE qw(xget_session xcreate_res xput_res_data);
 
my $prj; my $xprj;
my $idir;
my $ifile = "outliers.log";
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
	$_ = shift;
	last if /^--$/;
	if (/^-o/) { $ifile = shift; chomp($ifile);} #archivo de outliers
	if (/^-d/) { $idir = shift; chomp($idir);} # directorio con los datos de FS
	if (/^-x/) { $xprj = shift; chomp($xprj);} #nombre del proyecto en XNAT
	if (/^-p/) { $prj = shift; chomp($prj);} #nombre local del proyecto
}
 
if ($prj and not $xprj) {
	my %pdata = load_project($prj);
	$xprj = $pdata{'XNAME'};
}
 
$idir = $prj.'_fsresults' unless $idir;
die "Should supply XNAT project name or define it at local project config!\n" unless $xprj;
my %xconf = xget_session();
 
my %good_boy = ("rating" => "good", "notes" => "None");
# put outliers file into array buffer 
open IDF, $ifile or die "No such file as $ifile";
chomp(my @outliers = <IDF>);
close IDF;
 
# put FS dir into array
opendir IDD, $idir or die "Couldn't open dir $idir";
my @experiments = readdir IDD;
close IDD;
 
# grep subjects into outliers array 
foreach my $experiment (@experiments) {
	unless (grep /$experiment/, @outliers) {
		xcreate_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc');
		xput_res_data($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc', 'rating.json', \%good_boy);
		print "$experiment\n";
	}
}

y podemos subir rapidamente un archivo json compatible con los QC antiguos.

External Surface Method

VM

Ahora vamos a ejecutar el primer metodo de QC para producir las imagenes de la segmentacion cortical. Este metodo necesita ejecutar los visores de Freesurfer (con GLX), que no funcionan por VNC directamente al server. Asi qeu vamos a eejcutarlo en una VM. Asi que la VM debe estar ejecutandose,

[root@brick03 ~]# virsh list
 Id   Name        State
---------------------------
 1    FSLVm7_64   running

y la redireccion de puertos para el no_vnc activa,

[root@brick03 ~]# ./novnc_ws.sh 
WebSocket server settings:
  - Listen on :6080
  - Web server. Web root: /usr/share/novnc
  - SSL/TLS support
  - Backgrounding (daemon)

y desde la maquina local ha de hacerse un tunel ssh al 6080,

alias vnc_tunnel="~/vpn_ace.sh; ssh -C -L 6080:localhost:6080 osotolongo@172.26.2.43"

Tras esto, se supone que accediendo a la url http://localhost:6080/vnc.html tenemos acceso a la maquina virtual,

Ejecutando

Nota: Antes de ejecutar el script para preparar las imagenes en la VM, debemos crear el directorio fsqcdir y asegurarnos que el propietario es el user de la VM.

Por alguna razon hay prblemas para salvar algunas imagenes en baja resolucion. Aprovechando que tengo problemas con la visualizacion de los TIF, voy a editar el script original de ENIGMA,

En la pagina web voy a utilizar PNGs y no TIFFs porque no puedo verlos con Firefox
Solo voy a crear las imagenes TIFF de alta resolucion con tksurfer
Voy a usar ImageMagick (implica que lo necesitas instalado) para convertir las imagenes TIFF a formato PNG, tanto en alta como en baja resolucion

El script modificado queda entonces,

fsqc.sh

#!/bin/bash
idir=$1
[[ -z ${idir} ]] && { echo "Should supply input directory"; exit 1; }
[[ -d ${idir} ]] || { echo "Directory not found"; exit 404; }
 
export SUBJECTS_DIR=$(realpath ${idir})
 
[[ -d fsqcdir ]] || { echo "you chould make a 777 fsqcdir directory first"; exit 1;}
cd fsqcdir
 
echo "<html>" >  index.html
echo "<table>" >> index.html
for sub in $(ls ${SUBJECTS_DIR});
do
        echo "<tr>" >> index.html
        echo "<td><a href=\"file:"$sub".lh.lat.hr.png\"><img src=\""$sub".lh.lat.lr.png\"></a></td>" >> index.html
        echo "<td><a href=\"file:"$sub".lh.med.hr.png\"><img src=\""$sub".lh.med.lr.png\"></a></td>" >> index.html
        echo "<td><a href=\"file:"$sub".rh.lat.hr.png\"><img src=\""$sub".rh.lat.lr.png\"></a></td>" >> index.html
        echo "<td><a href=\"file:"$sub".rh.med.hr.png\"><img src=\""$sub".rh.med.lr.png\"></a></td>" >> index.html
        echo "</tr>" >> index.html
        echo "<tr>"  >> index.html
        echo "<td colspan=4><center>"$sub"</center><br></td>" >> index.html
        echo "</tr>" >> index.html
 
 
        echo "labl_import_annotation \"aparc.annot\"" > tmp.tcl
        echo "scale_brain 1.35" >>tmp.tcl
        echo "redraw" >>tmp.tcl
        echo "save_tiff "$sub".lh.lat.hr.tif" >>tmp.tcl
        echo "rotate_brain_y 180.0"  >>tmp.tcl
        echo "redraw" >>tmp.tcl
        echo "save_tiff "$sub".lh.med.hr.tif" >>tmp.tcl
        echo "exit 0" >>tmp.tcl
        tksurfer $sub lh pial -tcl tmp.tcl
 
        echo "labl_import_annotation \"aparc.annot\"" > tmp.tcl
        echo "scale_brain 1.35" >>tmp.tcl
        echo "redraw" >>tmp.tcl
        echo "save_tiff "$sub".rh.lat.hr.tif" >>tmp.tcl
        echo "rotate_brain_y 180.0" >>tmp.tcl
        echo "redraw" >>tmp.tcl
        echo "save_tiff "$sub".rh.med.hr.tif" >>tmp.tcl
        echo "rotate_brain_y -180.0" >>tmp.tcl
        echo "exit 0" >>tmp.tcl
        tksurfer $sub rh pial -tcl tmp.tcl
 
        rm tmp.tcl
# Here I use image magick, so you need to install it    
        convert ${sub}.lh.lat.hr.tif ${sub}.lh.lat.hr.png
        convert ${sub}.lh.lat.hr.tif -resize 300x300 ${sub}.lh.lat.lr.png
        convert ${sub}.lh.med.hr.tif ${sub}.lh.med.hr.png
        convert ${sub}.lh.med.hr.tif -resize 300x300 ${sub}.lh.med.lr.png
        convert ${sub}.rh.lat.hr.tif ${sub}.rh.lat.hr.png
        convert ${sub}.rh.lat.hr.tif -resize 300x300 ${sub}.rh.lat.lr.png
        convert ${sub}.rh.med.hr.tif ${sub}.rh.med.hr.png
        convert ${sub}.rh.med.hr.tif -resize 300x300 ${sub}.rh.med.lr.png
 
done;
echo "</table>" >> index.html
echo "</html>"  >> index.html

Y, en la VM, hacemos algo como,

[fsluser@fslvm7_64 mri_face] qc_scripts/fsqc.sh mriface_fsresults/

y voila, se lanza el procesamiento de los directorios de FS,

el output final de todo este laaaargo proceso es la pagina fsqcdir/index.html, y, mas importante aun, las imagenes para cada experimento, por ejemplo,

[osotolongo@brick03 mri_face]$ ls fsqcdir/XNAT05_E00065*
fsqcdir/XNAT05_E00065.lh.lat.hr.png  fsqcdir/XNAT05_E00065.rh.lat.hr.png
fsqcdir/XNAT05_E00065.lh.lat.hr.tif  fsqcdir/XNAT05_E00065.rh.lat.hr.tif
fsqcdir/XNAT05_E00065.lh.lat.lr.png  fsqcdir/XNAT05_E00065.rh.lat.lr.png
fsqcdir/XNAT05_E00065.lh.med.hr.png  fsqcdir/XNAT05_E00065.rh.med.hr.png
fsqcdir/XNAT05_E00065.lh.med.hr.tif  fsqcdir/XNAT05_E00065.rh.med.hr.tif
fsqcdir/XNAT05_E00065.lh.med.lr.png  fsqcdir/XNAT05_E00065.rh.med.lr.png

Al final de todo basta con un,

~ # firefox fsqcdir/index.html

Internal Surface Method

Para producir las imagenes restantes necesitamos,

Y ahora hemos de acomodar y preparar todo. El paquete de ENIGMA lo dejo en /old_nas/mri_face/ENIGMA_scripts/, y edito el archivo de matlab para reflejar la ubicacion de mis imagenes y archivos.

FS_directory='/old_nas/mri_face/mriface_fsresults/';
QC_output_directory='/old_nas/mri_face/ENIGMA_QC/';
ENIGMA_QC_folder='/old_nas/mri_face/ENIGMA_scripts/';

Ahora ya puedo abrir MATLAB. Aqui selecciono los directorios /old_nas/mri_face/ENIGMA_scripts/ y /old_nas/mri_face/qc_scripts/ y los añado al path con click derecho. Este ultimo directorio es donde tengo el script que voy a correr (cortical_qc.m). Y entonces voy al prompt de matlab y escribo,

>> cortical_qc

y si todo va bien me tengo que esperar un rato largo a que esten todas las imagenes para evaluar el QC.

Al final de este proceso hemos de construir la pagina web para poder mirar los slices de QC,

[osotolongo@brick03 mri_face]$ chmod +x ENIGMA_scripts/make_ENIGMA_QC_webpage.sh 
[osotolongo@brick03 mri_face]$ ENIGMA_scripts/make_ENIGMA_QC_webpage.sh /old_nas/mri_face/ENIGMA_QC/

Uploading QC results and images

Ahora, lo primero que se ha de hacer es guardar las imagenes para evaluar el QC cmo un recurso del experimento en XNAT.

Si reutilizo las funciones de mi biblioteca esto puede hacerse con relativa facilidad,

put_qc_images.pl

#!/usr/bin/perl
 
# Copyright 2023 O. Sotolongo <asqwerty@gmail.com>
 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
 
# Script purpose is to put the resulting images of ENIGMA QC into
# XNAT as eperiment resources
 
# We are going to just read the contents of a directory 
# (two of them really), pack the rigth images and put the gziped
# file inside the fsqc resource,
 
 
use strict; use warnings;
use NEURO4 qw(load_project);
use XNATACE qw(xget_session xcreate_res xput_res_file);
use Cwd qw(cwd);
use Data::Dump qw(dump);
use File::Basename;
use File::Temp qw(:mktemp tempdir);
 
my $prj; my $xprj;
my $edir; my $idir;
my $ifile = "outliers.list";
my $tmp_dir = $ENV{TMPDIR};
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
	$_ = shift;
	last if /^--$/;
	if (/^-e/) { $edir = shift; chomp($edir);} # directorio con el QC de FS
	if (/^-q/) { $idir = shift; chomp($idir);} # directorio con el QC de MATLAB
	if (/^-x/) { $xprj = shift; chomp($xprj);} #nombre del proyecto en XNAT
	if (/^-p/) { $prj = shift; chomp($prj);} #nombre local del proyecto
}
 
if ($prj and not $xprj) {
	my %pdata = load_project($prj);
	$xprj = $pdata{'XNAME'};
}
die "Should supply XNAT project name or define it at local project config!\n" unless $xprj;
 
$idir = cwd()."/ENIGMA_QC" unless $idir;
$edir = cwd()."/fsqcdir" unless $edir;
 
my %xconf = xget_session();
 
# No voy a complicarme mucho, voy a basar esto en el contenido del directorio de matlab
my @experiments = map { basename($_) } glob($idir.'/XNAT*');
foreach my $experiment (@experiments) {
	my $selection = $idir.'/'.$experiment.' '.$edir.'/'.$experiment.'.*';
	my $tzf = mktemp($tmp_dir.'/enigmaqc.XXXXX');
	my $pre = 'tar czf '.$tzf.' '.$selection;
	print "$experiment\n";
	system($pre);
	xcreate_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc');
	xput_res_file($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc', 'enigma_qc_images.tar.gz', $tzf);
}

y queda un archivo compactado con todas las imagenes de ambos QC,

The rest of them

Ahora, solo quedan 289 MRI por evaluar, lo que haré es construir una web conjunta para evaluar el QC de estas imagenes. Entonces , leyendo el archivo outliers.list puedo escoger solo las imagenes que me faltan por evaluar,

El script es bastante sencillo,

make_fsqc_web.html

 #!/bin/bash
list='outliers.list'
if [ -f fsqc.csv ]; then rm fsqc.csv; fi
echo "<html><body>" > fsqc.html
while read -r line; do
        echo "${line},fail,None" >> fsqc.csv
        echo "<table cellspacing=\"1\" style=\"width:100%;background-color:#000;\"><tr>" >> fsqc.html
        echo "<td colspan=4><font color=\"white\" size=\"5\"><center>${line}</center><br></font></td>" >> fsqc.html
        echo "</tr>" >> fsqc.html
        echo "<tr>" >> fsqc.html
        echo "<td><a href=\"file:fsqcdir/${line}.lh.lat.hr.png\"><img src=\"fsqcdir/${line}.lh.lat.lr.png\"></a></td>" >> fsqc.html
        echo "<td><a href=\"file:fsqcdir/${line}.lh.med.hr.png\"><img src=\"fsqcdir/${line}.lh.med.lr.png\"></a></td>" >> fsqc.html
        echo "<td><a href=\"file:fsqcdir/${line}.rh.lat.hr.png\"><img src=\"fsqcdir/${line}.rh.lat.lr.png\"></a></td>" >> fsqc.html
        echo "<td><a href=\"file:fsqcdir/${line}.rh.med.hr.png\"><img src=\"fsqcdir/${line}.rh.med.lr.png\"></a></td>" >> fsqc.html
        echo "</tr>" >> fsqc.html
        echo "<tr>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_10.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_10.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_20.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_20.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_40.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_40.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_50.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_50.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "</tr> <tr>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_60.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_60.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_80.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_80.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_85.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_85.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Axial_90.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Axial_70_90.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "</tr> <tr>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_10.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_10.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_20.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_20.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_40.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_40.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_50.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_50.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "</tr> <tr>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_60.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_60.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_80.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_80.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_85.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_85.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "<td><a href=\"file:ENIGMA_QC/${line}/Cortical_set_Coronal_90.png\"><img src=\"ENIGMA_QC/${line}/Cortical_set_Coronal_70_90.png\" width=\"100%\" ></a></td>" >> fsqc.html
        echo "</tr> </table>" >> fsqc.html
done < ${list}
echo "</body></html>" >> fsqc.html

y una vez ejecutado,

[osotolongo@brick03 mri_face]$ ./make_fsqc_web.sh

nos quedan, la pagina web para revisar,

y un archivo CSV, donde debemos escribir los resultados,

[osotolongo@brick03 mri_face]$ head fsqc.csv
XNAT05_E00042,fail,None
XNAT05_E00043,fail,None
XNAT05_E00047,fail,None
XNAT05_E00049,fail,None
XNAT05_E00050,fail,None
XNAT05_E00051,fail,None
XNAT05_E00057,fail,None
XNAT05_E00059,fail,None
XNAT05_E00060,fail,None
XNAT05_E00062,fail,None

Por defecto, es solo una lista de los outliers, indicando la MRI como fail y las notas a ecribir como None. Durante la posterior revision ha de editarse,

[osotolongo@brick03 mri_face]$ head fsqc.csv
XNAT05_E00042,good,None
XNAT05_E00043,good,None
XNAT05_E00047,fail,inner fail, too much lessions?
XNAT05_E00049,good,None
XNAT05_E00050,fail,None
XNAT05_E00051,good,None
XNAT05_E00057,good,None
XNAT05_E00059,good,None
XNAT05_E00060,good,None
XNAT05_E00062,good,None

Y una vez revisados todos los outliers y editado el CSV, hay que subir el json correspondiente al recurso del experimento.

El script para subir los resultados es muy similar al de los no outliers

check_bad_qc.pl

#!/usr/bin/perl
 
# Copyright 2023 O. Sotolongo <asqwerty@gmail.com>
 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
 
# Este es un script de apoyo para el QC de FS. Si ya tengo los outliers,
# puedo marcar los que no lo son como buenos en XNAT
 
use strict; use warnings;
use NEURO4 qw(load_project);
use XNATACE qw(xget_session xcreate_res xput_res_data);
 
my $prj; my $xprj;
my $idir;
my $ifile = "fsqc.csv";
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
        $_ = shift;
        last if /^--$/;
        if (/^-i/) { $ifile = shift; chomp($ifile);} #archivo de outliers
        if (/^-x/) { $xprj = shift; chomp($xprj);} #nombre del proyecto en XNAT
        if (/^-p/) { $prj = shift; chomp($prj);} #nombre local del proyecto
}
 
if ($prj and not $xprj) {
        my %pdata = load_project($prj);
        $xprj = $pdata{'XNAME'};
}
 
die "Should supply XNAT project name or define it at local project config!\n" unless $xprj;
my %xconf = xget_session();
 
# read fsqc results fron input file
open IDF, $ifile or die "No such file as $ifile";
while (<IDF>) {
        my ($experiment, $rating, $notes) = /^(.*),(.*),(.*)$/;
        if ($experiment and $rating) {
                $notes = 'None' unless $notes;
                my %bad_boy = ("rating" => "$rating", "notes" => "$notes");
                xcreate_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc');
                xput_res_data($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'fsqc', 'rating.json', \%bad_boy);
                print "$experiment, $rating : $notes\n";
        }
}
close IDF;

y la ejecucion igual de sencilla,

[osotolongo@brick03 mri_face]$ ./check_bad_qc.pl -p mriface

Los datos deberian quedar entonces almacenados como recurso para todos los experimentos y disponibles para bajar, revisar y manipular posteriormente,

Bajando los resutados a una tabla

No pain. Lo hemos hecho backwards compatible.

[osotolongo@brick03 mri_face]$ xnat_get_fsqc.pl -x mriface

[osotolongo@brick03 mri_face]$ head mriface_fsqc_data.csv
Subject,Date,FSQC,Notes
19990095,2022-10-18,Good,None
20040526,2022-07-05,Good,None
20050456,2022-03-06,Good,None
20050604,2022-09-02,Good,None
20070391,2022-10-25,Good,None
20070711,2023-02-28,Good,None
20070841,2023-01-25,Good,None
20071018,2022-04-10,Good,None
20080792,2022-11-04,Good,None

Detritus Wiki

Table of Contents