MOPEAD

Conversion DICOM (full)

Hay nu cambio de protocolo de transmision de los DCMs. Asi que vamos a convertir el DCM completo y luego tomar los nifti que necesitamos y ordenarlos. Algo similar a como hariamos en formato BIDS pero para nuestro pipeline. Tenemos el protocolo, conocemos la estructura de archivos, sabemos que problemas puede haber, asi que no deberia ser complicado.

MRI

Pruebas iniciales

Las imágenes se suben por FTP desde Corachan. De momento yo voy a procesar las T1. Voy a intentar una primero a ver que tal se convierte.

enga que pa luego es tarde

$ cd /nas/osotolongo/mopead/raw
$ mkdir tmp
$ for x in `find /nas/corachan/MOPEAD/52JAYALW -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep t1_mprage` ]]; then dcm2nii -o tmp/ $x; fi; done
Chris Rorden's dcm2nii :: 4AUGUST2014 64bit BSD License
reading preferences file /home/osotolongo/.dcm2nii/dcm2nii.ini
Data will be exported to tmp/
Validating 1 potential DICOM images.
Found 1 DICOM images.
Converting 1/1  volumes: 1
00057101->s005a001.nii
 brightest voxel was 1015: data will be saved as 16-bit signed integer.
GZip...s005a001.nii.gz
$ fslview /nas/osotolongo/mopead/raw/tmp/s005a001.nii.gz &

y todo parece OK

Así que vamos a intentar hacerlo en bulk

$ make_std.pl mopead
$ cd /nas/osotolongo/data/mopead
$ mkdir tmp
$ mkdir converted
$ for x in /nas/corachan/MOPEAD/*; do for y in `find $x -type f`; do if [[ `dckey -k "SeriesDescription" $y 2>&1 | grep t1_mprage` ]]; then dcm2nii -o tmp/ $y; mv tmp/*.nii.gz converted/${x}_t1.nii.gz; fi; done; done
$ ls converted/
52JAYALW_t1.nii.gz  AAVEYS8L_t1.nii.gz  AVTNFTVE_t1.nii.gz  D2A774IW_t1.nii.gz  PI8NVHLE_t1.nii.gz

Y parece que funciona. La cuestion ahora es pasar estos archivos a la convencion de nombres del pipeline de neuroimagen. Para el proyecto MOPEAD vamos a usar mop0000s0005.nii.gz para las MRI.

DB y preprocesamiento T1

El primer paso seria entonces en construir y mantener automaticamente una DB de las imagenes y al mismo tiempo asignarles nuevos nombres. Como se supone que el FTP de Corachan sea de solo lectura (Es decir, no se borraran imagenes), voy a tomarlo como base para construir la lista. De aqui creo un fichero de IDs y lo actualizo cada vez que sea necesario.

Hasta ahi es bastante simple

El siguiente paso es convertir todo a nifti. para ello voy a hacer un directorio en $std{'DATA'}./niftis/ y convierto ahi todo.

Basicamente correr dcm2niix

Ahora toca buscar y extraer las T1, leer la DB y copiarlas a mri/ .

Tampoco parece muy dificil

update_t1.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study);
use Data::Dump qw(dump);
use File::Find::Rule;
 
my $study = "mopead";
my %std = load_study($study);
my $src_dir = '/nas/corachan/MOPEAD/';
my $conv_file = $std{'DATA'}.'/converted_t1.csv';
my $ids_file = $std{'DATA'}.'/ids.csv';
my $names_file =$std{'DATA'}.'/ni_names.csv';
my $proj_file = $std{'DATA'}.'/mopead.csv';
my @other_exts = ("nii.gz","json");
my $series = 5;
my $tag = 't1_mprage_sag_p2_iso';
 
#Leo la DB
my %idsinfo;
open IDS, "<$ids_file" or die $!;
while(<IDS>){
        chomp;
        my ($key, $value) = split(/,\s?/);
        $idsinfo{$key} = $value;
}
close IDS;
#Leo los que ya he descomprimido
my @orig_str = qx/ls $src_dir/;
chomp @orig_str;
#Leo los que se han convertido
my @conv_str;
open CONV, "<$conv_file" or die $!;
chomp (@conv_str = <CONV>);
close CONV;
#Y a ver cuantos faltan por convertir!
my %conv_str = map { $_ => 1 } @conv_str;
my @not_conv = grep !$conv_str{$_}, @orig_str;
# y a convertirlos!
foreach my $guy (@not_conv){
        if(exists($idsinfo{$guy})){
                my @conv_files = find(file => 'name' => "*$tag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                my $count = $series;
                foreach my $nii_file (@conv_files){
                        my $order = 'fslval '.$nii_file.' dim1';
                        my $xinfo = qx/$order/;
                        if ($xinfo == '176'){
                                my $file_name = $nii_file;
                                foreach my $ext (@other_exts){
                                        ($file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                        $order = sprintf "cp %s %s/mop%04ds%04d.%s", $file_name, $std{'MRI'}, $idsinfo{$guy}, $count, $ext;
                                        print "$order\n";
                                        system($order);
                                }
                        $count++;
                        }
                }
        }
}
#actualizo los convertidos
push(@conv_str, @not_conv);
open CONV, ">$conv_file" or die $!;
print CONV "$_\n" foreach @conv_str;
close CONV;

So far, so good. Ya la estructura estaría preparada para importar los sujetos en freesurfer.

$ pfsl2fs.pl mopead
$ ls -l /home/data/subjects/mopead_mop00*/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10541099 Sep  7 10:35 /home/data/subjects/mopead_mop0001/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11368986 Sep  7 10:35 /home/data/subjects/mopead_mop0002/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10591570 Sep  7 10:35 /home/data/subjects/mopead_mop0003/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10859697 Sep  7 10:35 /home/data/subjects/mopead_mop0004/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10492571 Sep  7 10:35 /home/data/subjects/mopead_mop0005/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11305758 Sep  7 10:35 /home/data/subjects/mopead_mop0006/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10985299 Sep  7 10:35 /home/data/subjects/mopead_mop0007/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 12155077 Sep  7 10:35 /home/data/subjects/mopead_mop0008/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 12134606 Sep  7 10:35 /home/data/subjects/mopead_mop0009/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10235152 Sep  7 10:35 /home/data/subjects/mopead_mop0010/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 12949235 Sep  7 10:35 /home/data/subjects/mopead_mop0011/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11922977 Sep  7 10:35 /home/data/subjects/mopead_mop0012/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11349419 Sep  7 10:35 /home/data/subjects/mopead_mop0013/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10657028 Sep  7 10:35 /home/data/subjects/mopead_mop0014/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10488166 Sep  7 10:35 /home/data/subjects/mopead_mop0015/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11327625 Sep  7 10:35 /home/data/subjects/mopead_mop0016/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11650366 Sep  7 10:35 /home/data/subjects/mopead_mop0017/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10562121 Sep  7 10:35 /home/data/subjects/mopead_mop0018/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10545429 Sep  7 10:35 /home/data/subjects/mopead_mop0019/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10924535 Sep  7 10:35 /home/data/subjects/mopead_mop0020/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 11250253 Sep  7 10:35 /home/data/subjects/mopead_mop0021/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10130792 Sep  7 10:35 /home/data/subjects/mopead_mop0022/mri/orig/001.mgz
-rw-rw-r-- 1 osotolongo osotolongo 10216987 Sep  7 10:35 /home/data/subjects/mopead_mop0023/mri/orig/001.mgz

Preprocesamiento T2

Esto, ya hecho la conversión de los nifti y el script para T1 es bastante sencillo. La unica diferencia es que hay que buscar la etiqueta t2_space_dark-fluid_sag_p2_iso y guardar las imagenes con otro nombre, en un directorio distinto. voy a esribirlos como ….s0006.nii.gz, dentro de mri2. Despues se puede cambiar facilmente.

El codigo aqui

update_t2.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study);
use Data::Dump qw(dump);
use File::Find::Rule;
 
my $study = "mopead";
my %std = load_study($study);
my $src_dir = '/nas/corachan/MOPEAD/';
my $conv_file = $std{'DATA'}.'/converted_t2.csv';
my $ids_file = $std{'DATA'}.'/ids.csv';
my $names_file =$std{'DATA'}.'/ni_names.csv';
my $proj_file = $std{'DATA'}.'/mopead.csv';
my @other_exts = ("nii.gz","json");
my $series = 6;
my $tag = "t2_space_dark-fluid_sag_p2_iso";
 
#Leo la DB
my %idsinfo;
open IDS, "<$ids_file" or die $!;
while(<IDS>){
        chomp;
        my ($key, $value) = split(/,\s?/);
        $idsinfo{$key} = $value;
}
close IDS;
#Leo los que ya he descomprimido
my @orig_str = qx/ls $src_dir/;
chomp @orig_str;
#Leo los que se han convertido
my @conv_str;
open CONV, "<$conv_file" or die $!;
chomp (@conv_str = <CONV>);
close CONV;
#Y a ver cuantos faltan por convertir!
my %conv_str = map { $_ => 1 } @conv_str;
my @not_conv = grep !$conv_str{$_}, @orig_str;
# y a convertirlos!
foreach my $guy (@not_conv){
        if(exists($idsinfo{$guy})){
                my @conv_files = find(file => 'name' => "*$tag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');       
                my $count = $series;
                foreach my $nii_file (@conv_files){
                        my $order = 'fslval '.$nii_file.' dim1';
                        my $xinfo = qx/$order/;
                        if ($xinfo == '192'){
                                my $file_name = $nii_file;
                                foreach my $ext (@other_exts){
                                        ($file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                        my $order = sprintf "cp %s %s/mop%04ds%04d.%s", $file_name, $std{'MRI2'}, $idsinfo{$guy}, $count, $ext;
                                        print "$order\n";
                                        system($order);
                                }
                                $count++;
                        }
                }
        }
}
#actualizo los convertidos
push(@conv_str, @not_conv);
open CONV, ">$conv_file" or die $!;
print CONV "$_\n" foreach @conv_str;
close CONV;

Paralelizando Freesurfer

Una nota sobre el cluster

El schedule manager que utiliza el cluster es slurm. El cluster se ha reconfigurado de tal manera que utiliza las CPU como recursos consumibles a la hora de asignar tareas. Tambien se ha creado un recurso consumible de GPU que puede utilizarse. Hay dos particiones posibles, devel con todas las maquinas y cuda con las maquinas que incluyen tarjetas GPU. En total tenemos 64*3+32 CPU y 3 slots GPU.

En algunos casos (como con bedpostx, eddy o probtrackx) vale la pena utilizar la GPU pues la velocidad de los calculos es mas de 100 veces superior o las funciones que estan implementadas en la GPU no existen en la version de CPU. En otros casos (como Freesurfer) la ventaja es dudosa, ya que la velocidad de procesamiento es solo unas 10 veces mayor. Aunque la GPU podria sobrecargarse con ejecuciones recurrentes no es aconsejable ya que es dificil de controlar a priori cualquier desbordamiento de memoria o algo parecido. Asi, que en algunos casos se ignorara la GPU y en otros casos se lanzaran solo 3 tareas simultaneas.

En cualquier caso, cada tarea se lanzara por separado y solo se enviara un email cuando alguna falle.Para enviar un email cuando todas las tareas terminan he creado un script sencillo,

emailme.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study breciahi);
 
my $study = shift;
my $script = shift;
die "Must supply study name!\n" unless $study;
 
breciahi($study, $script);

esto habria que lanzarlo con un script para sbatch, y habria que insertarlo al final del codigo de cualquier script,

my $sjobs = join(',',@jobs_list);
my $orderfile = $outdir.'/fs_recon_end.sh';
open ORD, ">$orderfile";
print ORD '#!/bin/bash'."\n";
print ORD '#SBATCH -J fs_recon_end'."\n";
print ORD '#SBATCH -o '.$outdir.'/fs_recon_end-%j'."\n";
print ORD 'srun —-nodelist=detritus '.$ENV{'PIPEDIR'}.'/bin/emailme.pl '.$study.' '.basename($ENV{_})."\n";
close ORD;
my $order = 'sbatch --depend=afterany:'.$sjobs.' '.$orderfile;
exec($order);

Nota: Siempre habria que guardar en un array los JobID de todos los procesos que se han lanzado para pasarlos como dependencias.

Ahora si, Freesurfer

La idea general recorrer la lista de los sujetos, lanzar recon-all como una tarea de slurm y guardar el jobid,

my @jobs_list;
foreach my $pkey (sort keys %plist){
	my $subj = $study."_".$plist{$pkey}.$pkey;
	my $gpu;
	if(check_subj($subj)){
		$gpu = ""; #because it dont use gpu, YET!
		my $order = "recon-all ".$stage." -subjid ".$subj." ".$gpu;
		print CNF "$order\n";
		my $orderfile = $outdir.'/'.$subj.'fs_orders.sh';
		open ORD, ">$orderfile";
		print ORD '#!/bin/bash'."\n";
		print ORD '#SBATCH -J fs_recon_'.$subj."\n";
		print ORD '#SBATCH --time=72:0:0'."\n"; #si no ha terminado en X horas matalo
		print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
		print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
		print ORD '#SBATCH -o '.$outdir.'/fs_recon-slurm-'.$subj.'-%j'."\n";
		print ORD "srun $order\n";
                close ORD;
		my $jobid = `sbatch $orderfile`;
		$jobid = ( split ' ', $jobid )[ -1 ];
		push @jobs_list, $jobid;
		$debug ? print DBG "$order\n" :0;
	}
}

Aqui el codigo completo

precon.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Basename qw(basename);
use Data::Dump qw(dump);
 
use NEURO qw(populate check_subj load_study print_help achtung shit_done get_pair check_or_make);
 
my $attach = 0;
my $cfile;
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-a/) { $attach = 1;} 
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/recon.hlp'; exit;}
}
 
my $study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/recon.hlp'; exit;}
my $debug = 1;
my $max_processes = 24;                                                                                                                                                  
my $max_gpu = 12;
my $nodes = 3;
my $stage = "-autorecon-all";
y %std = load_study($study);
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
#open debug file
my $logfile = "$std{'DATA'}/.debug.controlled.log";
$debug ? open DBG, ">$logfile" :0;
#open slurm file
my $orderfile = "$std{'DATA'}/mri_orders.sh";
my $conffile = "$std{'DATA'}/mri_orders.conf";
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
 
#get subjects from database or file
my %plist;
my %iplist = populate($db);
 
 
if ($cfile){
        my %cuts = get_pair($data_dir."/".$cfile);
        foreach my $cut (keys %cuts){
                if(grep {/$cut/} %iplist){
                        $plist{$cut} = $iplist{$cut};
                }
        }
}else{
        %plist = %iplist;
}
my @jobs_list;
foreach my $pkey (sort keys %plist){
        my $subj = $study."_".$plist{$pkey}.$pkey;
        my $care; my $gpu;
        if(check_subj($subj)){
                $gpu = ""; #because it dont use gpu, YET!
                my $order = "recon-all ".$stage." -subjid ".$subj." ".$gpu;
                my $orderfile = $outdir.'/'.$subj.'fs_orders.sh';
                open ORD, ">$orderfile";
                print ORD '#!/bin/bash'."\n";
                print ORD '#SBATCH -J fs_recon_'.$subj."\n";
                print ORD '#SBATCH --time=72:0:0'."\n"; #si no ha terminado en X horas matalo
                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
                print ORD '#SBATCH -o '.$outdir.'/fs_recon-slurm-'.$subj.'-%j'."\n";
                print ORD "srun $order\n";
                close ORD;
                my $jobid = `sbatch $orderfile`;
                $jobid = ( split ' ', $jobid )[ -1 ];
                push @jobs_list, $jobid;
                $debug ? print DBG "$order\n" :0;
        }
}
$debug ? close DBG:0;   
my $sjobs = join(',',@jobs_list);
my $orderfile = $outdir.'/fs_recon_end.sh';
open ORD, ">$orderfile";
print ORD '#!/bin/bash'."\n";
print ORD '#SBATCH -J fs_recon_end'."\n";
print ORD '#SBATCH -o '.$outdir.'/fs_recon_end-%j'."\n";
print ORD 'srun -w detritus '.$ENV{'PIPEDIR'}.'/bin/emailme.pl '.$study.' '.basename($ENV{_})."\n";
close ORD;
my $order = 'sbatch --depend=afterany:'.$sjobs.' '.$orderfile;
exec($order);

Ejecuto el script y veo que los jobs se esta ejecutando,

$ precon.pl mopead
Submitted batch job 6003
 
$ squeue | grep fs_recon
              6003     devel fs_recon osotolon PD       0:00      1 (Dependency)
              5980     devel fs_recon osotolon  R       0:33      1 brick02
              5981     devel fs_recon osotolon  R       0:33      1 brick02
              5982     devel fs_recon osotolon  R       0:33      1 brick02
              5983     devel fs_recon osotolon  R       0:33      1 brick02
              5984     devel fs_recon osotolon  R       0:33      1 brick02
              5985     devel fs_recon osotolon  R       0:33      1 brick02
              5986     devel fs_recon osotolon  R       0:33      1 brick02
              5987     devel fs_recon osotolon  R       0:33      1 brick02
              5988     devel fs_recon osotolon  R       0:33      1 brick02
              5989     devel fs_recon osotolon  R       0:33      1 brick02
              5990     devel fs_recon osotolon  R       0:33      1 brick02
              5991     devel fs_recon osotolon  R       0:33      1 brick02
              5992     devel fs_recon osotolon  R       0:33      1 brick02
              5993     devel fs_recon osotolon  R       0:33      1 brick02
              5994     devel fs_recon osotolon  R       0:33      1 brick02
              5995     devel fs_recon osotolon  R       0:33      1 brick02
              5996     devel fs_recon osotolon  R       0:33      1 brick02
              5997     devel fs_recon osotolon  R       0:33      1 brick02
              5998     devel fs_recon osotolon  R       0:33      1 brick02
              5999     devel fs_recon osotolon  R       0:33      1 brick02
              6000     devel fs_recon osotolon  R       0:33      1 brick02
              6001     devel fs_recon osotolon  R       0:33      1 brick02
              6002     devel fs_recon osotolon  R       0:33      1 brick02

Nota: El job 6003, en este caso, es el email de aviso, que depende de todos los demas.

$ cat slurm/fs_recon_end.sh 
#!/bin/bash
#SBATCH -J fs_recon_end
#SBATCH -o /nas/osotolongo/data/mopead/slurm/fs_recon_end-%j
srun -w detritus /nas/software/neuro.dev/bin/emailme.pl mopead precon.pl

Y ahi lo dejo. Faltaria comprobar los resultados de freesurfer. En este caso, cuando un script falla manda un email, señalando el sujeto y el numero de la tarea.

Los logs estan en el directorio slurm/ , por lo que bastaria buscar el log correspondiente para ver que ha pasado.

$ ls slurm/fs*
slurm/fs_recon_end.sh                     slurm/fs_recon-slurm-mopead_mop0012-5991
slurm/fs_recon-slurm-mopead_mop0001-5980  slurm/fs_recon-slurm-mopead_mop0013-5992
slurm/fs_recon-slurm-mopead_mop0002-5981  slurm/fs_recon-slurm-mopead_mop0014-5993
slurm/fs_recon-slurm-mopead_mop0003-5982  slurm/fs_recon-slurm-mopead_mop0015-5994
slurm/fs_recon-slurm-mopead_mop0004-5983  slurm/fs_recon-slurm-mopead_mop0016-5995
slurm/fs_recon-slurm-mopead_mop0005-5984  slurm/fs_recon-slurm-mopead_mop0017-5996
slurm/fs_recon-slurm-mopead_mop0006-5985  slurm/fs_recon-slurm-mopead_mop0018-5997
slurm/fs_recon-slurm-mopead_mop0007-5986  slurm/fs_recon-slurm-mopead_mop0019-5998
slurm/fs_recon-slurm-mopead_mop0008-5987  slurm/fs_recon-slurm-mopead_mop0020-5999
slurm/fs_recon-slurm-mopead_mop0009-5988  slurm/fs_recon-slurm-mopead_mop0021-6000
slurm/fs_recon-slurm-mopead_mop0010-5989  slurm/fs_recon-slurm-mopead_mop0022-6001
slurm/fs_recon-slurm-mopead_mop0011-5990  slurm/fs_recon-slurm-mopead_mop0023-6002
 
$ grep ERROR slurm/fs_recon-slurm-mopead_mop0022-6001 
ERROR: talairach_afd: Talairach Transform: transforms/talairach.xfm ***FAILED*** (p=0.0195, pval=0.0000 < threshold=0.0050)
ERROR: talairach_afd: Talairach Transform: transforms/talairach.xfm ***FAILED*** (p=0.0000, pval=0.0000 < threshold=0.0050)
ERROR: talairach_afd: Talairach Transform: transforms/talairach.xfm ***FAILED*** (p=0.0195, pval=0.0000 < threshold=0.0050)
ERROR: Talairach failed!
recon-all -s mopead_mop0022 exited with ERRORS at Tue Nov 20 22:35:48 CET 2018

Nota: Por si hay que arreglar algo me dejo apuntado mi garbage collector para slurm y freesurfer

$ for x in `squeue | awk {'print $1'} | grep -v JOBID`; do scancel $x; done
$ for x in `find /nas/data/subjects/mopead* -name IsRunning.lh+rh`; do rm $x; done

DTI

Entendiendo la estructura

