Table of Contents
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.
El siguiente paso es convertir todo a nifti. para ello voy a hacer un directorio en $std{'DATA'}./niftis/ y convierto ahi todo.
Ahora toca buscar y extraer las T1, leer la DB y copiarlas a mri/ .
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.
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; } }
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
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
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.
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.
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.
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
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
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
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;
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)
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.
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,
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.
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
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,
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.
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.
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 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
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);
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 +---------------------------------------------------------------++ | +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+