El DTI Viene con etiqueta DTIep2d_diff_mddw_48dir_p3_AP

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep DTIep2d_diff_mddw_48dir_p3_AP` ]]; then dirname $x; fi;done | uniq
/nas/corachan/MOPEAD/52JAYALW/00000000/00000019
/nas/corachan/MOPEAD/52JAYALW/00000000/00000020
/nas/corachan/MOPEAD/52JAYALW/00000000/00000021
/nas/corachan/MOPEAD/52JAYALW/00000000/00000022
/nas/corachan/MOPEAD/52JAYALW/00000000/00000023
/nas/corachan/MOPEAD/52JAYALW/00000000/00000024
/nas/corachan/MOPEAD/52JAYALW/00000000/00000025
/nas/corachan/MOPEAD/52JAYALW/00000000/DTI

El directorio DTI parece contener la misma informacion que el resto de los archivos, asi que lo voy a quitar de la busqueda.

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f | grep -v DTI`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep DTIep2d_diff_mddw_48dir_p3_AP` ]]; then dirname $x; fi;done | uniq
/nas/corachan/MOPEAD/52JAYALW/00000000/00000019
/nas/corachan/MOPEAD/52JAYALW/00000000/00000020
/nas/corachan/MOPEAD/52JAYALW/00000000/00000021
/nas/corachan/MOPEAD/52JAYALW/00000000/00000022
/nas/corachan/MOPEAD/52JAYALW/00000000/00000023
/nas/corachan/MOPEAD/52JAYALW/00000000/00000024
/nas/corachan/MOPEAD/52JAYALW/00000000/00000025

Me quedan un total de siete imagenes a convertir:

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f | grep -v DTI`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep DTIep2d_diff_mddw_48dir_p3_AP` ]]; then dirname $x; fi;done | uniq > tmp_list.dirs
$ for x in `cat tmp_list.dirs`; do dcm2niix -z y -o tmp $x; done
$ rm tmp_list.dirs
$ ls -lh tmp
total 81M
-rw-rw-r-- 1 osotolongo osotolongo  252 Sep 27 10:31 00000019_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_20.bval
-rw-rw-r-- 1 osotolongo osotolongo 1.4K Sep 27 10:31 00000019_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_20.bvec
-rw-rw-r-- 1 osotolongo osotolongo 1.3K Sep 27 10:31 00000019_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_20.json
-rw-rw-r-- 1 osotolongo osotolongo  40M Sep 27 10:31 00000019_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_20.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo  508 Sep 27 10:31 00000020_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_21.json
-rw-rw-r-- 1 osotolongo osotolongo  39M Sep 27 10:31 00000020_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_21.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo 1.3K Sep 27 10:31 00000021_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_22.json
-rw-rw-r-- 1 osotolongo osotolongo 509K Sep 27 10:31 00000021_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_22.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo 1.3K Sep 27 10:31 00000022_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_23.json
-rw-rw-r-- 1 osotolongo osotolongo 364K Sep 27 10:31 00000022_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_23.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo 1.2K Sep 27 10:31 00000023_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_24.json
-rw-rw-r-- 1 osotolongo osotolongo 462K Sep 27 10:31 00000023_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_24.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo 1.3K Sep 27 10:31 00000025_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_27.json
-rw-rw-r-- 1 osotolongo osotolongo 767K Sep 27 10:31 00000025_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_27.nii.gz

Ahora, para estimar ls distorsiones, usando topup, necesito el B0 adquirido en sentido contrario,que esta bajo la etiqueta DTIep2d_diff_mddw_4b0_PA.Luego,

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f | grep -v DTI`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep DTIep2d_diff_mddw_4b0_PA` ]]; then dcm2niix -z y -o tmp $x; fi;done
$ ls tmp
00000026_DTIep2d_diff_mddw_4b0_PA_20180705100054_28.json  00000026_DTIep2d_diff_mddw_4b0_PA_20180705100054_28.nii.gz

Escogiendo el DTI

De todos los archivos convertidos,la informacion real está en 00000019_DTIep2d_diff_mddw_48dir_p3_AP_20180705100054_20, el resto son B0s adicionales. Vamos a comenzar procesando este. La idea es seguir el mismo esquema pero ahora se ha de encontrar el archivo con dim4 mayor que 1 y copiarlo al directorio dti. Ojo, que se han de copiar ademas los archivos .bval, .bvec y .json. Nota: Ya estan todos convertidos en el directorio niftis/ .

Además, copio el B0 P»A

Aqui el codigo

update_dti.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study);
use File::Find::Rule;
use Data::Dump qw(dump);
 
my $study = "mopead";
my $prefix = "mop";
my %std = load_study($study);
my $src_dir = '/nas/corachan/MOPEAD/';
my $conv_file = $std{'DATA'}.'/converted_dti.csv';
my $ids_file = $std{'DATA'}.'/ids.csv';
my $names_file =$std{'DATA'}.'/ni_names.csv';
my $proj_file = $std{'DATA'}.'/mopead.csv';
my @other_exts = ("nii.gz","bval","bvec","json");
my @b0_exts = ("nii.gz","json");
my $tag = "DTIep2d_diff_mddw_48dir_p3_AP";
my $b0tag = "DTIep2d_diff_mddw_4b0_PA";
 
#Leo la DB
my %idsinfo;
open IDS, "<$ids_file" or die $!;
while(<IDS>){
        chomp;
        my ($key, $value) = split(/,\s?/);
        $idsinfo{$key} = $value;
}
close IDS;
#Leo los que ya he descomprimido
my @orig_str = qx/ls $src_dir/;
chomp @orig_str;
#Leo los que se han convertido
my @conv_str;
open CONV, "<$conv_file" or die $!;
chomp (@conv_str = <CONV>);
close CONV;
#Y a ver cuantos faltan por convertir!
my %conv_str = map { $_ => 1 } @conv_str;
my @not_conv = grep !$conv_str{$_}, @orig_str;
# y a convertirlos!
my @newones;
foreach my $guy (@not_conv){
        if(exists($idsinfo{$guy})){
                print "Getting DTI data A>>P\n";
                my @conv_files = find(file => 'name' => "*$tag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                foreach my $nii_file (@conv_files){
                        my $order = 'fslval '.$nii_file.' dim4';
                        my $xinfo = qx/$order/;
                        (my $checkv = $nii_file) =~ s/nii.gz$/bval/;
                        my $dest = sprintf "%s/%s%04ds0020.nii.gz", $std{'DTI'}, $prefix, $idsinfo{$guy};
                        if($xinfo>1 && -e $checkv && ! -e $dest){
                                print "Choosing and moving files\n";
                                my $file_name = $nii_file;
                                foreach my $ext (@other_exts){
                                        ($file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                        $order = sprintf "cp %s %s/%s%04ds0020.%s", $file_name, $std{'DTI'}, $prefix, $idsinfo{$guy}, $ext;
                                        print "$order\n";
                                        system($order);
                                }
                        push @newones, $guy;                    
                        }
                }
                #Convierto los B0 en sentido contrario (P>>A)
                print "Getting P>>A B0\n";
                #system($order);
                @conv_files = find(file => 'name' => "*$b0tag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                foreach my $nii_file (@conv_files){
                        foreach my $ext (@b0_exts){
                                (my $file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                my $order = sprintf "cp %s %s/%s%04d_p2a_b0.%s", $file_name, $std{'DTI'}, $prefix, $idsinfo{$guy}, $ext;
                                print "$order\n";
                                system($order);
                        }
                }
        }
}
#actualizo los convertidos
push(@conv_str, @newones);
open CONV, ">$conv_file" or die $!;
print CONV "$_\n" foreach @conv_str;
close CONV;

Preprocesamiento

El primer paso es aplicar topup. Para ello vamos a crear la estructura y los archivos correctos.

$ mkdir working/.tmp_mop0001
$ fslroi dti/mop0001s0020.nii.gz working/.tmp_mop0001/a2p_b0 0 1 
$ fslmerge -t working/.tmp_mop0001/a2p_p2a_b0 working/.tmp_mop0001/a2p_b0 working/.tmp_mop0001/mop0001_p2a_b0

Ahora corremos topup. Para ello también necesitaremos el archivo acqparams.txt que se construye con los datos de la adquisicion. Como este se supone que sea el mismo para todo el protocolo lo podemos hacer comun al proyecto. Nota: El archivo debe estar en el directorio del proyecto

$ cat acqparams.txt 
0 1 0 0.10656
0 -1 0 0.10656

topup se supone que tiene ya todo lo necesario para correr. (esto demora un ratejo)

$ topup --imain=working/.tmp_mop0001/a2p_p2a_b0 --datain=acqparams.txt --out=working/.tmp_mop0001/topup_results --iout=working/.tmp_mop0001/hifi_b0

Esta ultima imagen la aprovecharé para extraer el cerebro, (despues necesitare la mascara)

$ fslmaths working/.tmp_mop0001/hifi_b0 -Tmean working/.tmp_mop0001/hifi_b0
$ bet working/.tmp_mop0001/hifi_b0 working/.tmp_mop0001/hifi_b0_brain -m

Y ya casi está. Ahora necesito un archivo que diga a eddy que línea del acqparams.txt es la correcta. Esto se puede hacer comun al proyecto.

$ indx=""
$ for ((i=1; i<=54; i+=1)); do indx="$indx 1"; done
$ echo $indx > index.txt
$ cat index.txt 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

El proximo paso es utilizar eddy. Vamos a usar eddy_cuda, que utiliza las capacidades GPU y deberia funcionar mas rapido. Esta nueva version,utiliza directamente los datos de topup para realizar las correcciones y mueve los vectores, entre otras cosas. La sintaxis es,

${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=${out}/../dti_index.txt --bvecs=${td}/bvecs --bvals=${td}/bvals --topup=${td}/topup_results --out=${td}/data

Este paso no es rapido.

Ahora el dtifit, usando la mascara que habiamos hecho

$ dtifit --data=working/.tmp_mop0001/data --out=working/.tmp_mop0001/dti --mask=working/.tmp_mop0001/hifi_b0_brain --bvecs=working/.tmp_mop0001/bvecs --bvals=working/.tmp_mop0001/bvals
$ ls working/.tmp_mop0001/dti*
working/.tmp_mop0001/dti_FA.nii.gz  working/.tmp_mop0001/dti_MD.nii.gz  working/.tmp_mop0001/dti_V2.nii.gz
working/.tmp_mop0001/dti_L1.nii.gz  working/.tmp_mop0001/dti_MO.nii.gz  working/.tmp_mop0001/dti_V3.nii.gz
working/.tmp_mop0001/dti_L2.nii.gz  working/.tmp_mop0001/dti_S0.nii.gz
working/.tmp_mop0001/dti_L3.nii.gz  working/.tmp_mop0001/dti_V1.nii.gz

y tenemos,los vectores, la FA, MD, etc. Podemos copiarlo todo al directorio de trabajo para borrar el dichoso .tmp_XXXX

$ for x in working/.tmp_mop0001/dti*; do imcp ${x} working/mop0001_$(basename $x); done;
$ imcp working/.tmp_mop0001/data working/mop0001_dti_data
$ imcp working/.tmp_mop0001/hifi_b0_brain working/mop0001_dti_brain_mask
$ ls working/
mop0001_dti_brain_mask.nii.gz  mop0001_dti_L1.nii.gz  mop0001_dti_MD.nii.gz  mop0001_dti_V1.nii.gz
mop0001_dti_data.nii.gz        mop0001_dti_L2.nii.gz  mop0001_dti_MO.nii.gz  mop0001_dti_V2.nii.gz
mop0001_dti_FA.nii.gz          mop0001_dti_L3.nii.gz  mop0001_dti_S0.nii.gz  mop0001_dti_V3.nii.gz

dti_preproc.sh

#!/bin/sh                                                                                                                                                                
 
Usage() {                                                                                                                                                                
    echo ""                                                                                                                                                              
    echo "Usage: dti_proc.sh  <name> <T1> <input> <output_dir>"                                                                                                          
    echo ""                                                                                                                                                              
    echo "You must have FSL installed in order to run this script"                                                                                                       
    echo ""                                                                                                                                                              
    exit 1                                                                                                                                                               
}                                                                                                                                                                        
 
[ "$4" = "" ] && Usage                                                                                                                                                   
debug=1                                                                                                                                                                  
in=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                         
shift                                                                                                                                                                    
t1=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                         
shift                                                                                                                                                                    
b_in=$1                                                                                                                                                                  
shift                                                                                                                                                                    
out=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                        
shift                                                                                                                                                                    
td=${out}'/.tmp_'${b_in}                                                                                                                                                 
p2a=$(dirname ${in}.nii.gz)'/'${b_in}'_p2a_b0'                                                                                                                           
if [ ! -d "$td" ]; then                                                                                                                                                  
        mkdir $td                                                                                                                                                        
fi                                                                                                                                                                       
if [ ! -d "$out" ]; then                                                                                                                                                 
        mkdir $out                                                                                                                                                       
fi                                                                                                                                                                       
echo "DTI preproccessing begins ..."                                                                                                                                     
echo [`date`]                                                                                                                                                            
echo                                                                                                                                                                     
echo "Preparing topup"                                                                                                                                                   
${FSLDIR}/bin/imcp ${p2a} ${td}/                                                                                                                                         
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1                                                                                                                                
${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${td}/${b_in}_p2a_b0
echo "Running topup"            
${FSLDIR}/bin/topup --imain=${td}/a2p_p2a_b0 --datain=${out}/../acqparams.txt --out=${td}/topup_results --iout=${td}/hifi_b0                                             
${FSLDIR}/bin/fslmaths ${td}/hifi_b0 -Tmean ${td}/hifi_b0                                                                                                                
${FSLDIR}/bin/bet ${td}/hifi_b0 ${td}/hifi_b0_brain -m                                                                                                                   
#${FSLDIR}/bin/applytopup --imain=${in} --inindex=1 --datain=${out}/../acqparams.txt --topup=${td}/topup_results --out=${td}/${b_in}_hifi --method=jac                    
echo "Copying files for ${b_in} to ${td}/"                                                                                                                               
echo                                                                                                                                                                     
cp ${in}.bval ${td}/bvals
cp ${in}.bvec ${td}/bvecs 
echo "Doing correction on ${td}/${b_in}_hifi"                                                                                                                            
echo [`date`]
echo                                                                                                                                                                     
#${FSLDIR}/bin/eddy_correct ${td}/${b_in}_hifi ${td}/data 0
${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=${out}/../dti_index.txt --bvecs=${td}/bvecs --bvals=${td}/bvals --topup=${td}/topup_results --out=${td}/data 
echo "Running dtifit on ${td}/data"                                                                                                                                      
echo [`date`]
echo                                                                                                                                                                     
${FSLDIR}/bin/dtifit --data=${td}/data --out=${td}/dti --mask=${td}/hifi_b0_brain --bvecs=${td}/bvecs --bvals=${td}/bvals                                                
echo "I will copy all output files to ${out}/${b_in}_XXXXX"                                                                                                              
echo                                                                                                                                                                     
for x in ${td}/dti*; do ${FSLDIR}/bin/imcp ${x} ${out}/${b_in}_$(basename $x); done;                                                                                     
${FSLDIR}/bin/imcp ${td}/hifi_b0_brain ${out}/${b_in}_dti_brain_mask                                                                                                     
${FSLDIR}/bin/imcp ${td}/data ${out}/${b_in}_dti_data                                                                                                                    
echo [`date`]
exit 0

Ahora viene la parte buena. Hay que poner esto en un atlas o algo asi para sacar la info

Corregistro

Esto de aqui es el concept proof. Funciona pero he de probar un metodo mas exacto

En principio, si logramos registrar la plantilla ${FSLDIR}/data/standard/FMRIB58_FA_1mm a nuestro DTI y guardar la matriz de transformacion, podriamos aplicar esta matriz a algun atlas que esté en el espacio MNI.

$ ANTS 3 -m CC[working/mop0001_dti_FA.nii.gz, ${FSLDIR}/data/standard/FMRIB58_FA_1mm, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o working/mop0001_dti_ants_elast_b0_fmrib.nii.gz
$ WarpImageMultiTransform 3 ${FSLDIR}/data/standard/FMRIB58_FA_1mm working/.tmp_mop0001/FMRIB58_FA_t.nii.gz -R working/mop0001_dti_FA.nii.gz working/mop0001_dti_ants_elast_b0_fmribWarp.nii.gz working/mop0001_dti_ants_elast_b0_fmribAffine.txt

Pos no muy bien,

Vamos a intentar algo distinto. Si uso el T1 puedo obtener una tranformacion mas fiable para el espacio MNI, aunque añadiendo un paso intermedio. Vamos a darle mas info al algoritmo de registro.

Como diría Jack, el destripador, Vayamos por partes.

¿Podemos llevar satisfactoriamente el template de MNI a subject space?

$ fslreorient2std mri/mop0001s0005 working/mop0001_t1_reoriented
$ ANTS 3 -m CC[working/mop0001_t1_reoriented.nii.gz, ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1.nii.gz
$ WarpImageMultiTransform 3 ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz working/.tmp_mop0001/MNI152_warped_tmp.nii.gz -R working/mop0001_t1_reoriented.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Affine.txt

casi perfecto. Veamos que pasa con la plantilla de FA.

$ WarpImageMultiTransform 3 ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz working/.tmp_mop0001/FMRIB58_FA_warped_tmp.nii.gz -R working/mop0001_t1_reoriented.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Affine.txt

No está mal.

¿Se puede pasar esta plantilla a custom space (DTI)?

$ ANTS 3 -m CC[working/.tmp_mop0001/hifi_b0.nii.gz, working/mop0001_t1_reoriented.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0.nii.gz
$ WarpImageMultiTransform 3 working/.tmp_mop0001/MNI152_warped_tmp.nii.gz working/.tmp_mop0001/MNI152_custom_tmp.nii.gz -R working/.tmp_mop0001/hifi_b0.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Affine.txt

Buena pinta. Y para la FA,

$ WarpImageMultiTransform 3 working/.tmp_mop0001/FMRIB58_FA_warped_tmp.nii.gz working/.tmp_mop0001/FMRIB58_FA_custom_tmp.nii.gz -R working/.tmp_mop0001/hifi_b0.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Affine.txt

que no esta mal. Veamos ahora las dos FAs juntas.

que se ve bastante mejor que antes.

Resumiendo, lo que hemos de hacer es,

$ fslreorient2std mri/mop0001s0005 working/mop0001_t1_reoriented
$ ANTS 3 -m CC[working/mop0001_t1_reoriented.nii.gz, ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1.nii.gz
$ ANTS 3 -m CC[working/.tmp_mop0001/hifi_b0.nii.gz, working/mop0001_t1_reoriented.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0.nii.gz
$ WarpImageMultiTransform 3 ${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz working/.tmp_mop0001/JHU_labels_warped_tmp.nii.gz -R working/mop0001_t1_reoriented.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_mni_t1Affine.txt
$ WarpImageMultiTransform 3 working/.tmp_mop0001/JHU_labels_warped_tmp.nii.gz working/.tmp_mop0001/JHU_labels_custom_tmp.nii.gz -R working/.tmp_mop0001/hifi_b0.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Warp.nii.gz working/.tmp_mop0001/mop0001_dti_ants_elast_t1_b0Affine.txt

y tendremos el atlas JHU-labels registrado sobre nuestro archivo de FA.

Nota: Esto es solo un concept proof. El procedimiento sirve para pasar cualquier atlas del espacio MNI al DTI del sujeto. Lo he hecho con el JHU porque de aquí es muy sencillo extraer los valores.

Update: He cambiado el corregistro a la nueva version de ANTs. Sobre como usar ANTs. El paso de T1 a B0 se hace ahora usando epireg. Ver explicacion aqui y las notas en el procedimiento de FACEHBI.

echo "Calculating transformation from MNI to T1"
echo [`date`]
echo
antsRegistrationSyN.sh -d 3 -f ${out}/${b_in}_t1_reoriented.nii.gz -m ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -o ${td}/${b_in}_mni_t1_warp
antsApplyTransforms -d 3 -i ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_mni_t1_warped.nii.gz
 
echo "Calculating transformation from T1 to B0"                                                       
echo [`date`]
echo
${FSLDIR}/bin/epi_reg --epi=${out}/${b_in}_dti_data --t1=${out}/${b_in}_t1_reoriented --t1brain=${out}/${b_in}_brain --out=${td}/${b_in}_tmp_diff2std
${FSLDIR}/bin/convert_xfm -omat ${td}/${b_in}_tmp_std2diff.mat -inverse ${td}/${b_in}_tmp_diff2std.mat
${FSLDIR}/bin/flirt -in ${out}/${b_in}_t1_reoriented.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_t1_b0 -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm
${FSLDIR}/bin/flirt -in ${td}/${b_in}_mni_t1_warped.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_mni_to_b0 -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm

Esto es mucho mas lento pero los resultados deberian ser mucho mejores.

El script nuevo es mucho mas complicado y agrupa los procesos de preprocesado y corregistro.

dti_proc_epi.sh

#!/bin/sh 
 
Usage() {  
    echo "" 
    echo "Usage: dti_proc.sh  <study> <file_name> <T1> <name> <output_dir>"  
    echo ""  
    echo "You must have FSL installed in order to run this script" 
    echo "" 
    exit 1 
} 
 
[ "$5" = "" ] && Usage 
debug=1
study=$1
shift  
in=`${FSLDIR}/bin/remove_ext $1` 
shift  
t1=`${FSLDIR}/bin/remove_ext $1` 
shift  
b_in=$1  
shift  
out=`${FSLDIR}/bin/remove_ext $1` 
shift  
td=${out}'/.tmp_'${b_in}  
p2a=$(dirname ${in}.nii.gz)'/'${b_in}'_p2a_b0' 
if [ ! -d "$td" ]; then
        mkdir $td
fi
if [ ! -d "$out" ]; then 
        mkdir $out
fi
 
echo "DTI preproccessing begins on ${b_in} ..." 
echo [`date`] 
echo 
 
echo "Preparing topup"
${FSLDIR}/bin/imcp ${p2a} ${td}/                                                                   
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1                                                             
${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${td}/${b_in}_p2a_b0
 
echo "Running topup"
${FSLDIR}/bin/topup --imain=${td}/a2p_p2a_b0 --datain=${out}/../acqparams.txt --out=${td}/topup_results --iout=${td}/hifi_b0
${FSLDIR}/bin/fslmaths ${td}/hifi_b0 -Tmean ${td}/hifi_b0                                             
${FSLDIR}/bin/bet ${td}/hifi_b0 ${td}/hifi_b0_brain -m
 
echo "Copying files for ${b_in} to ${td}/"  
echo  
cp ${in}.bval ${td}/bvals
cp ${in}.bvec ${td}/bvecs 
 
echo "Doing correction on ${td}/${b_in}_hifi"  
echo [`date`]
echo  
#${FSLDIR}/bin/eddy_correct ${td}/${b_in}_hifi ${td}/data 0
${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=${out}/../dti_index.txt --bvecs=${td}/bvecs --bvals=${td}/bvals --topup=${td}/topup_results --out=${td}/data 
 
echo "Running dtifit on ${td}/data" 
echo [`date`]
echo 
${FSLDIR}/bin/dtifit --data=${td}/data --out=${td}/dti --mask=${td}/hifi_b0_brain --bvecs=${td}/bvecs --bvals=${td}/bvals                                                
 
echo "I will copy all output files to ${out}/${b_in}_XXXXX" 
echo  
for x in ${td}/dti*; do ${FSLDIR}/bin/imcp ${x} ${out}/${b_in}_$(basename $x); done; 
${FSLDIR}/bin/imcp ${td}/hifi_b0_brain ${out}/${b_in}_dti_brain_mask  
${FSLDIR}/bin/imcp ${td}/data ${out}/${b_in}_dti_data  
 
echo "I need the T1 image"  
${FSLDIR}/bin/fslreorient2std ${t1} ${out}/${b_in}_t1_reoriented  
 
echo "I get the extracted brain"
${FREESURFER_HOME}/bin/mri_vol2vol --mov ${SUBJECTS_DIR}/${study}_${b_in}/mri/brain.mgz --targ ${SUBJECTS_DIR}/${study}_${b_in}/mri/rawavg.mgz --regheader --o ${td}/${b_in}_tmp_brain_in_rawavg.mgz
${FREESURFER_HOME}/bin/mri_convert --in_type mgz --out_type nii ${td}/${b_in}_tmp_brain_in_rawavg.mgz ${td}/${b_in}_tmp_brain.nii.gz
${FSLDIR}/bin/fslreorient2std ${td}/${b_in}_tmp_brain ${out}/${b_in}_brain
 
echo "Calculating transformation from MNI to T1"
echo [`date`]
echo
antsRegistrationSyN.sh -d 3 -f ${out}/${b_in}_t1_reoriented.nii.gz -m ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -o ${td}/${b_in}_mni_t1_warp
antsApplyTransforms -d 3 -i ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_mni_t1_warped.nii.gz
 
echo "Calculating transformation from T1 to B0"                                                       
echo [`date`]
echo
${FSLDIR}/bin/epi_reg --epi=${out}/${b_in}_dti_data --t1=${out}/${b_in}_t1_reoriented --t1brain=${out}/${b_in}_brain --out=${td}/${b_in}_tmp_diff2std
${FSLDIR}/bin/convert_xfm -omat ${td}/${b_in}_tmp_std2diff.mat -inverse ${td}/${b_in}_tmp_diff2std.mat
${FSLDIR}/bin/flirt -in ${out}/${b_in}_t1_reoriented.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_t1_b0 -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm
${FSLDIR}/bin/flirt -in ${td}/${b_in}_mni_t1_warped.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_mni_to_b0 -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm
 
echo "Applying transformations to JHU Atlas (ICBM-labels)"
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_labels_tmp.nii.gz
${FSLDIR}/bin/flirt -in ${td}/${b_in}_JHU_labels_tmp.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_JHU_labels.nii.gz -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm
 
echo "Applying transformations to JHU Atlas (ICBM-tracts)"
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_tracts_tmp.nii.gz
${FSLDIR}/bin/flirt -in ${td}/${b_in}_JHU_tracts_tmp.nii.gz -ref ${td}/hifi_b0.nii.gz -out ${out}/${b_in}_JHU_tracts.nii.gz -init ${td}/${b_in}_tmp_std2diff.mat -applyxfm
 
${FSLDIR}/bin/imcp ${td}/hifi_b0 ${out}/${b_in}_hifi_b0
echo "Game over :-P"
echo [`date`]
exit 0

Procesamiento en paralelo

Pues ahora que tengo resuelto el problema lo que necesito es automatizar estas tareas de manera desatendida. Los procesos de preprocesamiento y corregistro los he puesto en el mismo script ya que utilizan la misma estructura y de cualquier manera deben ejecutarse secuencialmente.

Pero hay un problema!: Cuando no tenemos ningun B0 P»A no se puede aplicar topup y hay que vivir sin las correciones. Así que voy a preparar un script para esta situacion.

script sin topup

dti_proc_nocorr.sh

#!/bin/sh                                                                                                                                                                
 
Usage() {                                                                                                                                                                
    echo ""                                                                                                                                                              
    echo "Usage: dti_proc_nocorr.sh  <name> <T1> <input> <output_dir>"                                                                                                          
    echo ""                                                                                                                                                              
    echo "You must have FSL installed in order to run this script"                                                                                                       
    echo ""                                                                                                                                                              
    exit 1                                                                                                                                                               
}                                                                                                                                                                        
 
[ "$4" = "" ] && Usage                                                                                                                                                   
debug=1                                                                                                                                                                  
in=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                         
shift                                                                                                                                                                    
t1=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                         
shift                                                                                                                                                                    
b_in=$1                                                                                                                                                                  
shift                                                                                                                                                                    
out=`${FSLDIR}/bin/remove_ext $1`                                                                                                                                        
shift                                                                                                                                                                    
td=${out}'/.tmp_'${b_in}                                                                                                                                                 
#p2a=$(dirname ${in}.nii.gz)'/'${b_in}'_p2a_b0'                                                                                                                           
if [ ! -d "$td" ]; then                                                                                                                                                  
        mkdir $td                                                                                                                                                        
fi                                                                                                                                                                       
if [ ! -d "$out" ]; then                                                                                                                                                 
        mkdir $out                                                                                                                                                       
fi                                                                                                                                                                       
echo "DTI preproccessing begins ..."                                                                                                                                     
echo [`date`]                                                                                                                                                            
echo                                                                                                                                                                     
#echo "Preparing topup"                                                                                                                                                   
#${FSLDIR}/bin/imcp ${p2a} ${td}/                                                                                                                                         
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1                                                                                                                                
#${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${td}/${b_in}_p2a_b0
#echo "Running topup"            
#${FSLDIR}/bin/topup --imain=${td}/a2p_p2a_b0 --datain=${out}/../acqparams.txt --out=${td}/topup_results --iout=${td}/hifi_b0                                             
${FSLDIR}/bin/fslmaths ${td}/a2p_b0 -Tmean ${td}/hifi_b0                                                                                                                
${FSLDIR}/bin/bet ${td}/hifi_b0 ${td}/hifi_b0_brain -m                                                                                                                   
#${FSLDIR}/bin/applytopup --imain=${in} --inindex=1 --datain=${out}/../acqparams.txt --topup=${td}/topup_results --out=${td}/${b_in}_hifi --method=jac                    
echo "Copying files for ${b_in} to ${td}/"                                                                                                                               
echo                                                                                                                                                                     
#cp ${in}.bval ${td}/bvals
#cp ${in}.bvec ${td}/bvecs
echo "Doing correction on ${td}/${b_in}_hifi"                                                                                                                            
echo [`date`]
echo                                                                                                                                                                     
#${FSLDIR}/bin/eddy_correct ${td}/${b_in}_hifi ${td}/data 0                                           
${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=dti_index.txt --bvecs=${in}.bvec --bvals=${in}.bval --out=${td}/data
echo "Running dtifit on ${td}/data"                                                                                                                                      
echo [`date`]
echo                                                                                                                                                                     
${FSLDIR}/bin/dtifit --data=${td}/data --out=${td}/dti --mask=${td}/hifi_b0_brain --bvecs=${in}.bvec  --bvals=${in}.bval
echo "I will copy all output files to ${out}/${b_in}_XXXXX"                                                                                                              
echo                                                                                                                                                                     
for x in ${td}/dti*; do ${FSLDIR}/bin/imcp ${x} ${out}/${b_in}_$(basename $x); done;                                                                                     
${FSLDIR}/bin/imcp ${td}/hifi_b0_brain ${out}/${b_in}_dti_brain_mask                                                                                                     
${FSLDIR}/bin/imcp ${td}/data ${out}/${b_in}_dti_data                                                                                                                    
echo "I need the T1 image"                                                                                                                                               
${FSLDIR}/bin/fslreorient2std ${t1} ${out}/${b_in}_t1_reoriented                                                                                                             
echo "Calculating transformation from MNI to T1"                                                                                                                         
echo [`date`]
echo
ANTS 3 -m CC[${out}/${b_in}_t1_reoriented.nii.gz, ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o ${td}/${b_in}_dti_ants_elast_mni_t1.nii.gz                                                                                                                                                       
echo "Calculating transformation from T1 to B0"                                                                                                                          
echo [`date`]
echo
ANTS 3 -m CC[${td}/hifi_b0.nii.gz, ${out}/${b_in}_t1_reoriented.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o ${td}/${b_in}_dti_ants_elast_t1_b0.nii.gz       
echo "Applying transformations to JHU Atlas"                                                                                                                             
WarpImageMultiTransform 3 ${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz ${td}/JHU_labels_warped_tmp.nii.gz -R ${out}/${b_in}_t1_reoriented.nii.gz ${td}/${b_in}_dti_ants_elast_mni_t1Warp.nii.gz ${td}/${b_in}_dti_ants_elast_mni_t1Affine.txt                                                                                            
WarpImageMultiTransform 3 ${td}/JHU_labels_warped_tmp.nii.gz ${td}/JHU_labels_custom_tmp.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${b_in}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${b_in}_dti_ants_elast_t1_b0Affine.txt                                                                                                                                
echo "Game over :-P"
echo [`date`]
exit 0

Lo siguiente es reorganizar radicalmente el procesamiento. Hasta la versión anterior se ejecutaba todo en la misma maquina, haciendo forks. Ahora vamos a mover el procesamiento al schedule manager del cluster, por lo que debemos crear un archivo de configuracion para sbatch y ejecutar la tareas. Pero esto ya lo hemos hecho antes, asi que no deberia haber problemas.

Nota: En un primer paso hemos lanzado todos los trabajos en la misma tarea pero se ha visto que las GPUno son capaces de controlar esto, asi que ahora cada tarea se lanza independientemente y la politica de recursos de slurm los maneja. Es mas lento pero no fallan.

A ver si asi funciona OK

dti_reg.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
my $debug=1;
my $corr=1;
my $legacy=0;
my $old=0;
my $chop=0;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-nocorr/) {$corr=0;} 
    if (/^-legacy/) {$legacy=1;}
    if (/^-old/) {$old=1;}
    if (/^-chop/) {$old=1; $chop=1;}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
#my $pipe_dir = "/nas/software/neuro.dev/bin/";
# Redirect ouput to logfile (do it only when everything is fine)
my $logfile = "$std{'DATA'}/.debug_dti_register.log";
open STDOUT, ">$logfile" or die "Can't redirect stdout";
open STDERR, ">&STDOUT" or die "Can't dup stdout";
$debug ? open DBG, ">$logfile" :0;
 
#Run this script ONLY on "Detritus"
#or change the paths to the correct ones
 
my $w_dir=$std{'WORKING'};
my $dti_dir=$std{'DTI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
 
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
 
my %plist = populate($db);
my @ok_dtis;
 
print "Collecting needed files\n";
 
my %dtis;
if ($cfile){
        my %cuts = get_pair($data_dir."/".$cfile);
        foreach my $cut (keys %cuts){
                if(grep {/$cut/} %plist){
                        $dtis{$cut} = $plist{$cut};
                }
        }
}else{
        %dtis = %plist;
}
 
my $count = 0;
my @jobs_list;
foreach my $subject (sort keys %dtis){
        if($subject){
                my @names = find(file => 'name' => "$dtis{$subject}$subject*.nii.gz", in => $dti_dir);
                my @t1_names = find(file => 'name' => "$dtis{$subject}$subject*.nii.gz", in => $mri_dir);
                if(@names && @t1_names){
                        my $order;
                        unless($old){
                                if($corr){
                                        if($legacy){
                                                $order = $pipe_dir."/bin/dti_proc_legacy.sh ".$names[0]." ".$t1_names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                        }else{
                                                $order = $pipe_dir."/bin/dti_proc_epi.sh ".$study." ".$names[0]." ".$t1_names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                        }
                                }else{
                                        if($legacy){
                                                $order = $pipe_dir."/bin/dti_proc_nocorr_legacy.sh ".$names[0]." ".$t1_names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                        }else{
                                                $order = $pipe_dir."/bin/dti_proc_nocorr.sh ".$names[0]." ".$t1_names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                        }
                                }
                        }else{
                                if($chop){
                                        $order = $pipe_dir."/bin/dti_proc_x.sh ".$names[0]." ".$t1_names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                }else{
                                        $order = $pipe_dir."/bin/dti_proc_deprecated.sh ".$names[0]." ".$dtis{$subject}.$subject." ".$w_dir;
                                }
                        }
                        my $orderfile = $outdir.'/'.$subject.'dti_orders.sh';
                        open ORD, ">$orderfile";
                        print ORD '#!/bin/bash'."\n";
                        print ORD '#SBATCH -J dti_reg-'.$study."\n";
                        print ORD '#SBATCH --time=15:0:0'."\n"; #si no ha terminado en X horas matalo
                        print ORD '#SBATCH -c 4'."\n";
                        print ORD '#SBATCH --mail-type=FAIL,STAGE_OUT'."\n"; #no quieres que te mande email de todo
                        print ORD '#SBATCH -o '.$outdir.'/dti_reg-slurm-%j'."\n";
                        unless($legacy or $old){ print ORD '#SBATCH --gres=gpu:1'."\n"; }
                        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
                        unless($legacy or $old){ print ORD '#SBATCH -p cuda'."\n"; }
                        print ORD "srun $order\n";
                        close ORD;
                        my $jobid = `sbatch $orderfile`;
                        $jobid = ( split ' ', $jobid )[ -1 ];
                        push @jobs_list, $jobid;
                        $debug ? print DBG "$order\n" :0;
                }
        }
}
$debug ? close DBG:0;   
my $sjobs = join(',',@jobs_list);
my $orderfile = $outdir.'/dti_reg_end.sh';
open ORD, ">$orderfile";
print ORD '#!/bin/bash'."\n";
print ORD '#SBATCH -J dti_reg-'.$study."\n";
print ORD '#SBATCH --mail-type=END'."\n"; #email cuando termine
print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
print ORD '#SBATCH -o '.$outdir.'/dti_reg_end-%j'."\n";
print ORD ":\n";
close ORD;
my $order = 'sbatch --depend=afterany:'.$sjobs.' '.$orderfile;
exec($order);

Backward compatibility

¿Que hago si no tengo el P»A? Pues no se podría utilizar topup, pero si el nuevo eddy. A la hora de hacer el script hemos de incluir el switch -nocorr para que se salte el topup.

¿Y si eddy me falla? Pues usamos la version antigua. Para esto se incluye el switch -legacy como argumento del script.

Ejemplo: Para analizar los DTI de FACEHBI habría que utlizar

$ dti_reg.pl -nocorr -legacy facehbi

Configurar PATHs

Para que todo funcione OK, los bricks tienen que poder acceder a FSL y al pipeline. FSL da muchos problemas dada la diferencia de plataformas entre los brick y detritus por lo que no he tenido más remedio que instalarlo localmente en cada nodo . El pipeline lo he puesto en el NAS pero, para que todo funcione correctamente, en el .bash_profile debe existir la línea,

export PIPEDIR=/nas/software/neuro.dev

de nada.

Haciendo máscaras

Una vez que hemos puesto un atlas sobre nuestra imagen de FA, ya podemos hacer una mascara para cada etiqueta del atlas y seleccionar el valor medio de FA en esa región. Ejemplo: El atlas JHU los he desglosado en el archivo $PIPEDIR/lib/jhu.labels.list

labels_0;Unclassified
labels_1;Middle_cerebellar_peduncle
labels_2;Pontine_crossing_tract_(a_part_of_MCP)
labels_3;Genu_of_corpus_callosum
labels_4;Body_of_corpus_callosum
...
...
...
labels_45;Uncinate_fasciculus_R
labels_46;Uncinate_fasciculus_L
labels_47;Tapetum_R
labels_48;Tapetum_L

Así que lo primero será hacer una máscara para cada una de las ROIs del archivo JHU_labels_custom_tmp.nii.gz correspondiente a cada sujeto. Por ejemplo, la orden,

$ fslmaths working/.tmp_mop0001/JHU_labels_custom_tmp.nii.gz -uthr 1 -thr 1 -div 1 working/.tmp_mop0001/JHU_labels_1_tmp.nii.gz

hará una máscara binaria del Middle_cerebellar_peduncle ~~(lo que demonios sea eso)~~ del sujeto 0001. Luego sólo habrá que aplicar esta máscara sobre el FA del sujeto para calcular el valor medio y la SD de la FA en esa ROI

$ fslstats working/mop0001_dti_FA -k working/.tmp_mop0001/JHU_labels_1_tmp -M -S
0.446217 0.177036

Se pone en una tabla y ya está!

Bueno, hay que hacerlo para cada sujeto y para cada una de las ROI.

Voy a echarle un poco de ganas

dti_metrics.pl

#!/usr/bin/perl                                                                                       
 
use strict; use warnings;                                                                             
 
use File::Find::Rule;                                                                                 
use NEURO qw(print_help get_pair load_study achtung shit_done);                                       
use Data::Dump qw(dump);                                                                              
use File::Basename qw(basename);                                                                      
 
my $attach = 1;                                                                                       
my $with_sd = 0;                                                                                      
my $atlas = 1;                                                                                        
@ARGV = ("-h") unless @ARGV;                                                                          
while (@ARGV and $ARGV[0] =~ /^-/) {                                                                  
    $_ = shift;                                                                                       
    last if /^--$/;                                                                                   
    if (/^-l/) { $attach = 0;}                                                                                                                                                
    if (/^-sd/) { $with_sd = 1;}                                                                      
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_metrics.hlp'; exit;}                            
}
my $study = shift;                                                                                    
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_metrics.hlp'; exit;}                           
my %std = load_study($study);                                                                         
 
my $w_dir=$std{'WORKING'};                                                                            
my $data_dir=$std{'DATA'};                                                                            
my $roi_list = $ENV{'PIPEDIR'}.'/lib/jhu.labels.list';                                                                                                                                      
my $roi_dir = $ENV{'PIPEDIR'}.'/lib/jhu/';                                                            
my $jhu = "JHU_labels_custom_tmp";                                                                    
my %masks = get_pair $roi_list;                                                                       
my @dti_fas = find(file => 'name' => "*_dti_FA.nii.gz", in => $w_dir);                                                         
my %csv;                                                                                              
 
my $ofile = $data_dir."/".$study."_dti.csv";                                                          
open OF, ">$ofile";                                                                                   
 
print OF "Subject"; 
foreach my $rmask (sort keys %masks){                                                                 
        if($with_sd){                                                                                 
                print OF ";$masks{$rmask}","_FA_Mean;","$masks{$rmask}","_FA_STD";                    
        }else{                                                                                        
                print OF ";$masks{$rmask}","_FA";                                                     
        }                                                                                             
}                                                                                                     
print OF "\n";                                                                                        
 
foreach my $dti (sort @dti_fas){                                                                      
        (my $a,my $subject) = $dti =~ /(\w{1,3})(\d{1,6})_dti_FA/;                                                                               
        if($subject){                                                                                 
                foreach my $rmask (sort keys %masks){                                                 
                        my ($m_index) = $rmask =~ /.*_(\d+)$/;                                        
                        if($m_index){                                                                 
                                my $order = "fslmaths ".$w_dir."/.tmp_".$a.$subject."/".$jhu." -uthr ".$m_index." -thr ".$m_index." -div ".$m_index." ".$w_dir."/.tmp_".$a.$subject."/JHU_labels_".$rmask."_tmp";                                                                                                 
                        print "$order\n";                                                             
                        system($order);                                                               
                        } else {
                        my $order = "fslmaths ".$w_dir."/.tmp_".$a.$subject."/".$jhu." -uthr ".$m_index." -thr ".$m_index." ".$w_dir."/.tmp_".$a.$subject."/JHU_labels_".$rmask."_tmp";                     
                        print "$order\n";                                                                                                                                                                   
                        system($order);                                                               
                        }                                                                             
                }                                                                                     
                print OF "$subject";                                                                  
                foreach my $rmask (sort keys %masks){                                                 
                        my $order = "fslstats ".$dti." -k ".$w_dir."/.tmp_".$a.$subject."/JHU_labels_".$rmask."_tmp -M -S";                                                                                 
                        print "$order\n";                                                             
                        (my $mean, my $std) = map{/(\d+\.\d*)\s*(\d+\.\d*)/} qx/$order/;              
                        if($with_sd){                                                                 
                                print OF ";$mean",";$std";                                            
                        }else{                                                                        
                                print OF ";$mean";                                                    
                        }
                }                                                                                     
                print OF "\n";                                                                        
        }                                                                                             
}                                                                                                     
 
close OF;                                                                                             
 
my $zfile = $ofile.'.gz';                                                                             
system("gzip -c $ofile > $zfile");                                                                    
 
if ($attach){                                                                                         
        shit_done basename($ENV{_}), $study, $zfile;                                                  
}else{                                                                                                
        achtung basename($ENV{_}), $ofile, $study;                                                    
}

y tachán!

Nota: Esto se puede paralelizar pero no vale la pena. Demorará un ratejo cuando sean 200 o 500 pero en general es rápido. Demora más en mostrar en pantalla lo que hace que en hacerlo (opción: se pueden quitar los print una vez que este terminado).

Tractografía (How to)

bedpostx

El primer paso para hacer la tractografía es correr bedpostx sobre los resultados. Primero organizaremos el directorio de trabajo.

$ mkdir working/.tmp_mop0002/bedpostx
$ imcp working/mop0002_dti_data working/.tmp_mop0002/bedpostx/data
$ imcp working/mop0002_dti_brain_mask working/.tmp_mop0002/bedpostx/nodif_brain_mask
$ cp working/.tmp_mop0002/bvecs working/.tmp_mop0002/bedpostx/bvecs
$ cp working/.tmp_mop0002/bvals working/.tmp_mop0002/bedpostx/bvals
$ bedpostx working/.tmp_mop0002/bedpostx

y a esperar que son 15 horas aprox. El resultado es un nuevo directorio con todos los archivos,

$ ls working/.tmp_mop0002/bedpostx.bedpostX/
bvals                        dyads3.nii.gz                mean_ph1samples.nii.gz   merged_ph1samples.nii.gz
bvecs                        dyads3_thr0.05_modf3.nii.gz  mean_ph2samples.nii.gz   merged_ph2samples.nii.gz
commands.txt                 dyads3_thr0.05.nii.gz        mean_ph3samples.nii.gz   merged_ph3samples.nii.gz
dyads1_dispersion.nii.gz     logs                         mean_S0samples.nii.gz    merged_th1samples.nii.gz
dyads1.nii.gz                mean_dsamples.nii.gz         mean_th1samples.nii.gz   merged_th2samples.nii.gz
dyads2_dispersion.nii.gz     mean_d_stdsamples.nii.gz     mean_th2samples.nii.gz   merged_th3samples.nii.gz
dyads2.nii.gz                mean_f1samples.nii.gz        mean_th3samples.nii.gz   monitor
dyads2_thr0.05_modf2.nii.gz  mean_f2samples.nii.gz        merged_f1samples.nii.gz  nodif_brain_mask.nii.gz
dyads2_thr0.05.nii.gz        mean_f3samples.nii.gz        merged_f2samples.nii.gz  xfms
dyads3_dispersion.nii.gz     mean_fsumsamples.nii.gz      merged

Estos son los archivos para correr probtrackx.

probtrackx

Lo proximo que necesito es una mascara de alguna ROI que se use como seed para calcular las fibras. Si quiero usar una region ya definida me voy a traer la segmentacion de freesurfer primero,

Nota: Ya tenemos un script que convierte el aseg.mgz de freesurfer en el espacio del sujeto , así que vamos a usarlo.

Nota: Usaremos el script get_aparc.sh que convierte el aparc+aseg.mgz al espacio del sujeto. El get_aseg.sh convierte el aseg.mgz pero este archivo solo contiene las regiones subcorticales y el cortex como un todo mientras que en el aparc estan parceladas las ROIs del cortex.

$ get_aparc.sh mopead mop0002 /nas/data/mopead/working/
Number of labels: 0
Annot File:      (null)
Template Volume: /nas/data/subjects/mopead_mop0002/mri/rawavg.mgz
Outut Volume: /nas/data/mopead/working//mop0002_tmp_aseg.mgz
Registration File: (null)
Fill Threshold: 0
Label Vox Vol:  1
ProjType:       (null)
ProjTypeId:     0
ProjStart:      0
ProjStop:       0
ProjDelta:      0.1
Subject:  (null)
Hemi:     (null)
UseNewASeg2Vol:  0
DoLabelStatVol  0
LabelCodeOffset  0
setenv SUBJECTS_DIR /nas/data/subjects
$Id: mri_label2vol.c,v 1.46 2014/12/08 21:11:54 greve Exp $
Template RAS-to-Vox: --------
-0.83333  -0.00000  -0.00000   87.99999;
-0.00000  -0.00000  -0.94815   119.99999;
-0.00000   0.94815  -0.00000   127.99999;
-0.00000  -0.00000  -0.00000   1.00000;
Template Voxel Volume: 1.33485
nHits Thresh: 0
Computing registration based on header
RegMat: --------
-0.99932  -0.03498  -0.01203  -0.00003;
-0.01575   0.10824   0.99400   0.00005;
 0.03347  -0.99351   0.10872  -0.00003;
 0.00000   0.00000   0.00000   1.00000;
Label RAS-to-Vox: --------
 0.83276   0.02915   0.01002   88.00002;
-0.03173   0.94199  -0.10308   120.00002;
-0.01494   0.10263   0.94246   128.00005;
 0.00000   0.00000   0.00000   1.00000;
mri_label2vol done
mri_convert.bin --in_type mgz --out_type nii /nas/data/mopead/working//mop0002_tmp_aseg.mgz /nas/data/mopead/working//mop0002_tmp_aseg.nii.gz 
$Id: mri_convert.c,v 1.226 2016/02/26 16:15:24 mreuter Exp $
reading from /nas/data/mopead/working//mop0002_tmp_aseg.mgz...
TR=2200.00, TE=0.00, TI=0.00, flip angle=0.00
i_ras = (0.999316, 0.034981, 0.0120294)
j_ras = (-0.033469, 0.993509, -0.108721)
k_ras = (-0.0157545, 0.108244, 0.994)
writing to /nas/data/mopead/working//mop0002_tmp_aseg.nii.gz...
$ ls /nas/data/mopead/working//mop0002_aseg.nii.gz
/nas/data/mopead/working//mop0002_aseg.nii.gz

Ahora debería pasar este archivo al espacio DTI. Pero esto ya lo habiamos calculado antes para pasar la T1, así que debería funcionar,

$ WarpImageMultiTransform 3 /nas/data/mopead/working/mop0002_aseg.nii.gz /nas/data/mopead/working/.tmp_mop0002/mop0002_aseg_warped.nii.gz -R /nas/data/mopead/working/.tmp_mop0002/hifi_b0.nii.gz /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Warp.nii.gz /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Affine.txt
FIELD: /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Warp.nii.gz
AFFINE: /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Affine.txt
moving_image_filename: /nas/data/mopead/working/mop0002_aseg.nii.gz components 1
output_image_filename: /nas/data/mopead/working/.tmp_mop0002/mop0002_aseg_warped.nii.gz
reference_image_filename: /nas/data/mopead/working/.tmp_mop0002/hifi_b0.nii.gz
[0/2]: FIELD: /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Warp.nii.gz
[1/2]: AFFINE: /nas/data/mopead/working/.tmp_mop0002/mop0002_dti_ants_elast_t1_b0Affine.txt
User Linear interpolation 
  We check the syntax of your call .... 
 syntax probably ok. 
output origin: [-117.606, 109.699, -28.4168]
output size: [112, 112, 64]
output spacing: [2.00893, 2.00893, 2]
output direction: 0.999047 0.0436588 2.13526e-08
0.0436588 -0.999047 -9.77694e-07
2.13526e-08 -9.77694e-07 1

que nos alcanzara para hacer una mascara de alguna ROI aproximada. Digamos que queremos los hipocampos,

$ fslmaths working/.tmp_mop0002/mop0002_aseg_warped -uthr 17 -thr 17 -div 17 working/.tmp_mop0002/mop0002_tmp_mask1 
$ fslmaths working/.tmp_mop0002/mop0002_aseg_warped -uthr 53 -thr 53 -div 53 working/.tmp_mop0002/mop0002_tmp_mask2 
$ fslmaths working/.tmp_mop0002/mop0002_tmp_mask1 -add working/.tmp_mop0002/mop0002_tmp_mask2 working/.tmp_mop0002/mop0002_tmp_mask

Y ahora vamos a intentarlo,

$ probtrackx2 --opd --forcedir -s /nas/data/mopead/working/.tmp_mop0002/bedpostx.bedpostX/merged -m /nas/data/mopead/working/mop0002_dti_brain_mask -x /nas/data/mopead/working/.tmp_mop0002/mop0002_tmp_mask --dir='/nas/data/mopead/working/.tmp_mop0002/ptout'
Log directory is: /nas/data/mopead/working/.tmp_mop0002/ptout
Running in seedmask mode
load seeds
done.
Load bedpostx samples
1_1
1_2
1_3
2_1
2_2
2_3
3_1
3_2
3_3

nfibres  : 3
nsamples : 50

Done loading samples.
Volume seeds
volume 1

time spent tracking: 206 seconds

save results
finished

Algo ha hecho, pero todavia no se que,

$ ls working/.tmp_mop0002/ptout/
fdt_paths.nii.gz  probtrackx.log  waytotal

Un vistazo con FSLeyes,

Falta ponerlo bonito, escoger las ROI adecuadas para seeds, organizar un poco y limpiar los temporales, pero la idea general esta.

Tractografía (Allons-y!)

Después de varias pruebas he visto que bedpostx puede demorar mas de 24 horas en hacer su tarea. Afortunadamente, la nueva version de FSL viene con bedpostx_gpu. Como su nombre indica, utiliza la GPU para calcular y termina en unos 15 minutos.

Por partes, primero metemos todo en un script,

Por ahora sin calcular las mascaras

dti_bedtrack_cuda_nomasks.sh

#!/bin/sh                                                                                                                                                                
 
Usage() {                                                                                                                                                                
    echo ""                                                                                                                                                              
    echo "Usage: dti_bedtrack.sh  <project> <name> <working dir>"                                                                                                          
    echo ""                                                                                                                                                              
    echo "You must have FSL installed in order to run this script"                                                                                                       
    echo ""                                                                                                                                                              
    exit 1                                                                                                                                                               
}                                                                                                                                                                        
 
[ "$3" = "" ] && Usage
debug=1
prj=$1
shift
pollo=$1
shift
w_dir=$1
shift
td=${w_dir}'/.tmp_'${pollo}
bd=${td}'/bedpostx'
if [ ! -d "$bd" ]; then                                                                                                                                                  
        mkdir $bd                                                                                                                                                        
fi 
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_data ${bd}/data
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_brain_mask ${bd}/nodif_brain_mask
cp ${td}/bvecs ${bd}/bvecs
${FSLDIR}/bin/bedpostx_gpu ${bd}
${PIPEDIR}/bin/get_aseg.sh ${prj} ${pollo} ${w_dir}
WarpImageMultiTransform 3 ${w_dir}/${pollo}_aseg.nii.gz ${td}/${pollo}_aseg_warped.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt
########################################################
######## Falta calcular las mascaras para seeds ########
########################################################
probtrackx2 --opd --forcedir -s ${bd}.bedpostX/merged -m ${w_dir}/${pollo}_dti_brain_mask -x ${td}/${pollo}_tmp_mask --dir=\'${td}/probtrack_out\'

probtrackx necesita una mascara o grupo de mascaras para tener un punto de partida del cual buscar las fibras que conectan, ya sea entre ROIs o que entran y salen de una ROI especifica (se pueden hacer muchas cosas pero por ahora me quedo con lo simple).

Esto tiene un problema logistico

Lo mas sencillo sería construir una mascar predeterminada y hacerlo siempre igual, la meto dentro del proceso a paralelizar y ya esta. Pero entonces habria que cambiar el script cada vez que se quiere cambiar de mascara, o hacer un script para cada mascara. No sirve
También se podrian hacer las mascaras antes del proceso de paralelizacion, en el manager (dti_noseque.pl). Pero entonces habria que hacerlas en serie, o hacer dos procesos de paralelizacion con dependencias. Ineficiente, No sirve
Para poner la decision dentro de los procesos paralelizados, tendriamos que modificar el shell script para leer un archivo externo con una lista de las mascaras y hacer las mascaras on the fly. Puedo entonces hacer un archivo texto con la lista de mascaras que le paso a probtrackx y que las vaya cargando. Me gusta

Voy a usar la segmentacion de freesurfer para hacer las mascaras, asi que lo mas sencillo es utilizar propias etiquetas de freesurfer, o mejor, los valores que asigna freesurfer a la segmentacion y que son siempre los mismos (FreeSurferColorLUT.txt o ASegStatsLUT.txt).

Entonces, para correr el proceso habria que hacer una archivo masks.list en el directorio del projecto que contenga los codigos de las ROI que queremos introducir en probtrackx. Ejemplo, los hipocampos,

$ cat masks.list 
17
53

Debo entonces añadir un par de lineas al script y cambiar la llamada a probtrackx,

###########################################
for x in `cat $list`; do
        ${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${x} -thr ${x} -div ${x} ${td}/${pollo}_mask_${x};
        echo "${td}/${pollo}_mask_${x}.nii.gz" >> ${td}/${pollo}_mask.list;
done
###########################################
probtrackx2 --opd --forcedir -s ${bd}.bedpostX/merged -m ${w_dir}/${pollo}_dti_brain_mask -x ${td}/${pollo}_mask.list --dir=${td}/probtrack_out
rm ${td}/${pollo}_mask.list;

Asi es como queda entonces

dti_bedtrack_cuda.sh

#!/bin/sh                                                                                                                                                                
 
Usage() {                                                                                                                                                                
    echo ""                                                                                                                                                              
    echo "Usage: dti_bedtrack.sh  <project> <name> <working dir>"                                                                                                                                                                                                                                                                                 
    echo ""                                                                                                                                                              
    echo "You must have FSL installed in order to run this script"                                                                                                       
    echo ""                                                                                                                                                              
    exit 1                                                                                                                                                               
}                                                                                                                                                                        
 
[ "$3" = "" ] && Usage                                                                                                                                                   
debug=1                                                                                                                                                                  
prj=$1                                                                                                                                                                   
shift                                                                                                                                                                    
pollo=$1                                                                                                                                                                 
shift                                                                                                                                                                    
w_dir=$1                                                                                                                                                                 
shift                                                                                                                                                                    
td=${w_dir}'/.tmp_'${pollo}                                                                                                                                              
bd=${td}'/bedpostx'                                                                                                                                                      
list=${w_dir}/../masks.list                                                                                                                                              
if [ ! -d "$bd" ]; then                                                                                                                                                  
        mkdir $bd                                                                                                                                                        
fi                                                                                                                                                                       
echo "Copying files"                                                                                                                                                     
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_data ${bd}/data                                                                                                                 
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_brain_mask ${bd}/nodif_brain_mask                                                                                               
cp ${td}/bvecs ${bd}/bvecs                                                                                                                                               
cp ${td}/bvals ${bd}/bvals                                                                                                                                               
echo "Making bedpostx"                                                                                                                                                   
echo [`date`]                                                                                                                                                            ${FSLDIR}/bin/bedpostx_gpu ${bd}                                                                                                                                         
echo "So far, so good"                                                                                                                                                   
echo [`date`]                                                                                                                                                            
echo "Getting aseg"                                                                                                                                                      
${PIPEDIR}/bin/get_aseg.sh ${prj} ${pollo} ${w_dir}                                                                                                                      
WarpImageMultiTransform 3 ${w_dir}/${pollo}_aseg.nii.gz ${td}/${pollo}_aseg_warped.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt 
for x in `cat $list`; do                                                                                                                                                 
        ${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${x} -thr ${x} -div ${x} ${td}/${pollo}_mask_${x};                                                       
        echo "${td}/${pollo}_mask_${x}.nii.gz" >> ${td}/${pollo}_mask.list;                                                                                              
done                                                                                                                                                                     
echo "Doing probtrackx"                                                                                                                                                  
probtrackx2 --opd --forcedir -s ${bd}.bedpostX/merged -m ${w_dir}/${pollo}_dti_brain_mask -x ${td}/${pollo}_mask.list --dir=${td}/probtrack_out                          
rm ${td}/${pollo}_mask.list;                                                                                                                                             
echo "Done"                                                                                                                                                              
echo [`date`]

Un vistazo,

Organizamos la fiesta. El problema aquí es que cuando haya muchos sujetos que procesar al mismo tiempo no vamos a poder lanzarlos indiscriminadamente (uno por CPU) pues la tarjeta GPU no será capaz de procesar todo, las tareas fallarán y se perderá mucho tiempo revisando. Así que voy a lanzarlas por separado y dejar que Slurm se encargue de decirme si algo va mal.

Aqui el codigo, despues tendre que cambiarlo (ver abajo)

dti_track.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
my $debug=1;
my $corr=1;
my $legacy=0;
my $cuda=1;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-slow/) {$cuda=0;} 
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
 
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
 
# Redirect ouput to logfile (do it only when everything is fine)
my $logfile = "$std{'DATA'}/.debug_dti_tracks.log";
open STDOUT, ">$logfile" or die "Can't redirect stdout";
open STDERR, ">&STDOUT" or die "Can't dup stdout";
$debug ? open DBG, ">$logfile" :0;
 
#Run this script ONLY on "Detritus"
#or change the paths to the correct ones
 
my $w_dir=$std{'WORKING'};
my $dti_dir=$std{'DTI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
 
#open slurm file
#my $orderfile = "$std{'DATA'}/dti_track_orders.sh";
#my $conffile = "$std{'DATA'}/dti_track_orders.conf";
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
#my $jobname = "track_$study";
 
my %plist = populate($db);
my @ok_dtis;
 
print "Collecting needed files\n";
 
my %dtis;
 
if ($cfile){
	my %cuts = get_pair($data_dir."/".$cfile);
	foreach my $cut (keys %cuts){
		if(grep {/$cut/} %plist){
			$dtis{$cut} = $plist{$cut};
		}
	}
}else{
	%dtis = %plist;
}
#open CNF, ">$conffile";
#my $count = 0;
foreach my $subject (sort keys %dtis){
	#Compulsory: el DTI debe haber pasado el preproc
	my $dti_fa = $w_dir.'/'.$dtis{$subject}.$subject.'_dti_FA.nii.gz';
	print "$dti_fa\n";
	#Compulsory: Necesito tambien la segmentacion de freesurfer de este pollo
	my $aseg = $subj_dir.'/'.$study.'_'.$dtis{$subject}.$subject.'/mri/aseg.mgz';
	print "$aseg\n";
	if((-e $dti_fa) && (-e $aseg)){
		my $order;
		if($cuda){
			#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
			$order = $pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
		}else{
			#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
			$order = $pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
		}
		#$count++;
		#print CNF "$order\n";
		my $orderfile = $outdir.'/'.$subject.'dti_orders.sh';
                open ORD, ">$orderfile";
		print ORD '#!/bin/bash'."\n";
                print ORD '#SBATCH -J dti_track_'.$subject."\n";
                print ORD '#SBATCH --time=1:0:0'."\n"; #si no ha terminado en X horas matalo
                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
		if($cuda){ print ORD '#SBATCH -p cuda'."\n";}
		print ORD "srun $order\n";
                close ORD;
                my $pid = fork();
                if ($pid == 0){
                        exec("sbatch $orderfile");
                }
                elsif (!defined($pid)){
                        die "could not fork";
		}
		$debug ? print DBG "$order\n" :0;
	}
}
#close CNF;
 
#open ORD, ">$orderfile";
#print ORD '#!/bin/bash'."\n";
#print ORD '#SBATCH --time=0'."\n";
#print ORD '#SBATCH -J '.$track_study."\n";
#print ORD '#SBATCH -n '.$count."\n";
#print ORD '#SBATCH -o '.$outdir.'/dti_track-slurm-%j'."\n";
#print ORD '#SBATCH --mail-type=ALL'."\n";                                                             
#print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
#if($cuda){print ORD '#SBATCH -p cuda'."\n";}
#print ORD "srun --exclusive --multi-prog ".$conffile."\n";
#close ORD;
 
$debug ? close DBG:0;   
#exec("sbatch -o $w_dir/dti_track.out $orderfile");

Nota: Dejo las lineas comentadas para darme cuenta de lo que he cambiado

Testing

[osotolongo@detritus mopead]$ dti_track.pl mopead
[osotolongo@detritus mopead]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
              1046      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1047      cuda dti_trac osotolon PD       0:00      1 (Priority)
              1048      cuda dti_trac osotolon PD       0:00      1 (Priority)
              1049      cuda dti_trac osotolon PD       0:00      1 (Priority)
              1050      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1051      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1052      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1053      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1054      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1055      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1056      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1057      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1058      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1059      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1060      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1044      cuda dti_trac osotolon  R       0:03      1 brick01
              1045      cuda dti_trac osotolon  R       0:03      1 detritus

y al otro día ahi están,

$ find /nas/data/mopead/ -name "fdt_paths.nii.gz"
/nas/data/mopead/working/.tmp_mop0002/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0003/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0004/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0005/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0006/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0007/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0008/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0009/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0010/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0012/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0013/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0016/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0017/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0018/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0019/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0021/probtrack_out/fdt_paths.nii.gz
/nas/data/mopead/working/.tmp_mop0023/probtrack_out/fdt_paths.nii.gz

Nota: No estan todos, para los que falló el proceso de segmentacion de freesurfer y/o el dtifit no puedo calcularlos. En este caso toca revisar e intentar arreglar los problemas.

Nota: bedpostx solo se ejecutará si el directorio de output no existe. Esto es útil porque podemos relanzar probtrackx sin necesidad de recalcular todo pero en caso de querer ejecutar de nuevo bedpostx hay que limpiar los directorios,

$ for x in `find /nas/data/mopead/working/ -name "*bedpostX*" -type d`; do rm -rf $x; done

Tractografía (Otros atlas?)

Es posible implementar este proceso para cualquier otro atlas.No obstante, hay que saber de antemano cual es la organizacion del atlas y que se necesita. En principio, los atlas deberian venir en el espacio MNI o afín. Dado que hemos usado la segmentación de Freesurfer,solo hemos hecho una transformación pero estos atlas exigen dos, del espacio MNI al espacio de sujeto (T1) y de ahi al espacio custom (DTI).

Vamos a intentar implementar la tractografia para el Atlas de UofM-JUH. Este atlas tiene definidas distintas zonas (nodos). Voy a intentar usar estos nodos para calcular la tractografia con ellos. Asi que voy a borrar todo lo demas, organizar esto e intentar proceder de manera similar a como hice para la segmentacion de FS.

La lista de redes es,

$ cat uofm.list 
DMN Default_Mode_Network
ECN Executive_Control_Network
SN Salience_Network
AN Auditory_Network
BGN Basal_Ganglia_Network
HVN Higher_Visual_Network
LN Language_Network
PN Precuneus_Network
PVN Primary_Visual_Network
SMN SensoriMotor_Network
VSN VisuoSpatial_Network

Algunas redes se dividen en otras,

$ ls
AN  BGN  DMN  ECN  HVN  LN  PN  PVN  SMN  SN  uofm.list  VSN
$ ls LN
mri_LN1.nii  mri_LN2.nii  mri_LN3.nii  mri_LN4.nii  mri_LN5.nii  mri_LN6.nii  mri_LN7.nii
$ ls DMN/
dorsal  ventral
$ ls DMN/dorsal/
mri_D1.nii  mri_D3.nii  mri_D5.nii  mri_D7.nii  mri_D9.nii
mri_D2.nii  mri_D4.nii  mri_D6.nii  mri_D8.nii

Pero como primer paso voy a ignorarlo y usar todos los nodos. Despues vere como tratar con las subdivisiones. La cosa sería hacer algo así después del bedpostx,

###########################################                                                                                                                              
echo "Getting nodes and making masks"                                                                                                                                    
        node=$(basename $x)                                                                                                                                              
        WarpImageMultiTransform 3 ${x} ${td}/${node%.nii.gz}_warped.nii.gz -R ${w_dir}/${pollo}_t1_reoriented.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Warp.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Affine.txt;                                                                                                                              
        WarpImageMultiTransform 3 ${td}/${node%.nii.gz}_warped.nii.gz ${td}/${pollo}_${node}.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt;                                                                                                                               
        echo "${td}/${pollo}_${node}.gz" >> ${td}/${pollo}_mask.list;                                                                                            
done                                                                                                                                                                     
###########################################

Pero para ello necesito pasarle al script el directorio donde buscar los nodos. Esto no debería ser dificil pero vamos a probarlo antes.

Lo meto en otro script,

dti_bedtrack_uofm.sh

#!/bin/sh                                                                                                                                                                
 
Usage() {                                                                                                                                                                
    echo ""                                                                                                                                                              
    echo "Usage: dti_bedtrack.sh  <project> <name> <working dir> <network dir>"                                                                                          
    echo ""                                                                                                                                                              
    echo "You must have FSL installed in order to run this script"                                                                                                       
    echo ""                                                                                                                                                              
    exit 1                                                                                                                                                               
}                                                                                                                                                                        
 
[ "$4" = "" ] && Usage                                                                                                                                                   
debug=1                                                                                                                                                                  
prj=$1                                                                                                                                                                   
shift                                                                                                                                                                    
pollo=$1                                                                                                                                                                 
shift                                                                                                                                                                    
w_dir=$1                                                                                                                                                                 
shift                                                                                                                                                                    
net_dir=$1                                                                                                                                                               
shift                                                                                                                                                                    
td=${w_dir}'/.tmp_'${pollo}                                                                                                                                              
bd=${td}'/bedpostx'                                                                                                                                                      
list=${w_dir}/../masks.list                                                                                                                                              
if [ ! -d "$bd" ]; then                                                                                                                                                  
mkdir $bd                                                                                                                                                                
fi                                                                                                                                                                       
echo "Copying files"                                                                                                                                                     
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_data ${bd}/data                                                                                                                 
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_brain_mask ${bd}/nodif_brain_mask                                                                                               
cp ${td}/bvecs ${bd}/bvecs                                                                                                                                               
cp ${td}/bvals ${bd}/bvals                                                                                                                                               
echo "Making bedpostx"                                                                                                                                                   
echo [`date`]                                                                                                                                                            
${FSLDIR}/bin/bedpostx_gpu ${bd}                                                                                                                                         
echo "So far, so good"                                                                                                                                                   
echo [`date`]
###########################################                                                                                                                              
echo "Getting nodes and making masks"                                                                                                                                    
for x in `find ${net_dir} -name "*.nii"`; do                                                                                                                             
        node=$(basename $x)                                                                                                                                              
        WarpImageMultiTransform 3 ${x} ${td}/${node%.nii.gz}_warped.nii.gz -R ${w_dir}/${pollo}_t1_reoriented.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Warp.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Affine.txt;                                                                                                                              
        WarpImageMultiTransform 3 ${td}/${node%.nii.gz}_warped.nii.gz ${td}/${pollo}_${node}.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt;                                                                                                                             
        echo "${td}/${pollo}_${node}.gz" >> ${td}/${pollo}_mask.list;                                                                                                    
done                                                                                                                                                                     
###########################################                                                                                                                              
echo "Doing probtrackx"                                                                                                                                                  
probtrackx2 --opd --forcedir -s ${bd}.bedpostX/merged -m ${w_dir}/${pollo}_dti_brain_mask -x ${td}/${pollo}_mask.list --dir=${td}/probtrack_out                          
rm ${td}/${pollo}_mask.list;                                                                                                                                             
echo "Done"                                                                                                                                                              
echo [`date`]

y alla vamos,

$ /opt/neuro.dev/bin/dti_bedtrack_nodes.sh mopead mop0002 /nas/data/mopead/working/ /opt/neuro.dev/lib/uofm/LN/

y hasta funciona, ese es el Language Network,

Pues ahora esto hay que organizarlo de tal manera que sea posible incluir cualquier red y cualquier atlas si se organiza igual.

Tractografía (tojunto)

Ahora tengo que modificar el script que lanza los procesos para que ejecute este nuevo batch usando algun switch. O sea, por defecto se hara el probtrack con la segmentacion de freesurfer pero si paso un determinado switch le puedo indicar otras opciones.

Voy a hacerlo directamente con la GPU, no vale la pena perder tiempo,

dti_track.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
my $debug=1;
my $corr=1;
my $legacy=0;
my $cuda=1;
my $atlas=0; 
my $network;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-slow/) {$cuda=0;} 
    if (/^-uofm/) {$atlas="uofm"; $network = shift; chomp($network); $network =~ s/\_/\//g;}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
 
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
 
# Redirect ouput to logfile (do it only when everything is fine)
my $logfile = "$std{'DATA'}/.debug_dti_tracks.log";
open STDOUT, ">$logfile" or die "Can't redirect stdout";
open STDERR, ">&STDOUT" or die "Can't dup stdout";
$debug ? open DBG, ">$logfile" :0;
 
#Run this script ONLY on "Detritus"
#or change the paths to the correct ones
 
my $w_dir=$std{'WORKING'};
my $dti_dir=$std{'DTI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
 
#open slurm file
#my $orderfile = "$std{'DATA'}/dti_track_orders.sh";
#my $conffile = "$std{'DATA'}/dti_track_orders.conf";
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
#my $jobname = "track_$study";
 
my %plist = populate($db);
my @ok_dtis;
 
print "Collecting needed files\n";
 
my %dtis;
 
if ($cfile){
	my %cuts = get_pair($data_dir."/".$cfile);
	foreach my $cut (keys %cuts){
		if(grep {/$cut/} %plist){
			$dtis{$cut} = $plist{$cut};
		}
	}
}else{
	%dtis = %plist;
}
#open CNF, ">$conffile";
#my $count = 0;
foreach my $subject (sort keys %dtis){
	#Compulsory: el DTI debe haber pasado el preproc
	my $dti_fa = $w_dir.'/'.$dtis{$subject}.$subject.'_dti_FA.nii.gz';
	print "$dti_fa\n";
	#Compulsory: Necesito tambien la segmentacion de freesurfer de este pollo
	my $aseg = $subj_dir.'/'.$study.'_'.$dtis{$subject}.$subject.'/mri/aseg.mgz';
	print "$aseg\n";
	if((-e $dti_fa) && (-e $aseg)){
		my $order;
		unless($atlas){
			if($cuda){
				#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
				$order = $pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
			}else{
				#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
				$order = $pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
			}
		}else{
			$order = $pipe_dir."/bin/dti_bedtrack_nodes.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir." ".$pipe_dir."/lib/".$atlas."/".$network;
		}
		#$count++;
		#print CNF "$order\n";
		my $orderfile = $outdir.'/'.$subject.'dti_orders.sh';
                open ORD, ">$orderfile";
		print ORD '#!/bin/bash'."\n";
                print ORD '#SBATCH -J dti_track_'.$subject."\n";
                print ORD '#SBATCH --time=1:0:0'."\n"; #si no ha terminado en X horas matalo
                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
                print ORD '#SBATCH -o '.$outdir.'/dti_track-slurm-%j'."\n";
		if($cuda){ print ORD '#SBATCH -p cuda'."\n";}
		print ORD "srun $order\n";
                close ORD;
                my $pid = fork();
                if ($pid == 0){
                        exec("sbatch $orderfile");
                }
                elsif (!defined($pid)){
                        die "could not fork";
		}
		$debug ? print DBG "$order\n" :0;
	}
}
#close CNF;
 
#open ORD, ">$orderfile";
#print ORD '#!/bin/bash'."\n";
#print ORD '#SBATCH --time=0'."\n";
#print ORD '#SBATCH -J '.$track_study."\n";
#print ORD '#SBATCH -n '.$count."\n";
#print ORD '#SBATCH -o '.$outdir.'/dti_track-slurm-%j'."\n";
#print ORD '#SBATCH --mail-type=ALL'."\n";                                                             
#print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
#if($cuda){print ORD '#SBATCH -p cuda'."\n";}
#print ORD "srun --exclusive --multi-prog ".$conffile."\n";
#close ORD;
 
$debug ? close DBG:0;   
#exec("sbatch -o $w_dir/dti_track.out $orderfile");

Probando, en un subconjunto pequeño para que sea mas rapido,

$ cat test_bpn.csv
0002;mop
0010;mop
0011;mop
0020;mop
$ dti_track.pl -cut test_bpn.csv -uofm DMN_dorsal mopead
$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
              1080      cuda dti_trac osotolon PD       0:00      1 (Resources)
              1078     debug dti_orde osotolon  R      14:18      1 brick01
              1079      cuda dti_trac osotolon  R       0:03      1 detritus
$ ps ax | grep probtrackx2 | grep -v grep  
 4982 ?        R     17:02 probtrackx2 --opd --forcedir -s /nas/osotolongo/data/mopead/working/.tmp_mop0002/bedpostx.bedpostX/merged -m /nas/osotolongo/data/mopead/working/mop0002_dti_brain_mask -x /nas/osotolongo/data/mopead/working/.tmp_mop0002/mop0002_mask.list --dir=/nas/osotolongo/data/mopead/working/.tmp_mop0002/probtrack_out

Y ahi quedan los tractos del Default Mode Network (dorsal) de mop0002,

Nota: Si pongo solo DMN me calcula los tractos de dorsal+ventral. Esto lo he hecho asi para dar un poco mas de versatilidad.

Single mask

Si quiero una ROI compuesta puedo añadir la opción -single y el script me une todas las mascaras en una sola,

###########################################
if [ ${mode} = 0 ]; then
	for x in `awk NF $list`; do
		${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${x} -thr ${x} -div ${x} ${td}/${pollo}_mask_${x};
		echo "${td}/${pollo}_mask_${x}.nii.gz" >> ${td}/${pollo}_mask.list;
	done;
else
	count=0
	for x in `awk NF $list`; do
		${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${x} -thr ${x} -div ${x} ${td}/${pollo}_mask_${x};
		if [ ${count} = 0 ]; then
			${FSLDIR}/bin/imcp ${td}/${pollo}_mask_${x} ${td}/${pollo}_mask_seed;
		else
			${FSLDIR}/bin/fslmaths ${td}/${pollo}_mask_seed -add ${td}/${pollo}_mask_${x} ${td}/${pollo}_mask_tmp;
			${FSLDIR}/bin/imcp ${td}/${pollo}_mask_tmp ${td}/${pollo}_mask_seed;
		fi;
		((count++));
	done;
	echo "${td}/${pollo}_mask_seed.nii.gz" > ${td}/${pollo}_mask.list;
fi
###########################################

Probando,

$ dti_bedtrack_cuda.sh mopead mop0002 /nas/data/mopead/working/ 1

Tractografía (Seed and Targets)

Una de las opciones de probtrackx es tomar una ROI y calcular los tractos que hay a un grupo de regiones (seed to targets). De esta manera es posible segmentar esta ROI según las conexione existentes. Para ello deberiamos utilizar dos archivos externos y leer independientemente cada uno.

A ver como queda esto,

dti_bedtrack_seek.sh

#!/bin/sh
 
Usage() { 
    echo "" 
    echo "Usage: dti_bedtrack.sh  <project> <name> <working dir> <network dir>"
    echo "" 
    echo "You must have FSL installed in order to run this script"
    echo ""
    exit 1
} 
 
[ "$3" = "" ] && Usage
debug=1
aseg_here=0
prj=$1
shift
pollo=$1
shift
w_dir=$1
shift
net_dir=$1
shift
td=${w_dir}'/.tmp_'${pollo}
bd=${td}'/bedpostx'
list=${w_dir}/../masks.list
seed_file=${w_dir}/../dti_track.seed
targets_file=${w_dir}/../dti_track.targets
if [ ! -d "$bd" ]; then
mkdir $bd 
fi 
echo "Copying files"
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_data ${bd}/data
${FSLDIR}/bin/imcp ${w_dir}/${pollo}_dti_brain_mask ${bd}/nodif_brain_mask
cp ${td}/bvecs ${bd}/bvecs
cp ${td}/bvals ${bd}/bvals
echo "Making bedpostx"
echo [`date`]
${FSLDIR}/bin/bedpostx_gpu ${bd}
echo "So far, so good"
echo [`date`]
###########################################
echo "Getting nodes and making masks"
if [ -e "${seed_file}" ]; then
	count=0;
	for seed in `awk NF ${seed_file}`;do
		if [[ ${seed} =~ ^[0-9]+$ ]]; then
			if [ ${aseg_here} = 0 ]; then
				${PIPEDIR}/bin/get_aseg.sh ${prj} ${pollo} ${w_dir};
				aseg_here=1
			fi;
			WarpImageMultiTransform 3 ${w_dir}/${pollo}_aseg.nii.gz ${td}/${pollo}_aseg_warped.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt;
			${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${seed} -thr ${seed} -div ${seed} ${td}/${pollo}_mask_${seed};
			if [ "${count}" -gt "0" ]; then
				${FSLDIR}/bin/fslmaths ${td}/${pollo}_mask_seed -add ${td}/${pollo}_mask_${seed} ${td}/${pollo}_mask_tmp;
				${FSLDIR}/bin/imcp ${td}/${pollo}_mask_tmp ${td}/${pollo}_mask_seed;
			else
				${FSLDIR}/bin/imcp ${td}/${pollo}_mask_${seed} ${td}/${pollo}_mask_seed;
			fi
			((count++))
		else
			x=`find ${net_dir} -name "${seed}"`;
			node=`${FSLDIR}/bin/remove_ext $(basename ${x})`;
			WarpImageMultiTransform 3 ${x} ${td}/${node}_warped.nii.gz -R ${w_dir}/${pollo}_t1_reoriented.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Warp.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Affine.txt;
			WarpImageMultiTransform 3 ${td}/${node}_warped.nii.gz ${td}/${pollo}_${node}.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Azffine.txt;
			if [ $count ]; then
				${FSLDIR}/bin/fslmaths ${td}/${pollo}_mask_seed -add ${td}/${pollo}_${node} ${td}/${pollo}_mask_tmp;
				${FSLDIR}/bin/imcp ${td}/${pollo}_mask_tmp ${td}/${pollo}_mask_seed;
			else
				${FSLDIR}/bin/imcp ${td}/${pollo}_${node} ${td}/${pollo}_mask_seed;
			fi
			((count++))
		fi;
	done;
	echo "${td}/${pollo}_mask_seed.nii.gz" > ${td}/${pollo}_seed.list;
else
	exit 1 "No seed file!";
fi
if [ -e "${targets_file}" ]; then
	for target in `awk NF ${targets_file}`; do
		if [[ ${target} =~ ^[0-9]+$ ]]; then
			if [ ! ${aseg_here} ]; then
				${PIPEDIR}/bin/get_aseg.sh ${prj} ${pollo} ${w_dir};
				aseg_here=1
			fi;
			WarpImageMultiTransform 3 ${w_dir}/${pollo}_aseg.nii.gz ${td}/${pollo}_aseg_warped.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt;
			${FSLDIR}/bin/fslmaths ${td}/${pollo}_aseg_warped -uthr ${target} -thr ${target} -div ${target} ${td}/${pollo}_mask_${target};
			echo "${td}/${pollo}_mask_${target}.nii.gz" >> ${td}/${pollo}_targets.list;
		else
			x=`find ${net_dir} -name "${target}"`;
			node=$(basename ${x});
			WarpImageMultiTransform 3 ${x} ${td}/${node%.nii.gz}_warped.nii.gz -R ${w_dir}/${pollo}_t1_reoriented.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Warp.nii.gz ${td}/${pollo}_dti_ants_elast_mni_t1Affine.txt;
			WarpImageMultiTransform 3 ${td}/${node%.nii.gz}_warped.nii.gz ${td}/${pollo}_${node}.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Warp.nii.gz ${td}/${pollo}_dti_ants_elast_t1_b0Affine.txt;
			echo "${td}/${pollo}_${node}.gz" >> ${td}/${pollo}_targets.list;
		fi;
	done;
else
	exit 1 "No targets file!";
fi
###########################################		
echo "Doing probtrackx"
probtrackx2 --opd --os2t --forcedir -s ${bd}.bedpostX/merged -m ${w_dir}/${pollo}_dti_brain_mask -x ${td}/${pollo}_seed.list --targetmasks=${td}/${pollo}_targets.list --dir=${td}/probtrack_out
rm ${td}/${pollo}_seed.list;
rm ${td}/${pollo}_targets.list;
echo "Done"
echo [`date`]

Lo he hecho de tal manera que lee el archivo dti_track.seed del directorio del proyecto y calcula una ROI compuesta con la lista que haya en el archivo. Los targets los adquiere del archivo dti_track.targets, en el mismo directorio. Ambos archivos son comunes a todos los sujetos.

ejemplo de los archivos para la segmentacion de FS

$ cat dti_track.seed
10
49
 
$ cat dti_track.targets 
11
12
13
17
18
19
20
30
31
50
51
52
53
54
55
56
58
60
62
63
85

Vamos a probarlo,

$ dti_bedtrack_seek.sh mopead mop0002 /nas/data/mopead/working/

Demora un rato y,ademas del archivo con todos los tractos, devuelve un archivo por cada target.

$ ls working/.tmp_mop0002/probtrack_out/
biggest.nii.gz                   seeds_to_mop0002_mask_20.nii.gz  seeds_to_mop0002_mask_56.nii.gz
fdt_paths.nii.gz                 seeds_to_mop0002_mask_30.nii.gz  seeds_to_mop0002_mask_58.nii.gz
probtrackx.log                   seeds_to_mop0002_mask_31.nii.gz  seeds_to_mop0002_mask_60.nii.gz
seeds_to_mop0002_mask_11.nii.gz  seeds_to_mop0002_mask_50.nii.gz  seeds_to_mop0002_mask_62.nii.gz
seeds_to_mop0002_mask_12.nii.gz  seeds_to_mop0002_mask_51.nii.gz  seeds_to_mop0002_mask_63.nii.gz
seeds_to_mop0002_mask_13.nii.gz  seeds_to_mop0002_mask_52.nii.gz  seeds_to_mop0002_mask_85.nii.gz
seeds_to_mop0002_mask_17.nii.gz  seeds_to_mop0002_mask_53.nii.gz  waytotal
seeds_to_mop0002_mask_18.nii.gz  seeds_to_mop0002_mask_54.nii.gz
seeds_to_mop0002_mask_19.nii.gz  seeds_to_mop0002_mask_55.nii.gz

Aqui podemos usar find_the_biggest para sacar la segmentacion,

$ find_the_biggest working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_* working/.tmp_mop0002/probtrack_out/biggest
number of inputs 21
Indices
1 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_11.nii.gz
2 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_12.nii.gz
3 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_13.nii.gz
4 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_17.nii.gz
5 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_18.nii.gz
6 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_19.nii.gz
7 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_20.nii.gz
8 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_30.nii.gz
9 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_31.nii.gz
10 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_50.nii.gz
11 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_51.nii.gz
12 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_52.nii.gz
13 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_53.nii.gz
14 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_54.nii.gz
15 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_55.nii.gz
16 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_56.nii.gz
17 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_58.nii.gz
18 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_60.nii.gz
19 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_62.nii.gz
20 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_63.nii.gz
21 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mask_85.nii.gz

y ahi queda el talamo segmentado por la cantidad de conexiones,

Vamos a intentar otra cosa,

$ ls /opt/neuro.dev/lib/uofm/LN/ > dti_track.targets
 
$ cat dti_track.targets 
mri_LN1.nii
mri_LN2.nii
mri_LN3.nii
mri_LN4.nii
mri_LN5.nii
mri_LN6.nii
mri_LN7.nii
 
$ dti_bedtrack_seek.sh mopead mop0002 /nas/data/mopead/working/ /opt/neuro.dev/lib/uofm/
 
$ find_the_biggest working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN* working/.tmp_mop0002/probtrack_out/biggest
number of inputs 7
Indices
1 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN1.nii.gz
2 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN2.nii.gz
3 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN3.nii.gz
4 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN4.nii.gz
5 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN5.nii.gz
6 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN6.nii.gz
7 working/.tmp_mop0002/probtrack_out/seeds_to_mop0002_mri_LN7.nii.gz

y tenemos la nueva segmentacion del talamo segun sus conexiones a la Language Network. Que no parece gran cosa pero es una prueba de concepto de que podemos usar varios atlas en el mismo experimento,

OK, ahora a meter todo esto en el wrapper.

dti_track.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
my $debug=1;
my $corr=1;
my $legacy=0;
my $cuda=1;
my $atlas=0; 
my $mode=0;
my $s2t=0;
my $network;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-slow/) {$cuda=0;}
    if (/^-single/) {$mode=1;}
    if (/^-s2t/) {$s2t=1;}
    if (/^-uofm/) {$atlas="uofm"; $network = shift; chomp($network); $network =~ s/\_/\//g;}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_reg.hlp'; exit;}
 
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
 
# Redirect ouput to logfile (do it only when everything is fine)
my $logfile = "$std{'DATA'}/.debug_dti_tracks.log";
open STDOUT, ">$logfile" or die "Can't redirect stdout";
open STDERR, ">&STDOUT" or die "Can't dup stdout";
$debug ? open DBG, ">$logfile" :0;
 
#Run this script ONLY on "Detritus"
#or change the paths to the correct ones
 
my $w_dir=$std{'WORKING'};
my $dti_dir=$std{'DTI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
 
#open slurm file
#my $orderfile = "$std{'DATA'}/dti_track_orders.sh";
#my $conffile = "$std{'DATA'}/dti_track_orders.conf";
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
#my $jobname = "track_$study";
 
my %plist = populate($db);
my @ok_dtis;
 
print "Collecting needed files\n";
 
my %dtis;
 
if ($cfile){
	my %cuts = get_pair($data_dir."/".$cfile);
	foreach my $cut (keys %cuts){
		if(grep {/$cut/} %plist){
			$dtis{$cut} = $plist{$cut};
		}
	}
}else{
	%dtis = %plist;
}
#open CNF, ">$conffile";
#my $count = 0;
foreach my $subject (sort keys %dtis){
	#Compulsory: el DTI debe haber pasado el preproc
	my $dti_fa = $w_dir.'/'.$dtis{$subject}.$subject.'_dti_FA.nii.gz';
	print "$dti_fa\n";
	#Compulsory: Necesito tambien la segmentacion de freesurfer de este pollo
	my $aseg = $subj_dir.'/'.$study.'_'.$dtis{$subject}.$subject.'/mri/aseg.mgz';
	print "$aseg\n";
	if((-e $dti_fa) && (-e $aseg)){
		my $order;
		unless($atlas){
			if($cuda){
				if($s2t){
					$order = $pipe_dir."/bin/dti_bedtrack_seek.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
				}else{
					#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
					$order = $pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir." ".$mode;
				}else{
					#$order = $count."\t".$pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir;
					$order = $pipe_dir."/bin/dti_bedtrack.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir." ".$mode;
				}
			}
		}else{
			if($s2t){
				$order = $pipe_dir."/bin/dti_bedtrack_seek.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir." ".$pipe_dir."/lib/".$atlas;
			}else{
				$order = $pipe_dir."/bin/dti_bedtrack_nodes.sh ".$study." ".$dtis{$subject}.$subject." ".$w_dir." ".$pipe_dir."/lib/".$atlas."/".$network." ".$mode;
			}
		}
		#$count++;
		#print CNF "$order\n";
		my $orderfile = $outdir.'/'.$subject.'dti_orders.sh';
                open ORD, ">$orderfile";
		print ORD '#!/bin/bash'."\n";
                print ORD '#SBATCH -J dti_track_'.$subject."\n";
                print ORD '#SBATCH --time=1:0:0'."\n"; #si no ha terminado en X horas matalo
                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
                print ORD '#SBATCH -o '.$outdir.'/dti_track-slurm-%j'."\n";
		if($cuda){ print ORD '#SBATCH -p cuda'."\n";}
		print ORD "srun $order\n";
                close ORD;
                my $pid = fork();
                if ($pid == 0){
                        exec("sbatch $orderfile");
                }
                elsif (!defined($pid)){
                        die "could not fork";
		}
		$debug ? print DBG "$order\n" :0;
	}
}
#close CNF;
 
#open ORD, ">$orderfile";
#print ORD '#!/bin/bash'."\n";
#print ORD '#SBATCH --time=0'."\n";
#print ORD '#SBATCH -J '.$track_study."\n";
#print ORD '#SBATCH -n '.$count."\n";
#print ORD '#SBATCH -o '.$outdir.'/dti_track-slurm-%j'."\n";
#print ORD '#SBATCH --mail-type=ALL'."\n";                                                             
#print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
#if($cuda){print ORD '#SBATCH -p cuda'."\n";}
#print ORD "srun --exclusive --multi-prog ".$conffile."\n";
#close ORD;
 
$debug ? close DBG:0;   
#exec("sbatch -o $w_dir/dti_track.out $orderfile");

seek and destroy

Alt metrics

Una vez que tenemos sacada las fibras de determinada red de interes, podriamos superponer esta red encima de los resultados de dtifit y calcular la FA y MD media en esas fibras.

La idea general es que probtrackx deja en el archivo probtrack_out/fdt_paths.nii.gz la informacion de las fibrasque existen entre los nodos de la red que le hemos pedido que busque. En principio, podemos conocer facilmente cual es el valor de la FA si usamos este archivo como mascara. Para eliminar un poco el ruido vamos a definir un threshold sobre el cual cortar los valores de la imagen y el resto lo binarizamos. Una vez hecho esto se prodece como con cualquier metrica que calculamos aqui.

Hala, primero a hacer la mascara,

track2mask.sh

#!/bin/sh
Usage() {
    echo ""
    echo "Usage: track2mask.sh  <file> <thr>"
    echo ""
    echo "You must have FSL installed in order to run this script"
    echo ""
    exit 1
}
 
[ "$1" = "" ] && Usage
pollo=`${FSLDIR}/bin/remove_ext $1`
shift
thr=$1
shift
max=($(${FSLDIR}/bin/fslstats ${pollo} -r | awk '{print $2}'))
cut=$(echo $max*$thr | bc)
${FSLDIR}/bin/fslmaths ${pollo} -thr ${cut} -bin ${pollo}_mask
exit 0

Esto es sencillo,tomamos un archivo , hallamos el valor maximo con fslstats, le sacamos un % determinado y cortamos y binarizamos con fslmaths. ¿escoger threshold?

Ahora, para tener el valor medio y SD de la FA en esta region habria que hacer algo asi,

$ track2mask.sh working/smc0209_probtrack_DMN/fdt_paths.nii.gz 0.25
$ fslstats working/smc0209_dti_FA -k working/smc0209_probtrack_DMN/fdt_paths_mask -M -S

Pos hacemos tablita con esto y ya ta.

A ponerlo bonito

dti_metrics_alt.pl

#!/usr/bin/perl
 
use strict; use warnings;
 
use File::Find::Rule;
use NEURO qw(print_help get_pair load_study achtung shit_done);
use Data::Dump qw(dump);
use File::Basename qw(basename);
 
my $attach = 1;
my $with_sd = 0;
my $threshold = 0.25;
#my $atlas;
my $net_path;
#my $netwstr;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-l/) { $attach = 0;}
    #if (/^-uofm/) {$atlas="uofm"; $network = shift; chomp($network); $netwstr = $network; $network =~ s/\_/\//g;}
    if (/^-thr/) {$threshold = shift;}
    if (/^-path/) {$net_path = shift;}
    if (/^-sd/) { $with_sd = 1;}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_metrics.hlp'; exit;}
}
 
my $study = shift;
$net_path = "out" unless $net_path;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_metrics.hlp'; exit;}
my %std = load_study($study);
 
my $w_dir=$std{'WORKING'};
my $data_dir=$std{'DATA'};
 
my @dti_fas = find(file => 'name' => "*_dti_FA.nii.gz", in => $w_dir);
my %csv;
 
my $ofile = $data_dir."/".$study."_dti_".$net_path.".csv";
open OF, ">$ofile";
print OF "Subject";
if($with_sd){
	print OF ";$net_path","_FA_Mean;","$net_path","_FA_STD",";$net_path","_MD_Mean;","$net_path","_MD_STD";
}else{
	print OF ";$net_path","_FA",";$net_path","_MD";
}
print OF "\n";
 
foreach my $dti (sort @dti_fas){
	(my $a,my $subject) = $dti =~ /(\w{1,4})(\d{1,6})_dti_FA/;
	my $tmpmd = $w_dir.'/'.$a.$subject."_dti_MD";
	if($subject){
		print OF "$subject";
		my $img_path = $w_dir.'/'.$a.$subject.'_probtrack_'.$net_path.'/fdt_paths.nii.gz'; 
		my $order = $ENV{'PIPEDIR'}.'/bin/track2mask.sh '.$img_path.' '.$threshold;
		print "$order\n";
		system($order);
		my $mask_path = $w_dir.'/'.$a.$subject.'_probtrack_'.$net_path.'/fdt_paths_mask.nii.gz'; 
		$order = "fslstats ".$dti." -k ".$mask_path." -M -S";
		print "$order\n";
		(my $mean, my $std) = map{/(\d+\.\d*)\s*(\d+\.\d*)/} qx/$order/;
		if($with_sd){
			print OF ";$mean",";$std";
		}else{
			print OF ";$mean";
		}
		$order = "fslstats ".$tmpmd." -k ".$mask_path." -M -S";
		print "$order\n";
		($mean, $std) = map{/(\d+\.\d*)\s*(\d+\.\d*)/} qx/$order/;
		if($with_sd){
			print OF ";$mean",";$std";
		}else{
			print OF ";$mean";
		}
	}
	print OF "\n";
}
 
close OF;
 
my $zfile = $ofile.'.gz';
system("gzip -c $ofile > $zfile");
 
if ($attach){
	shit_done basename($ENV{_}), $study, $zfile;
}else{
	achtung basename($ENV{_}), $ofile, $study;
}

trick or treat

La cosa es que esto hay que hacerlo con un poco de cuidado. La forma de proceder es esta (ejemplo con el facehbi),

$ dti_track.pl -uofm DMN facehbi
$ for x in `ls -d working/smc0*_probtrack_out`; do mv $x `echo $x | sed 's/out/DMN/'`;done
$ dti_metrics_alt.pl -path DMN facehbi

y miralo,

$ head facehbi_dti_DMN.csv
Subject;DMN_FA;DMN_MD
001;;
002;;
003;;
004;0.295167;0.001053
005;0.335169;0.000928
006;0.306471;0.001024
007;0.335922;0.000870
008;0.293183;0.000994
009;0.343917;0.000898

Estrictamente, no hay que mover los paths. Si no se le pasa ninguno, el script los busca en out, que es donde estan por defecto. Pero si no te organizas aqui te vas a volver loco. Tambien por defecto el threshold se toma del 25% pero se puede cambiar con la opcion -thr. El valor se pasa con base a 1, no a 100.

fMRI

Mirando los DICOMs

Los fMRI vienen con la etiqueta ep2d_bold_p2_resting_state,

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep ep2d_bold_p2_resting_state` ]]; then dirname $x; fi;done | uniq
/nas/corachan/MOPEAD/52JAYALW/00000000/00000011
/nas/corachan/MOPEAD/52JAYALW/00000000/00000012

Vaya, que hay 2! Vamos a convertirlos. Por partes para no meter la gamba,

$ for x in `find /nas/corachan/MOPEAD/52JAYALW/ -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep ep2d_bold_p2_resting_state` ]]; then dirname $x; fi;done | uniq > tmp_list.dirs
$ for x in `cat tmp_list.dirs`; do dcm2niix -z y -o tmp $x; done
Chris Rorden's dcm2niiX version v1.0.20180622 (JP2:OpenJPEG) (JP-LS:CharLS) GCC5.5.0 (64-bit Linux)
Found 200 DICOM file(s)
slices stacked despite varying acquisition numbers (if this is not desired please recompile)
Convert 200 DICOM as tmp/00000011_ep2d_bold_p2_resting_state_20180705100054_12 (64x64x38x200)
compress: "/usr/local/mricron/pigz_mricron" -n -f -6 "tmp/00000011_ep2d_bold_p2_resting_state_20180705100054_12.nii"
Conversion required 5.829411 seconds (0.630000 for core code).
Chris Rorden's dcm2niiX version v1.0.20180622 (JP2:OpenJPEG) (JP-LS:CharLS) GCC5.5.0 (64-bit Linux)
Found 200 DICOM file(s)
Unable to determine spatial orientation: 0020,0037 missing!
Warning: Unable to determine slice direction: please check whether slices are flipped
Warning: Bogus spatial matrix (perhaps non-spatial image): inspect spatial orientation
Convert 200 DICOM as tmp/00000012_ep2d_bold_p2_resting_state_20180705100054_13 (448x448x200x1)
compress: "/usr/local/mricron/pigz_mricron" -n -f -6 "tmp/00000012_ep2d_bold_p2_resting_state_20180705100054_13.nii"
Conversion required 12.211212 seconds (0.580000 for core code).
$ rm tmp_list.dirs
$ ls -lh tmp
total 63M
-rw-rw---- 1 osotolongo osotolongo 1.2K Oct 26 10:36 00000011_ep2d_bold_p2_resting_state_20180705100054_12.json
-rw-rw---- 1 osotolongo osotolongo  32M Oct 26 10:36 00000011_ep2d_bold_p2_resting_state_20180705100054_12.nii.gz
-rw-rw---- 1 osotolongo osotolongo  502 Oct 26 10:36 00000012_ep2d_bold_p2_resting_state_20180705100054_13.json
-rw-rw---- 1 osotolongo osotolongo  32M Oct 26 10:36 00000012_ep2d_bold_p2_resting_state_20180705100054_13.nii.gz

Vamos a mirarlos,

$ fslinfo tmp/00000011_ep2d_bold_p2_resting_state_20180705100054_12.nii.gz 
data_type      INT16
dim1           64
dim2           64
dim3           38
dim4           200
datatype       4
pixdim1        3.296880
pixdim2        3.296880
pixdim3        3.300000
pixdim4        2.130000
cal_max        0.0000
cal_min        0.0000
file_type      NIFTI-1+
$ fslinfo tmp/00000012_ep2d_bold_p2_resting_state_20180705100054_13.nii.gz 
data_type      INT16
dim1           448
dim2           448
dim3           200
dim4           1
datatype       4
pixdim1        1.000000
pixdim2        1.000000
pixdim3        1.000000
pixdim4        0.000000
cal_max        0.0000
cal_min        0.0000
file_type      NIFTI-1+

Vaya, pues no son iguales. A ver como se ven,

Pues el segundo pinta mal, no sirve para nada. Luego hay que seleccionar solo el primero,lo que es facil usando dim4>1.

Además, necesitaré los ep2d_fid_basic_bold_p2_AP y ep2d_fid_basic_bold_p2_PA para corregir las distorsiones. Con estos no parece que haya problemas.

Vamos alla

update_fmri.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study);
use File::Find::Rule;
use Data::Dump qw(dump);
 
my $study = "mopead";
my %std = load_study($study);
 
my $src_dir = '/nas/corachan/MOPEAD/';
#my $raw_dir = $std{'DATA'}.'/raw';
my $conv_file = $std{'DATA'}.'/converted_fmri.csv';
my $ids_file = $std{'DATA'}.'/ids.csv';
my $names_file =$std{'DATA'}.'/ni_names.csv';
my $proj_file = $std{'DATA'}.'/mopead.csv';
my @other_exts = ("nii.gz","json");
 
#Leo la DB
my %idsinfo;
open IDS, "<$ids_file" or die $!;
while(<IDS>){
        chomp;
        my ($key, $value) = split(/,\s?/);
        $idsinfo{$key} = $value;
}
close IDS;
#Leo los que ya he descomprimido
my @orig_str = qx/ls $src_dir/;
chomp @orig_str;
#Leo los que se han convertido
my @conv_str;
open CONV, "<$conv_file" or die $!;
chomp (@conv_str = <CONV>);
close CONV;
#Y a ver cuantos faltan por convertir!
my %conv_str = map { $_ => 1 } @conv_str;
my @not_conv = grep !$conv_str{$_}, @orig_str;
# y a convertirlos!
foreach my $guy (@not_conv){
	my $order = 'for x in `find '.$src_dir.'/'.$guy.' -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep ep2d_bold_p2_resting_state` ]]; then dirname $x; fi; done | uniq > '.$std{'DATA'}.'/tmp/tmp_list.dirs; for x in `cat '.$std{'DATA'}.'/tmp/tmp_list.dirs`; do dcm2niix -z y -o '.$std{'DATA'}.'/tmp $x; done; rm '.$std{'DATA'}.'/tmp/tmp_list.dirs';
	print "$order\n";
	system($order);
	my @conv_files = find(file => 'name' => "*.nii.gz", in => $std{'DATA'}.'/tmp/');
	#dump @conv_files;
	foreach my $nii_file (@conv_files){
		$order = 'fslinfo '.$nii_file;
		my %xinfo;
		foreach (qx/$order/){
			my ( $key, $value ) = /(\S+)\s+(\S+)\s*.*/;
			$xinfo{$key} =$value;
		}
		if($xinfo{"dim4"}>1){
			print "Choosing and moving files\n";
			my $file_name = $nii_file;
			foreach my $ext (@other_exts){
				($file_name = $nii_file) =~ s/nii.gz$/$ext/;
				$order = sprintf "mv %s %s/fmri/mop%04ds0012.%s", $file_name, $std{'DATA'}, $idsinfo{$guy}, $ext;
				print "$order\n";
				system($order);
			}			
		}
	}
	print "Cleaning house\n";
	$order = "rm $std{'DATA'}/tmp/*";
	system($order);
	#Convierto las imagenes para corregir (A>>P y P>>A)
	print "Getting A>>P\n";
	$order = 'for x in `find '.$src_dir.'/'.$guy.' -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep ep2d_fid_basic_bold_p2_AP` ]]; then dcm2niix -z y -o '.$std{'DATA'}.'/tmp $x; fi;done';
	system($order);
	$order = sprintf "mv %s/tmp/*.nii.gz %s/fmri/mop%04d_a2p.nii.gz", $std{'DATA'}, $std{'DATA'}, $idsinfo{$guy};
	print "$order\n";
        system($order);
	print "Getting P>>A\n";                                                                       
        $order = 'for x in `find '.$src_dir.'/'.$guy.' -type f`; do if [[ `dckey -k "SeriesDescription" $x 2>&1 | grep ep2d_fid_basic_bold_p2_PA` ]]; then dcm2niix -z y -o '.$std{'DATA'}.'/tmp $x; fi;done';
	system($order);
        $order = sprintf "mv %s/tmp/*.nii.gz %s/fmri/mop%04d_p2a.nii.gz", $std{'DATA'}, $std{'DATA'}, $idsinfo{$guy};                                                                                       
        print "$order\n";                                                                             
        system($order);
        print "Cleaning house\n";                                                          
        $order = "rm $std{'DATA'}/tmp/*";
	system($order); 
}
 
#actualizo los convertidos
push(@conv_str, @not_conv);
open CONV, ">$conv_file" or die $!;
print CONV "$_\n" foreach @conv_str;
close CONV;

A ver si funciona,

$ touch converted_fmri.csv
$ ./update_fmri.pl 
blah blah blah
$ ls fmri
mop0001_a2p.nii.gz   mop0006s0012.json    mop0012s0012.json    mop0019_a2p.nii.gz
mop0001_p2a.nii.gz   mop0006s0012.nii.gz  mop0012s0012.nii.gz  mop0019_p2a.nii.gz
mop0001s0012.json    mop0007_a2p.nii.gz   mop0013_a2p.nii.gz   mop0019s0012.json
mop0001s0012.nii.gz  mop0007_p2a.nii.gz   mop0013_p2a.nii.gz   mop0019s0012.nii.gz
mop0002_a2p.nii.gz   mop0007s0012.json    mop0014_a2p.nii.gz   mop0020_a2p.nii.gz
mop0002_p2a.nii.gz   mop0007s0012.nii.gz  mop0014_p2a.nii.gz   mop0020_p2a.nii.gz
mop0002s0012.json    mop0008_a2p.nii.gz   mop0014s0012.json    mop0020s0012.json
mop0002s0012.nii.gz  mop0008_p2a.nii.gz   mop0014s0012.nii.gz  mop0020s0012.nii.gz
mop0003_a2p.nii.gz   mop0008s0012.json    mop0015_a2p.nii.gz   mop0021_a2p.nii.gz
mop0003_p2a.nii.gz   mop0008s0012.nii.gz  mop0015_p2a.nii.gz   mop0021_p2a.nii.gz
mop0003s0012.json    mop0009_a2p.nii.gz   mop0015s0012.json    mop0021s0012.json
mop0003s0012.nii.gz  mop0009_p2a.nii.gz   mop0015s0012.nii.gz  mop0021s0012.nii.gz
mop0004_a2p.nii.gz   mop0009s0012.json    mop0016_a2p.nii.gz   mop0022_a2p.nii.gz
mop0004_p2a.nii.gz   mop0009s0012.nii.gz  mop0016_p2a.nii.gz   mop0022_p2a.nii.gz
mop0004s0012.json    mop0010_a2p.nii.gz   mop0017_a2p.nii.gz   mop0022s0012.json
mop0004s0012.nii.gz  mop0010_p2a.nii.gz   mop0017_p2a.nii.gz   mop0022s0012.nii.gz
mop0005_a2p.nii.gz   mop0010s0012.json    mop0017s0012.json    mop0023_a2p.nii.gz
mop0005_p2a.nii.gz   mop0010s0012.nii.gz  mop0017s0012.nii.gz  mop0023_p2a.nii.gz
mop0005s0012.json    mop0011_a2p.nii.gz   mop0018_a2p.nii.gz   mop0023s0012.json
mop0005s0012.nii.gz  mop0011_p2a.nii.gz   mop0018_p2a.nii.gz   mop0023s0012.nii.gz
mop0006_a2p.nii.gz   mop0012_a2p.nii.gz   mop0018s0012.json
mop0006_p2a.nii.gz   mop0012_p2a.nii.gz   mop0018s0012.nii.gz

Y listo para empezar a trabajar

seleccion de series fMRI a partir de los NIfTIs

Aqui se han de copiar tres archivos, ep2d_bold_p2_resting_state, ep2d_fid_basic_bold_p2_AP, ep2d_fid_basic_bold_p2_PA.

Veamos como organizarlo

update_fmri.pl

#!/usr/bin/perl
 
use strict; use warnings;
use NEURO qw(load_study check_or_make);
use File::Find::Rule;
use Data::Dump qw(dump);
 
my $study = "mopead";
my $prefix = "mop";
my %std = load_study($study);
#unless ($itype) {$itype='t1';}
#my %dtype = ("t1" => "mri", );
my $src_dir = '/nas/corachan/MOPEAD/';
#my $raw_dir = $std{'DATA'}.'/raw';
my $conv_file = $std{'DATA'}.'/converted_fmri.csv';
my $ids_file = $std{'DATA'}.'/ids.csv';
my $names_file =$std{'DATA'}.'/ni_names.csv';
my $proj_file = $std{'DATA'}.'/'.$study.'.csv';
my @other_exts = ("nii.gz","json");
my $series = 12;
my $tag = "ep2d_bold_p2_resting_state";
my $aptag = "ep2d_fid_basic_bold_p2_AP";
my $patag = "ep2d_fid_basic_bold_p2_PA";
check_or_make("$std{'DATA'}/fmri");
#Leo la DB
my %idsinfo;
open IDS, "<$ids_file" or die $!;
while(<IDS>){
        chomp;
        my ($key, $value) = split(/,\s?/);
        $idsinfo{$key} = $value;
}
close IDS;
#Leo los que ya he descomprimido
my @orig_str = qx/ls $src_dir/;
chomp @orig_str;
#Leo los que se han convertido
my @conv_str;
open CONV, "<$conv_file" or die $!;
chomp (@conv_str = <CONV>);
close CONV;
#Y a ver cuantos faltan por convertir!
my %conv_str = map { $_ => 1 } @conv_str;
my @not_conv = grep !$conv_str{$_}, @orig_str;
# y a convertirlos!
my @newones;
foreach my $guy (@not_conv){
        if(exists($idsinfo{$guy})){
                my $count = $series;
                my @conv_files = find(file => 'name' => "*$tag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                #dump @conv_files;
                foreach my $nii_file (@conv_files){
                        my $order = 'fslval '.$nii_file.' dim4';
                        my $xinfo = qx/$order/;
                        if($xinfo>1){
                                print "Choosing and moving files\n";
                                my $file_name = $nii_file;
                                foreach my $ext (@other_exts){
                                        ($file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                        $order = sprintf "cp %s %s/%s%04ds%04d.%s", $file_name, $std{'fMRI'}, $prefix, $idsinfo{$guy}, $count, $ext;
                                        print "$order\n";
                                        system($order);
                                }
                        push @newones, $guy;
                        $count++;                       
                        }
                }
                #Convierto lasimagenes para corregir (A>>P y P>>A)
                print "Getting A>>P\n";
                my $count = $series;
                @conv_files = find(file => 'name' => "*$aptag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                foreach my $nii_file (@conv_files){
                        foreach my $ext (@other_exts){
                                (my $file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                my $order = sprintf "cp %s %s/%s%04ds%04d_a2p.%s", $file_name, $std{'fMRI'}, $prefix, $idsinfo{$guy}, $count, $ext;
                                print "$order\n";
                                system($order);
                        }
                        $count++;
                }
                print "Getting P>>A\n";
                my $count = $series;
                @conv_files = find(file => 'name' => "*$patag*.nii.gz", in => $std{'DATA'}.'/niftis/'.$guy.'/');
                foreach my $nii_file (@conv_files){
                        foreach my $ext (@other_exts){
                                (my $file_name = $nii_file) =~ s/nii.gz$/$ext/;
                                my $order = sprintf "cp %s %s/%s%04ds%04d_p2a.%s", $file_name, $std{'fMRI'}, $prefix, $idsinfo{$guy}, $count, $ext;
                                print "$order\n";
                                system($order);
                        }
                        $count++;
                }
 
        }
}
 
#actualizo los convertidos
push(@conv_str, @newones);
open CONV, ">$conv_file" or die $!;
print CONV "$_\n" foreach @conv_str;
close CONV;

Veamos el output,

[osotolongo@detritus mopead]$ ls fmri/mop0001*
fmri/mop0001s0012_a2p.json    fmri/mop0001s0012.nii.gz      fmri/mop0001s0013_a2p.json    fmri/mop0001s0013.nii.gz      fmri/mop0001s0014_a2p.json    fmri/mop0001s0014.nii.gz
fmri/mop0001s0012_a2p.nii.gz  fmri/mop0001s0012_p2a.json    fmri/mop0001s0013_a2p.nii.gz  fmri/mop0001s0013_p2a.json    fmri/mop0001s0014_a2p.nii.gz  fmri/mop0001s0014_p2a.json
fmri/mop0001s0012.json        fmri/mop0001s0012_p2a.nii.gz  fmri/mop0001s0013.json        fmri/mop0001s0013_p2a.nii.gz  fmri/mop0001s0014.json        fmri/mop0001s0014_p2a.nii.gz

Nota: Tengo el problema de tener tres series distintas dond esolo deberia haber una y aparentemente no hay forma de diferenciarlas. Preguntar esto!

Independent Component Analysis (AKA ICA)

Para hacer un ICA en FSL se utiliza la herramienta FEAT. Esto es una interfaz que lanza varios script de FSL, creando los archivos necesarios para lanzar MELODIC. Esto es bueno porque acorta el procedimiento pero tiene la dificultad añadida de que la estructura de archivos no está documentada y practicamente hay que hacer ingenieria inversa para que funcione si no se lanza con la gui o con fsl_sub.

Para hacer un ICA individual hacemos algo asi,

$ ${FSLDIR}/bin/feat /path/design.fsf -D /output/path/ica -I 1 -init
$ ${FSLDIR}/bin/feat /path/design.fsf -D /output/path/ica -I 1 -prestats
$ ${FSLDIR}/bin/melodic -i /output/path/ica/filtered_func_data -o /output/path/ica/filtered_func_data.ica -v --nobet --bgthreshold=3 --tr=2.130000 --report -d 0 --mmthresh=0.5 --Ostats

El archivo design.fsf contiene todos los parametros necesarios. Para llenarlo se usa feat_gui. He guardado una plantilla en ${PIPEDIR}/lib/lone_ica_template.fsf. Los parametros a cambiar en cada sujeto son,

# Output directory
set fmri(outputdir) <output_dir> 

# TR(s)
set fmri(tr) <TR>

# Total volumes
set fmri(npts) <number_of_slices>

# Total voxels
set fmri(totalVoxels) <size_of_image>

# 4D AVW data or FEAT directory (1)
set feat_files(1) <data_input>

Aqui, <ouput_dir> es donde se van a guardar los resultados del registros y de MELODIC y <data_input> es el archivo fmri original. El resto de los datos han de ser inferidos del archivo a procesar (con fslinfo deberia funcionar).

Primero invocamos fslinfo para sacar y/o calcular los datos,

		system("$fsl'/bin/imcp' $names[0] $data");
                my $test_order = $fsl.'/bin/fslinfo '.$data;
                open(my $test_info, '-|', $test_order ) or die $!;
                my %test_data;
                while (<$test_info>){
                        if (/^dim1\s*(\d{1,3}).*/) {$test_data{'dim1'}=$1;}
                        if (/^dim2\s*(\d{1,3}).*/) {$test_data{'dim2'}=$1;}
                        if (/^dim3\s*(\d{1,3}).*/) {$test_data{'dim3'}=$1;}
                        if (/^dim4\s*(\d{1,3}).*/) {$test_data{'dim4'}=$1;}
                        if (/^pixdim4\s*(\d{1,2}\.\d{1,6})/) {$test_data{'tr'}=$1;}
                }               
                my $data_size = $test_data{'dim1'}*$test_data{'dim2'}*$test_data{'dim3'}*$test_data{'dim4'};

Y para cada sujeto se ha de rellenar la plantilla,

		my $dsg_file = $w_dir.'/'.$fmris{$subject}.$subject.'_design.fsf';
		open TPF,"<$template_file";
		open DSG,">$dsg_file";
		while(<TPF>){
			s/<output_dir>/\"$output_dir\"/;
			s/<data_input>/\"$data\"/;
                        s/<TR>/$test_data{'tr'}/;
                        s/<number_of_slices>/$test_data{'dim4'}/;
                        s/<size_of_image>/$data_size/;
			print DSG;
		}
		close DSG;
		close TPF;

y luego preparar los scripts en slurm, de manera que cada uno se ejecute si el siguiente ha terminado satisfactoriamente. Algo así,

sbatch /nas/osotolongo/data/mopead/working/mop0001_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:2169 /nas/osotolongo/data/mopead/working/mop0001_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:2170 /nas/osotolongo/data/mopead/working/mop0001_rs.ica/scripts/feat4.sh

Para esto voy a hacer tambien una plantilla de los scripts y los relleno para cada sujeto,

$ cat ${PIPEDIR}/lib/script_templates/feat1.template 
#!/bin/bash
#SBATCH -J feat1-<subject>
#SBATCH --time=00:10:00
#SBATCH --mail-type=FAIL,STAGE_OUT
#SBATCH --mail-user=<mailto>
#SBATCH -o <out_dir>/feat1-slurm-%j
${FSLDIR}/bin/feat <design_file> -D <data_dir> -I <number> -init
 
$ cat ${PIPEDIR}/lib/script_templates/feat2.template 
#!/bin/bash
#SBATCH -J feat2-<subject>
#SBATCH --time=00:20:00
#SBATCH --mail-type=FAIL,STAGE_OUT
#SBATCH --mail-user=<mailto>
#SBATCH -o <out_dir>/feat2-slurm-%j
${FSLDIR}/bin/feat <design_file> -D <data_dir> -I <number> -prestats
 
$ cat ${PIPEDIR}/lib/script_templates/feat4.template 
#!/bin/bash
#SBATCH -J feat4-<subject>
#SBATCH --time=00:30:00
#SBATCH --mail-type=FAIL,STAGE_OUT
#SBATCH --mail-user=<mailto>
#SBATCH -o <out_dir>/feat4-slurm-%j
${FSLDIR}/bin/melodic -i <data_dir>/filtered_func_data -o <data_dir>/filtered_func_data.ica -v --nobet --bgthreshold=3 --tr=<TR> --report -d 0 --mmthresh=0.5 --Ostats

y entonces los relleno y los lanzo on the fly,

		foreach my $feat (sort @feats){
			my $template = $pipe_dir.'/lib/script_templates/'.$feat.'.template';
			my $script = $scripts_dir.'/'.$feat.'.sh';
			open TPF, "<$template";
			open SF, ">$script";
			my $name = $fmris{$subject}.$subject;
			while(<TPF>){
				s/<subject>/$name/;
				s/<mailto>/$ENV{'USER'}/;
				s/<out_dir>/$output_dir/;
				s/<design_file>/$dsg_file/;
				s/<data_dir>/$data_dir/g;
				s/<number>/1/;
                                s/<TR>/$test_data{'tr'}/;
				print SF;
			}
			close SF;
			close TPF;
			unless($jobid){
				my $order = 'sbatch '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}else{
				my $order = 'sbatch --depend=afterok:'.$jobid.' '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}
		}
		$jobid="";

Y con esto bastaria pero FSL es muy puñetero con los paths, así que hay que hacer algunos ajustes, copiar archivos de un lado a otro, hacer directorios que se supone que haga feat pero no los hace, etc.

En fin, aqui el codigo completo

rs_ica_one.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_ica_one.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/rs_ica_one.hlp'; exit;}
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
my $fsl = $ENV{'FSLDIR'};
 
my $w_dir=$std{'WORKING'};
my $fmri_dir=$std{'fMRI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
my $template_file = $pipe_dir.'/lib/fsf/lone_ica_template.fsf';
 
print "Collecting needed files\n";
my %fmris;
my %plist = populate($db);
if ($cfile){
	my %cuts = get_pair($data_dir."/".$cfile);
	foreach my $cut (keys %cuts){
		if(grep {/$cut/} %plist){
			$fmris{$cut} = $plist{$cut};
		}
	}
}else{
	%fmris = %plist;
}
my @feats = ("feat1", "feat2", "feat4");
check_or_make($w_dir.'/.files');
system("cp $fsl'/doc/fsl.css' $w_dir'/.files'");
system("cp -r $fsl'/doc/images' $w_dir'/.files/images'");
chdir($w_dir);
foreach my $subject (sort keys %fmris){
	my @names = find(file => 'name' => $fmris{$subject}.$subject."s*.nii.gz", in => $fmri_dir);
	if(@names){
		my $output_dir = $w_dir.'/'.$fmris{$subject}.$subject.'_rsout';
		check_or_make($output_dir);
		my $data = $w_dir.'/'.$fmris{$subject}.$subject.'_rs';
		system("$fsl'/bin/imcp' $names[0] $data");		 
		my $test_order = $fsl.'/bin/fslinfo '.$data;
		open(my $test_info, '-|', $test_order ) or die $!;
		my %test_data;
		while (<$test_info>){
		        if (/^dim1\s*(\d{1,3}).*/) {$test_data{'dim1'}=$1;}
		        if (/^dim2\s*(\d{1,3}).*/) {$test_data{'dim2'}=$1;}
		        if (/^dim3\s*(\d{1,3}).*/) {$test_data{'dim3'}=$1;}
		        if (/^dim4\s*(\d{1,3}).*/) {$test_data{'dim4'}=$1;}
		        if (/^pixdim4\s*(\d{1,2}\.\d{1,6})/) {$test_data{'tr'}=$1;}
		}		
		my $data_size = $test_data{'dim1'}*$test_data{'dim2'}*$test_data{'dim3'}*$test_data{'dim4'};
		my $data_dir = $data.'.ica';
		my $logs_dir = $data_dir.'/logs';
		my $scripts_dir = $data_dir.'/scripts';
		check_or_make($data_dir);
		check_or_make($logs_dir);
		check_or_make($scripts_dir);
		my $dsg_file = $w_dir.'/'.$fmris{$subject}.$subject.'_design.fsf';
		open TPF,"<$template_file";
		open DSG,">$dsg_file";
		while(<TPF>){
			s/<output_dir>/\"$output_dir\"/;
			s/<data_input>/\"$data\"/;
			s/<TR>/$test_data{'tr'}/;
		        s/<number_of_slices>/$test_data{'dim4'}/;
        		s/<size_of_image>/$data_size/;	
			print DSG;
		}
		close DSG;
		close TPF;
		my $jobid;
		foreach my $feat (sort @feats){
			my $template = $pipe_dir.'/lib/script_templates/'.$feat.'.template';
			my $script = $scripts_dir.'/'.$feat.'.sh';
			open TPF, "<$template";
			open SF, ">$script";
			my $name = $fmris{$subject}.$subject;
			while(<TPF>){
				s/<subject>/$name/;
				s/<mailto>/$ENV{'USER'}/;
				s/<out_dir>/$output_dir/;
				s/<design_file>/$dsg_file/;
				s/<data_dir>/$data_dir/g;
				s/<number>/1/;
				s/<TR>/$test_data{'tr'}/;
				print SF;
			}
			close SF;
			close TPF;
			unless($jobid){
				my $order = 'sbatch '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}else{
				my $order = 'sbatch --depend=afterok:'.$jobid.' '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}
		}
		$jobid="";
	}
}

Para lanzar esto basta comprobar que los parametros de la plantilla esten OK (Se han escogido para el proyecto MOPEAD) y lanzar,

$ rs_ica_one.pl mopead

Esto deja un informe (que puede abrirse en un navegador) para cada sujeto con los resultados de MELODIC en un sitio como,

/home/osotolongo/data/mopead/working/mop0001_rs.ica/filtered_func_data.ica/report/00index.html

One ring to rule them all (Group ICA)

Para hacer un ICA de un grupo de sujetos necesitamos cambiar un poco las cosas. Primero, habra dos archivos design.fsf, uno para los registros individuales y otro para correr MELODIC en el grupo. Los que perseguimos es correr los dos primero scripts para cada sujeto y luego el ultimo para todos juntos. Luego este ultimo dependera de que todos los anteriores hayan salido OK. Si sale alguno mal se ha perdido el tiempo, hay que quitar el sujeto del grupo y volver a correr desde el principio. ~~No he encontrado todavia la forma de arreglar esto pero al menos se envia un email por cada sujeto que falla por lo que seria sencillo quitarlos de la lista.~~

La dificultad en rellenar los design.fsf individuales consiste en que hay que llenarlos con la info de todos los sujetos. Asi que primero voy a contarlos y capturar la info,

print "Counting available subjects\n";
foreach my $subject (sort keys %fmris){
	my @names = find(file => 'name' => $fmris{$subject}.$subject."s*.nii.gz", in => $fmri_dir);
	if(@names){
		my $sname = $fmris{$subject}.$subject;
		$subjects{$sname} = $names[0]; 
		$count++;
	}
}
my $pollos = $count;

Ahora voy a tomar el primer sujeto de la lista y voy obtener los datos para rellenar la plantilla.

my $test_subject = ( sort keys %subjects )[0];                                                                        
my $test_path = $subjects{$test_subject};                                                                             
my $test_order = $fsl.'/bin/fslinfo '.$test_path;                                                                     
open(my $test_info, '-|', $test_order ) or die $!;                                                                    
my %test_data;                                                                                                        
while (<$test_info>){                                                                                                 
        if (/^dim1\s*(\d{1,3}).*/) {$test_data{'dim1'}=$1;}                                                           
        if (/^dim2\s*(\d{1,3}).*/) {$test_data{'dim2'}=$1;}                                                           
        if (/^dim3\s*(\d{1,3}).*/) {$test_data{'dim3'}=$1;}                                                           
        if (/^dim4\s*(\d{1,3}).*/) {$test_data{'dim4'}=$1;}                                                           
        if (/^pixdim4\s*(\d{1,2}\.\d{1,6})/) {$test_data{'tr'}=$1;}                                                   
}                                                                                                                     
my $data_size = $test_data{'dim1'}*$test_data{'dim2'}*$test_data{'dim3'}*$test_data{'dim4'};

Este mismo procedimiento lo sigo para todos los demas sujetos y comparo. Esto me permite eliminar del procesamiento los sujetos cuyas imagenes tengan problemas y vayan a darme error al lanzar los scripts. De esta manera, no es necesario preocuparse por relanzar nuevamente en caso de fallo. Los problemas que tengalos guardo en el archivo fmri_image_errors.txt asi que pueden revisarse despues.

print "Checking images and excluding wrong subjects\n";                                                               
open(EDF, ">$error_data_file");                                                                                       
foreach my $subject (sort keys %subjects){                                                                            
        my $sub_path = $subjects{$subject};                                                                           
        my $order = $fsl.'/bin/fslinfo '.$sub_path;                                                                   
        open(my $sub_info, '-|', $order ) or die $!;                                                                  
        my %sub_data;                                                                                                 
        while (<$sub_info>){                                                                                          
                if (/^dim1\s*(\d{1,3}).*/) {$sub_data{'dim1'}=$1;}                                                    
                if (/^dim2\s*(\d{1,3}).*/) {$sub_data{'dim2'}=$1;}                                                    
                if (/^dim3\s*(\d{1,3}).*/) {$sub_data{'dim3'}=$1;}                                                    
                if (/^dim4\s*(\d{1,3}).*/) {$sub_data{'dim4'}=$1;}                                                    
                if (/^pixdim4\s*(\d{1,2}\.\d{1,6})/) {$sub_data{'tr'}=$1;}                                            
        }                                                                                                             
        unless(($sub_data{'dim1'} == $test_data{'dim1'}) && ($sub_data{'dim2'} == $test_data{'dim2'}) && ($sub_data{'dim3'} == $test_data{'dim3'}) && ($sub_data{'dim4'} == $test_data{'dim4'}) && ($sub_data{'tr'} == $test_data{'tr'})){  
                print EDF "\n$subject\ndim1: $sub_data{'dim1'}\ndim2: $sub_data{'dim2'}\ndim3: $sub_data{'dim3'}\ndim4: $sub_data{'dim4'}\ntr: $sub_data{'tr'}\n";                                                                          
                delete($subjects{$subject});                                                                          
        }                                                                                                             
}                                                                                                                     
close EDF;

Ahora bien, puede ocurrir que la imagen que tome de referencia sea la que tiene problemas. Esto me dejara con una lista deun solo sujeto. Asi que compruebo el tamaño d ela lista y decido si voy a proseguir o no.

my $subs_size = keys %subjects;                                                                                       
unless ($subs_size > 1){                                                                                              
        my ($bad) = %subjects;                                                                                        
        print "Error: bad subject taken as reference or individual group\n";                                          
        print "Delete subject $bad from list and try again\n";                                                        
        exit;                                                                                                         
}

A esta altura, todas las imagenes han sido comprobadas y tengo todos los datos necesarios. Voy a rellenar plantillas.

Las plantillas de los design.fsf son distintas ahora,

my $lone_template_file = $pipe_dir.'/lib/notalone_ica_template.fsf';
my %template_files = ( "begin" => $pipe_dir.'/lib/group_ica_template_p1.fsf', 
					"data" => $pipe_dir.'/lib/group_ica_template_targets.fsf', 
					"end" => $pipe_dir.'/lib/group_ica_template_p2.fsf');

La del grupo la he dividido en tres partes, la primera donde las variables importantes son,

# Output directory
set fmri(outputdir) <output_dir>

# TR(s)
set fmri(tr) <TR>

# Total volumes
set fmri(npts) <number_of_slices>

# Number of first-level analyses
set fmri(multiple) <number_of_subjects>

# Total voxels
set fmri(totalVoxels) <size_of_image>

La segunda que consiste en llenar la info para cada individuo, (todos aqui, uno a uno)

# 4D AVW data or FEAT directory (<number_of_sample>)
set feat_files(<number_of_sample>) <data>

y la tercera que es solo un bloque de texto fijo. Esto se hace mas o menos asi,

my $dsg_file = $w_dir.'/gica_design.fsf';
open DSG,">$dsg_file";
open TPF,"<$template_files{'begin'}";
 
print "Making global .fsf file\n";
while(<TPF>){
	s/<output_dir>/\"$output_dir\"/;
	s/<number_of_subjects>/$pollos/;
        s/<number_of_slices>/$test_data{'dim4'}/;
        s/<size_of_image>/$data_size/; 
        s/<TR>/$test_data{'tr'}/;
	print DSG;
}
close TPF;
$count = 1;
my $filtered_list = $gdata_dir.'/.filelist'; # que bonito!
open FCK, ">$filtered_list";
foreach my $subject (sort keys %subjects){
		my $idata = $w_dir.'/'.$subject.'_rs.ica';
		my $fck_line = $idata.'/reg_standard/filtered_func_data';
		print FCK "$fck_line\n";
		open TPF,"<$template_files{'data'}";
		while(<TPF>){
			s/<data>/$idata/;
			s/<number_of_sample>/$count/;
			print DSG;
		}
		$count++;
		close TPF;
}
close FCK;
open TPF,"<$template_files{'end'}";
while(<TPF>){
	print DSG;
}
close TPF;
close DSG;

Nota: un detalle aqui es que hay que ir haciendo la lista de outputs de los registros individuales y metiendolos en el archivo .filelist que hara falta luego para lanzar MELODIC.

Y ahora vamos a colectar la info y para ponerla en los archivos individuales,

print "Making individual .fsf files and scripts\n";
my $tranca="";
$count = 1;
foreach my $subject (sort keys %subjects){
                my $idata = $w_dir.'/'.$subject.'_rs';
                open TPF,"<$template_files{'data'}";
                while(<TPF>){
                        s/<data>/$idata/;
                        s/<number_of_sample>/$count/;
                        $tranca .= $_;
                }
                $count++;
                close TPF;
}

No hemos escrito nada aun pero ya tenemos el bloque que hay que escribir en cada caso.

Nota: Ojo con los paths. Aqui los input son _rs (que son los nifti) y en el global son _rs.ica (que son los directorios de salid de los dos primeros scripts).

Ahora, terminamos de llenar, para cada sujeto, sus archivos design.fsf,

		my $ioutput_dir = $w_dir.'/'.$subject.'_rsout';
		check_or_make($ioutput_dir);
		my $idata = $w_dir.'/'.$subject.'_rs';
		my $idata_dir = $idata.'.ica';
		my $ilogs_dir = $idata_dir.'/logs';
		my $iscripts_dir = $idata_dir.'/scripts';
		check_or_make($idata_dir);
		check_or_make($ilogs_dir);
		check_or_make($iscripts_dir);
		my $idsg_file = $w_dir.'/'.$subject.'_rs.ica/design.fsf';
		open TPF,"<$lone_template_file";
		open DSG,">$idsg_file";
		while(<TPF>){
			s/<output_dir>/\"$ioutput_dir\"/;
			s/<full_data>/$tranca/;
                        s/<TR>/$test_data{'tr'}/;
                        s/<number_of_slices>/$test_data{'dim4'}/;
                        s/<size_of_image>/$data_size/;
                        s/<number_of_subjects>/$pollos/;
			print DSG;
		}
		close DSG;
		close TPF;		
		system("$fsl'/bin/imcp' $subjects{$subject} $idata");

Nota: Si no se pone este archivo exactamente en el directorio de entrada y con ese nombre, MELODIC no funcionara.

Ahora ya lanzamos los scripts, que son iguales (pero solo los dos primeros)

		my $jobid;
		foreach my $feat (sort @feats){
			my $template = $pipe_dir.'/lib/script_templates/'.$feat.'.template';
			my $script = $iscripts_dir.'/'.$feat.'.sh';
			open TPF, "<$template";
			open SF, ">$script";
			while(<TPF>){
				s/<subject>/$subject/;
				s/<mailto>/$ENV{'USER'}/;
				s/<out_dir>/$ioutput_dir/;
				s/<design_file>/$idsg_file/;
				s/<data_dir>/$idata_dir/g;
				s/<number>/$count/;
				print SF;
			}
			close SF;
			close TPF;
			unless($jobid){
				my $order = 'sbatch '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}else{
				my $order = 'sbatch --depend=afterok:'.$jobid.' '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
				push @jobs_list, $jobid;
			}
		}
		$jobid="";
		$count++;

y ya podriamos lanzar el ultimo porque hemos guardado en @jobs_list todos los job_id que necesitamos, pero este script es un poco distinto,

$ cat ${PIPEDIR}/lib/script_templates/feat4_ica.template 
#!/bin/bash
#SBATCH -J feat4-ica
#SBATCH --time=00:30:00
#SBATCH --mail-type=FAIL,STAGE_OUT,END
#SBATCH --mail-user=<mailto>
#SBATCH -o <out_dir>/feat4-slurm-%j
${FSLDIR}/bin/feat <design_file> -D <output_dir> -gica

Es mas sencillo de llenar pero mucho mas complejo si da problemas

print "Making global script\n";
my $logs_dir = $gdata_dir.'/logs';
my $scripts_dir = $gdata_dir.'/scripts';
check_or_make($logs_dir);
check_or_make($scripts_dir);
 
my $template = $pipe_dir.'/lib/script_templates/'.$gfeat.'.template';
my $script = $scripts_dir.'/'.$gfeat.'.sh';
open TPF, "<$template";
open SF, ">$script";
while(<TPF>){
	s/<mailto>/$ENV{'USER'}/;
	s/<out_dir>/$output_dir/;
	s/<design_file>/$dsg_file/;
	s/<output_dir>/$gdata_dir/g;
	print SF;
}
close SF;
close TPF;
my $sjobs = join(',',@jobs_list);
my $order = 'sbatch --depend=afterok:'.$sjobs.' '.$script;
print "$order\n";
exec($order);

todo junto aqui

rs_ica_group.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Slurp qw(read_file);
use File::Find::Rule;
use File::Basename qw(basename);
use Data::Dump qw(dump);
use File::Copy::Recursive qw(dirmove);
use File::Copy qw(copy);
 
use NEURO qw(populate check_subj load_study print_help get_pair check_or_make);
 
my $study;
my $cfile;
 
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/rs_ica_one.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/rs_ica_one.hlp'; exit;}
my %std = load_study($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
my $fsl = $ENV{'FSLDIR'};
 
my $w_dir=$std{'WORKING'};
my $fmri_dir=$std{'fMRI'};
my $mri_dir=$std{'MRI'};
my $db = $std{'DATABASE'};
my $data_dir=$std{'DATA'};
my $lone_template_file = $pipe_dir.'/lib/fsf/notalone_ica_template.fsf';
my %template_files = ( "begin" => $pipe_dir.'/lib/fsf/group_ica_template_p1.fsf', 
					"data" => $pipe_dir.'/lib/fsf/group_ica_template_targets.fsf', 
					"end" => $pipe_dir.'/lib/fsf/group_ica_template_p2.fsf');
 
 
print "Collecting needed files\n";
my %fmris;
my %plist = populate($db);
if ($cfile){
	my %cuts = get_pair($data_dir."/".$cfile);
	foreach my $cut (keys %cuts){
		if(grep {/$cut/} %plist){
			$fmris{$cut} = $plist{$cut};
		}
	}
}else{
	%fmris = %plist;
}
my $count = 0;
my %subjects;
print "Counting available subjects\n";
foreach my $subject (sort keys %fmris){
	my @names = find(file => 'name' => $fmris{$subject}.$subject."s*.nii.gz", in => $fmri_dir);
	if(@names){
		my $sname = $fmris{$subject}.$subject;
		$subjects{$sname} = $names[0]; 
		$count++;
	}
}
my $pollos = $count;
print "Getting info from images\n";
my $test_subject = ( sort keys %subjects )[0];
my $test_path = $subjects{$test_subject};
my $test_order = $fsl.'/bin/fslinfo '.$test_path;
open(my $test_info, '-|', $test_order ) or die $!;
my %test_data;
while (<$test_info>){
	if (/^dim1\s*(\d{1,3}).*/) {$test_data{'dim1'}=$1;}
	if (/^dim2\s*(\d{1,3}).*/) {$test_data{'dim2'}=$1;}
	if (/^dim3\s*(\d{1,3}).*/) {$test_data{'dim3'}=$1;}
	if (/^dim4\s*(\d{1,3}).*/) {$test_data{'dim4'}=$1;}
	if (/^pixdim4\s*(\d{1,2}\.\d{1,6})/) {$test_data{'tr'}=$1;}
}
my $data_size = $test_data{'dim1'}*$test_data{'dim2'}*$test_data{'dim3'}*$test_data{'dim4'};
print "Copying FSL files and setting directories\n";
my @feats = ("feat1", "feat2");
my $gfeat = "feat4_ica";
check_or_make($w_dir.'/.files');
system("cp $fsl'/doc/fsl.css' $w_dir'/.files'");
system("cp -r $fsl'/doc/images' $w_dir'/.files/images'");
chdir($w_dir);
my @jobs_list;
my $gdata_dir = $w_dir.'/rs.gica';
my $output_dir = $w_dir.'/rsout_gica';
check_or_make($output_dir);
check_or_make($gdata_dir);
my $dsg_file = $w_dir.'/gica_design.fsf';
open DSG,">$dsg_file";
open TPF,"<$template_files{'begin'}";
 
print "Making global .fsf file\n";
while(<TPF>){
	s/<output_dir>/\"$output_dir\"/;
	s/<number_of_subjects>/$pollos/;
	s/<number_of_slices>/$test_data{'dim4'}/;
	s/<size_of_image>/$data_size/; 
	s/<TR>/$test_data{'tr'}/;
	print DSG;
}
close TPF;
$count = 1;
my $filtered_list = $gdata_dir.'/.filelist';
open FCK, ">$filtered_list";
foreach my $subject (sort keys %subjects){
		my $idata = $w_dir.'/'.$subject.'_rs.ica';
		my $fck_line = $idata.'/reg_standard/filtered_func_data';
		print FCK "$fck_line\n";
		open TPF,"<$template_files{'data'}";
		while(<TPF>){
			s/<data>/$idata/;
			s/<number_of_sample>/$count/;
			print DSG;
		}
		$count++;
		close TPF;
}
close FCK;
open TPF,"<$template_files{'end'}";
while(<TPF>){
	print DSG;
}
close TPF;
close DSG;
 
print "Making individual .fsf files and scripts\n";
my $tranca="";
$count = 1;
foreach my $subject (sort keys %subjects){
		my $idata = $w_dir.'/'.$subject.'_rs';
		open TPF,"<$template_files{'data'}";
		while(<TPF>){
			s/<data>/$idata/;
			s/<number_of_sample>/$count/;
			$tranca .= $_;
		}
		$count++;
		close TPF;
}
 
$count = 1;
foreach my $subject (sort keys %subjects){
		my $ioutput_dir = $w_dir.'/'.$subject.'_rsout';
		check_or_make($ioutput_dir);
		my $idata = $w_dir.'/'.$subject.'_rs';
		my $idata_dir = $idata.'.ica';
		my $ilogs_dir = $idata_dir.'/logs';
		my $iscripts_dir = $idata_dir.'/scripts';
		check_or_make($idata_dir);
		check_or_make($ilogs_dir);
		check_or_make($iscripts_dir);
		my $idsg_file = $w_dir.'/'.$subject.'_rs.ica/design.fsf';
		open TPF,"<$lone_template_file";
		open DSG,">$idsg_file";
		while(<TPF>){
			s/<output_dir>/\"$ioutput_dir\"/;
			s/<full_data>/$tranca/;
			s/<TR>/$test_data{'tr'}/;
			s/<number_of_slices>/$test_data{'dim4'}/;
			s/<size_of_image>/$data_size/;
			s/<number_of_subjects>/$pollos/;
			print DSG;
		}
		close DSG;
		close TPF;		
		system("$fsl'/bin/imcp' $subjects{$subject} $idata"); 
 
		my $jobid;
		foreach my $feat (sort @feats){
			my $template = $pipe_dir.'/lib/script_templates/'.$feat.'.template';
			my $script = $iscripts_dir.'/'.$feat.'.sh';
			open TPF, "<$template";
			open SF, ">$script";
			while(<TPF>){
				s/<subject>/$subject/;
				s/<mailto>/$ENV{'USER'}/;
				s/<out_dir>/$ioutput_dir/;
				s/<design_file>/$idsg_file/;
				s/<data_dir>/$idata_dir/g;
				s/<number>/$count/;
				s/<TR>/$test_data{'tr'}/;
				print SF;
			}
			close SF;
			close TPF;
			unless($jobid){
				my $order = 'sbatch '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
			}else{
				my $order = 'sbatch --depend=afterok:'.$jobid.' '.$script;
				print "$order\n";
				$jobid = `$order`;
				$jobid = ( split ' ', $jobid )[ -1 ];
				push @jobs_list, $jobid;
			}
		}
		$jobid="";
		$count++;
}
 
print "Making global script\n";
my $logs_dir = $gdata_dir.'/logs';
my $scripts_dir = $gdata_dir.'/scripts';
check_or_make($logs_dir);
check_or_make($scripts_dir);
 
my $template = $pipe_dir.'/lib/script_templates/'.$gfeat.'.template';
my $script = $scripts_dir.'/'.$gfeat.'.sh';
open TPF, "<$template";
open SF, ">$script";
while(<TPF>){
	s/<mailto>/$ENV{'USER'}/;
	s/<out_dir>/$output_dir/;
	s/<design_file>/$dsg_file/;
	s/<output_dir>/$gdata_dir/g;
	print SF;
}
close SF;
close TPF;
my $sjobs = join(',',@jobs_list);
my $order = 'sbatch --depend=afterok:'.$sjobs.' '.$script;
print "$order\n";
exec($order);

y vamos a probarlo,

$ cat test_fmri.csv 
0001;mop
0002;mop
0003;mop
0004;mop
 
$ rs_ica_group.pl -cut test_fmri.csv mopead
Collecting needed files
Counting available subjects
Copying FSL files and setting directories
Making global .fsf file
Making individual .fsf files and scripts
sbatch /nas/osotolongo/data/mopead/working/mop0001_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:2172 /nas/osotolongo/data/mopead/working/mop0001_rs.ica/scripts/feat2.sh
sbatch /nas/osotolongo/data/mopead/working/mop0002_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:2174 /nas/osotolongo/data/mopead/working/mop0002_rs.ica/scripts/feat2.sh
sbatch /nas/osotolongo/data/mopead/working/mop0003_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:2176 /nas/osotolongo/data/mopead/working/mop0003_rs.ica/scripts/feat2.sh
sbatch /nas/osotolongo/data/mopead/working/mop0004_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:2178 /nas/osotolongo/data/mopead/working/mop0004_rs.ica/scripts/feat2.sh
Making global script
sbatch --depend=afterok:2173,2175,2177,2179 /nas/osotolongo/data/mopead/working/rs.gica/scripts/feat4_ica.sh
Submitted batch job 2180
 
$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
              2173     debug feat2-mo osotolon PD       0:00      1 (Dependency)
              2175     debug feat2-mo osotolon PD       0:00      1 (Dependency)
              2177     debug feat2-mo osotolon PD       0:00      1 (Dependency)
              2179     debug feat2-mo osotolon PD       0:00      1 (Dependency)
              2180     debug feat4-ic osotolon PD       0:00      1 (Dependency)
              2172     debug feat1-mo osotolon  R       0:17      1 brick01
              2174     debug feat1-mo osotolon  R       0:14      1 brick02
              2176     debug feat1-mo osotolon  R       0:11      1 brick03
              2178     debug feat1-mo osotolon  R       0:11      1 detritus

Cuando termina manda un email,

y deja un informe en working/rs.gica/groupmelodic.ica/report/00index.html.

FSLNets

FSLNets es un conjunto de scripts escritos en MATLAB con los cuales se puede extraer la red de un grupo de imagenes fMRI.

Se han de realizar algunos pasos previos.

El tratamiento con melodic ya se ha hecho en el ICA grupal. Los archivos estan guardados en working/rs.gica/.filelist.

[osotolongo@detritus mopead]$ cat working/rs.gica/.filelist
/nas/data/mopead/working/mop0001_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0002_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0003_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0004_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0005_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0006_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0007_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0008_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0009_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0010_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0012_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0014_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0015_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0017_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0018_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0019_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0020_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0021_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0022_rs.ica/reg_standard/filtered_func_data
/nas/data/mopead/working/mop0023_rs.ica/reg_standard/filtered_func_data

La sintaxis de dual_regression es incorrecta. La nueva version usa otro orden para los parametros. ver ayuda!

***NOTE*** ORDER OF COMMAND-LINE ARGUMENTS IS DIFFERENT FROM PREVIOUS VERSION

Usage: dual_regression <group_IC_maps> <des_norm> <design.mat> <design.con> <n_perm> <output_directory> <input1> <input2> <input3> .........
e.g.   dual_regression groupICA.gica/groupmelodic.ica/melodic_IC 1 design.mat design.con 500 grot `cat groupICA.gica/.filelist`

<group_IC_maps_4D>            4D image containing spatial IC maps (melodic_IC) from the whole-group ICA analysis
<des_norm>                    0 or 1 (1 is recommended). Whether to variance-normalise the timecourses used as the stage-2 regressors
<design.mat>                  Design matrix for final cross-subject modelling with randomise
<design.con>                  Design contrasts for final cross-subject modelling with randomise
<n_perm>                      Number of permutations for randomise; set to 1 for just raw tstat output, set to 0 to not run randomise at all.
<output_directory>            This directory will be created to hold all output and logfiles
<input1> <input2> ...         List all subjects' preprocessed, standard-space 4D datasets

<design.mat> <design.con>     can be replaced with just
-1                            for group-mean (one-group t-test) modelling.
If you need to add other randomise options then edit the line after "EDIT HERE" in the dual_regression script

Ha de hacerse algo así,

[osotolongo@detritus mopead]$ dual_regression working/rs.gica/groupmelodic.ica/melodic_IC 1 -1 500 working/groupICA100.dr `cat working/rs.gica/.filelist`
creating common mask
doing the dual regressions
......

Es extremadamente lento y no se como paralelizar randomize en slurm. Viene solo preparado para SGE. Debo mirar aqui.

Despues se hace slices_summary,

[osotolongo@detritus mopead]$ slices_summary working/rs.gica/groupmelodic.ica/melodic_IC 4 $FSLDIR/data/standard/MNI152_T1_2mm working/groupICA100.sum -1
Generating single-slice summary images instead of 3-slice
Upper colour threshold is 12.576081 
......

De aqui me voy a ir a brick02 que es donde he instalado octave.

[osotolongo@brick02 ~]$ octave
GNU Octave, version 3.8.2
Copyright (C) 2014 John W. Eaton and others.
..........
octave:1> group_maps='/nas/data/mopead/working/rs.gica/groupmelodic.ica/melodic_IC'
group_maps = /nas/data/mopead/working/rs.gica/groupmelodic.ica/melodic_IC
octave:2> addpath /nas/software/FSLNets
octave:3> addpath /nas/software/L1precision
octave:4> addpath /nas/software/pwling
octave:5> addpath(sprintf('%s/etc/matlab',getenv('FSLDIR')))
octave:6> ts_dir='/nas/data/mopead/working/groupICA100.dr';
octave:7> ts=nets_load(ts_dir,2.13,1);

A ver si ha leido los timeseries (esto es importante)

octave:8> ts
ts =

  scalar structure containing the fields:

    NtimepointsPerSubject =  195
    ts =   

     Columns 1 through 14:

       2.0802e-01  -4.0387e-01  -6.1757e-01  -1.5950e-01  -9.1144e-02   8.6173e-01   1.6692e+00   9.0748e-01  -2.2963e-01  -7.0069e-01  -7.3951e-01  -4.2056e-01   2.2683e-01  -9.9498e-01
      -1.6867e-01  -2.7437e-02  -6.0164e-01  -4.5625e-01  -3.7668e-01   1.1882e-01   9.6679e-01   9.1118e-02   7.7420e-01  -7.8032e-01  -6.3792e-01   9.9788e-01  -2.5313e-01  -5.7240e-01
      -1.0760e+00   1.7245e-01  -4.7770e-01  -1.2929e+00  -7.6983e-01  -7.0703e-01  -2.6424e-01  -2.5222e-01   1.8566e+00   2.6005e-01   3.3061e-01   6.4488e-01   2.2098e-01  -2.7346e-01
       5.2492e-01   4.2167e-02  -5.9531e-01   2.8384e-01  -1.4880e-01   3.7651e-01   2.4698e-01   1.9122e-01  -1.8680e-01   8.2724e-03  -3.9874e-01   4.4234e-01  -1.8046e-01   2.7766e-01
       9.5172e-01   2.7839e-01  -1.5722e+00   3.6075e-01   3.5335e-01   2.2732e-01   1.8472e+00   6.2682e-01  -7.0714e-01  -5.7604e-01  -1.3205e+00   1.2185e+00  -1.1469e+00  -8.8856e-01
      -1.3355e+00   2.7909e-01  -6.4247e-01  -1.1779e+00  -1.0891e+00  -5.1818e-01   4.9458e-02  -4.1333e-01   1.4599e+00   5.6708e-02  -2.6263e-01   1.0073e-01  -1.5722e-01  -2.6927e-01
       1.0286e+00   8.2392e-01  -9.2610e-01   6.3058e-01  -1.8016e-01  -1.7455e-01  -4.6052e-01  -2.7851e-01  -1.4289e-01   3.1517e-01  -7.3147e-01   6.0409e-01  -8.3970e-01   2.1175e-02
       5.8707e-01   3.6076e-01  -1.8519e+00  -1.6925e-01  -2.3339e-01   4.6212e-01   1.5051e+00   5.2714e-01  -4.7612e-01  -5.1707e-01  -7.0160e-01  -3.3002e-01  -1.5390e+00  -5.1593e-01
      -1.2313e+00  -1.0890e+00   2.8377e-01  -1.5782e+00  -1.3756e+00   1.2942e+00   1.3311e+00   1.2493e+00   1.1025e-01  -8.3362e-01   8.2314e-01  -2.0298e+00  -6.2277e-01   7.5722e-01
       2.6974e-01  -3.0517e-01   1.0004e+00   1.3926e-03   6.0781e-01   3.1420e-01  -1.5776e-01   1.8985e-01  -2.2337e-01  -1.9407e-01   2.7684e-01  -1.4240e+00  -5.0345e-01   6.1336e-03
      -4.2959e-01  -7.8180e-01   1.04
..................

Y de aqui, el espectro,

octave:9> ts_spectra=nets_spectra(ts);
+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|                                                                                                                                                                                                |
| 20 ++-----+-----+------+-----+-----++  35 ++-----+-----+------+-----+-----++  50 ++----+------+-----+------+----++    1 ++-----------+------------+------------+------------+-----------++     |
|    +-------------------------------++     +-------------------------------++     +------------------------------++      +---------------------------------------------------------------++     |
|    |                               ||     |                               ||     |                              ||      |  ||||| +                                                      ||     |
|    |                               ||     |                               ||     |                              ||      |  |||+| ++                                                     ||     |
|    |                               ||     |                               ||     |                              ||      |  |+||| ||                                                     ||     |
|    |                               ||     |                               ||     |                              ||      | ||+||+ ||+                                                    ||     |
|    |                               ||     |                               ||     |                              ||      | |+|||||||+ +                                                 +||     |
|    |                               ||     |                               ||     |                              ||      | |||+|+||||++                                          +    ++|||     |
|    |+++++++ +                      ||     |+++++++++                     +||     |+++++++++ +        + + +++ +++||  0.8 |+|+||++|+|+||                                          |    |+||+     |
| 15 |+++++++++++++++++++++++++++++++|+  30 |+++++++++++++++++++++++++++++++|+  45 |++++++++++++++++++++++++++++++|+      | ||+||+||||||                                       +  |   +|||||     |
|    |++++++++++  ++                +||     |++++++++++                    +||     |+++++++++++                  +||      | ||++||||+++|| +                                ++  |  ||  ||+|||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      | ||||+|++|||+| + +    +            +       +    ||+ + | |+ +|||||     |
|    |++++++++++++++   +++++  +++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      | +||| |||||+|| | |    +   +        |   ++  |    |+|+|||++| |||+||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      | |||+ ||+|+|+| + + +++| + |        || +|++ |++  ||+|+|++||++|+|||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||||++ ||||||++++|| +|++|| | + +  + |+ ++|| +||+ ||+||+|+||+||||||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||  0.6 |+|||+ |+|++||||||+ |||||+|| + |  | +|||+|| |+|||+|||||++||||||+|+     |
|    |++ ++++++++++++++++++++++++++++||     |++ ++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||+|+  +||||+|++|+||++||||+++||| +| ||+||++ ||||+||||||||+|||++|||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      |||||  || |||+|||+||||||||++||++ ||||+||+++|+||||++++||||++|+||+||     |
|    |+++++++++++++++++++++++++++++++||     |++ ++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||||   |+ |+||||+|+||+++|||+|+++++||++|+++|+++|++||+|+|+||+|||||||     |
| 10 |+++++++++++++++++++++++++++++++|+  25 |+++++++++++++++++++++++++++++++|+  40 |++++++++++++++++++++++++++++++|+      ||||   +||++|++++|||||+++++|+|+|+|+||+|++|++|+|+||+|+++|+||+|+||||     |
|    |++ ++++++++++++++++++++++++++++||     |++ ++++++++++++++++++++++++++++||     |++ +++++++++++++++++++++++++++||      ||||    ++++||++||||||+|+||||||+||++||++|||||+|||+|+|++|++||||||||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      |+||      + |+||+||||||||+||+|++||+|+||||+|++|+|+++++||+||+|+|+|||     |
|    |+  ++++++++++++++++++++++++++++||     |++ ++++++++++++++++++++++++++++||     |++ +++++++++++++++++++++++++++||      ||+         ||++|+|+++|+|+++++||||||++++++||+|+|+++++++||||+|+|+||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||  0.4 |+|         +|++||+||||++|++++++++||+||+++|||||+|++|+|||+++||||||+     |
|    |+  ++++++++++++++++++++++++++++||     |+  ++++++++++++++++++++++++++++||     |+  +++++++++++++++++++++++++++||      |||         |||++|||+++|||+||++|++|||++|+||+|+|++|++|+||++++||++||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      |||         | ++|+|||+++||||+|+|||++|+||++||++++||++||++|||+||++||     |
|    |+   +++++++++++++++++++++++++++||     |+  + ++++++++++++++++++++++++++||     |+  + +++++++++++++++++++++++++||      |+|         |   + ++|++|+++|+++++++|++++++++|+||++|+++|++||++|||||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      |||         |   +   +++++++++++||+++|+|+|+++++|+++||+++++++++|+ ||     |
|  5 |+++++++++++++++++++++++++++++++|+  20 |+++++++++++++++++++++++++++++++|+  35 |++++++++++++++++++++++++++++++|+      ||          +       +  + + ++++||+++| +++ ++|++++|+++++|+++||++|||     |
|    |+++++++++                     +||     |++++++++++    +  + +    ++ ++++||     |++++++++++++ +               +||      ||                         | ++++++++ + ++||++++|+|++|||| ++++ +||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||  0.2 |+                         +++   ||+ ++|+ ++++++++ | |++++|++|  |+     |
|    |+++++++++++++++++++++++++++++++||     |++++++++++++++   +++++  +++++++||     |+++++++++++++    ++++  +++++++||      ||                               ++++ ++   +  +++  +++   +  +   ||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||                                +                             ||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||                                                              ||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||                                                              ||     |
|    |+++++++++++++++++++++++++++++++||     |++ ++++++++++++++++++++++++++++||     |++ +++++++++++++++++++++++++++||      ||                                                              ||     |
|    |+++++++++++++++++++++++++++++++||     |+++++++++++++++++++++++++++++++||     |++++++++++++++++++++++++++++++||      ||                                                              ||     |
|    |      +     +      +     +     |+     |      +     +      +     +     |+     |     +      +     +      +    |+      ||           +            +            +            +           |+     |
|  0 +-------------------------------++  15 +-------------------------------++  30 +------------------------------++    0 +---------------------------------------------------------------++     |
+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

Detritus Wiki

Table of Contents

MOPEAD

Conversion DICOM (full)

MRI

Pruebas iniciales

DB y preprocesamiento T1

Preprocesamiento T2

Paralelizando Freesurfer

Una nota sobre el cluster

Ahora si, Freesurfer

DTI

Entendiendo la estructura

Escogiendo el DTI

Preprocesamiento

Corregistro

Procesamiento en paralelo

Backward compatibility

Configurar PATHs

Haciendo máscaras

Tractografía (How to)

bedpostx

probtrackx

Tractografía (Allons-y!)

Esto tiene un problema logistico

Testing

Tractografía (Otros atlas?)

Tractografía (tojunto)

Single mask

Tractografía (Seed and Targets)

Alt metrics

trick or treat

fMRI

Mirando los DICOMs

seleccion de series fMRI a partir de los NIfTIs

Independent Component Analysis (AKA ICA)

One ring to rule them all (Group ICA)

FSLNets