Pipeline de procesamiento de neuroimagen de ACE (acenip v0.6)

Porqué la nueva versión

Básicamente, estamos enfocando todo hacia XNAT.

Gran parte del procesamiento se ha desplazado a XNAT. Los procedimientos para Freesurfer, pgsWHM y PET-FBB son obsoletos.
El almacenamiento de resultados debe hacerse en XNAT. Los archivos intermedios deben poder borrarse. Los informes finales de procesamiento deben salir de XNAT.
Luego, el procesamiento no basico como DTI, fMRI o cosas mas complejas como el el analisis de hipocampo o los SBM si que han de hacerse fuera de XNAT. Hay que describir bien como organizar una estructura BIDS a partir de las imagenes guardadas en XNAT y documentar el procesamiento posterior.
Todo el procedimiento inicial alredededor de crear la base de datos del proyecto y el posterior procesamiento de resultados es obsoleto. La DB de trabajo debe generarse dinamicamente a partir de XNAT y puede ser para el proyecto completo o para un procesamiento parcial. Asi, el procesamiento es independiente por completo de la estructura de proyecto.

Una parte muy importante de este esfuerzo es la biblioteca de interaccion con la API de XNAT que permite construir facilmente scripts de interaccion con los proyectos solo ensamblando los bloques necesarios.

Dependencias

XNAT: http://www.xnat.org/
dcm2bids: https://pypi.org/project/dcm2bids/
FSL: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/
Freesurfer: https://surfer.nmr.mgh.harvard.edu/
ANTs: http://stnava.github.io/ANTs/
Convert3D: http://www.itksnap.org/pmwiki/pmwiki.php?n=Downloads.C3D
fmriprep: https://fmriprep.org/en/stable/
dicom3tools: https://www.dclunie.com/dicom3tools.html (no estrictamente necesario pero conveniente)

Flujo planeado

TO DO

Reescribir toda la documentacion para enfocar el flujo de trabajo basico hacia XNAT
Dar la posibilidad de Crear projectos locales al vuelo. Esto es, el proyecto local no necesita coincidir con un proyecto de XNAT, puede ser un subconjunto de este o incluso una interseccion entre varios proyectos de XNAT. (Esto ultimo puede ser complicado!)
Limpiar la libreria de interaccion con la API de XNAT, dejando solo las funciones necesarias
El codigo necestia ser mantenible. Hay reescribir los distintos procedimientos a un estilo unico.
XNAT to BIDS: Bajar directamente los archivos DICOM necesarios, convertir a estructura BIDS y borrar archivos intermedios

interactuando con XNAT

interfaz WEB

La interfaz web de XNAT en ACE esta en la direccion http://xnat.fundacioace.com:8088/. En principio el frontend se intenta mantener medianamente actualizado y toda la documentacion que puede encontrarse en https://wiki.xnat.org/ es aplicable.

En particular, para obtener resultados agilmente, es interesante la API: https://wiki.xnat.org/display/XAPI

xnatapic

Aunque el funcionamiento del pipeline es completamente independiente de xnatapic (https://github.com/asqwerty666/xnatapic), si que es una herramienta interesante para hace consultas rapidas a XNAT y crear reportrs independientes.

Aunque tambien podemos obtener muchos de estos reports desde la interfaz web de XNAT, es mucho mas rapido y sencillo usar xnatapic la mayoria de las veces. Por ejemplo, para obtener el listado de sujetos en un proyecto podemos ir a la interfaz web, buscar el proyecto y bajar un CSV con la lista de sujetos. O simplemente escribir,

[osotolongo@brick03 mri_face]$ xnatapic list_subjects --project_id mriface --label | awk -F',' '{print $2}'> mriface_sujetos.list
[osotolongo@brick03 mri_face]$ head mriface_sujetos.list
20211480
20211611
20210716
20211438
20211637
20211281
20200484
20211401
20211683
20211595

Esto nos permite ademas interacciones complejas con otras herramientas como la DB de Ekon.

Ejemplo, el siguiente script es una consulta a Ekon sobre las valores de tau medidos en puncion lumbar para un sujeto con un NHC determinado,

#!/bin/bash
nhc=$1
shift
sqlcmd -U osotolongo -P *********** -S 172.26.2.161 -s "," -W -Q "SELECT ficha.his_interno as Interno, pr.xdtform as FechaLCR, pr.xlcr_ttau as T_tau_LCR, pr.xlcr_ptau as P_tau_LCR, pr.xlcr_abeta40 as Ratio_Abeta42_40_LCR FROM [UNIT4_DATA].[imp].[his_pr_lcr] as pr JOIN  [UNIT4_DATA].[imp].[vh_pac_gral] as ficha ON (ficha.xpaciente_id = pr.xpaciente_id) WHERE ficha.his_interno = '"${nhc}"';"

En el caso concreto que nos interesa, este NHC coincide con el identificar de sujeto del proyecto (mriface), pero en cualquier caso no es dificil hacer una traslacion directa. Asi que basta con bajar la lista de sujetos del proyecto y directamente preguntar a Ekon por cada uno de ellos.

[osotolongo@brick03 mri_face]$ for x in `xnatapic list_subjects --project_id mriface --label | awk -F',' '{print $2}'`; do ./query_pl.sh ${x} | head -n 3 | tail -n 1 | grep -v "^$"; done > mriface_lcr.csv
[osotolongo@brick03 mri_face]$ head mriface_lcr.csv
20211480,2021-11-22 17:13:08.000,924,165,.0470
20211611,2021-12-10 10:33:12.000,270,39,.0560
20211438,2021-10-21 11:03:58.000,1410,268,.0370
20211637,2021-12-02 15:45:29.000,998,168,.0440
20211281,2021-11-22 16:55:27.000,813,142,.0420
20211401,2021-11-26 14:52:52.000,1301,208,.0510
20211683,2022-01-10 16:45:57.000,550,90,.0420
20211595,2021-12-17 13:25:06.000,532,88,.0500
20211627,2022-02-04 10:52:39.000,697,120,.0380
20211720,2021-12-22 12:38:02.000,460,92,.0410

MRI en XNAT

Hay dos pipelines incluidos en XNAT que realizan el tratamiento basico de las imagenes T1w y T2w. Para casi todos los proyectos esto es mas que suficiente. Sin embargo, el procesamiento mas avanzado, aqune puede resultar mas interesante, no se hace por defecto pues requiere interaccion con el investigador.

Freesurfer

Para el procesamiento general con FS no hay que hacer mas que incluir y configurar correctamente en el proyecto el pipeline RunFreesurfer. Si se configura correctamente, este pipeline se lanzara cada vez que se sube una MRI y los resultados se almacenan como un resource en el pexperimento correspondiente.

No obstante, en caso de ocurrir algun error, es posible relanzar el pipeline. Esto puede hacerse directamente en la web o usando xnatapic,

[osotolongo@brick03 mri_face]$ xnatapic run_pipeline --project_id unidad --pipeline RunFreesurfer --experiment_id XNAT_E00635 --dcmT1tag tfl3d1_16 --dcmT2tag spcir_220ns

pgs WMH

Ver WMH pgs

Este pipeline ejecuta un proceso muy intensivo y, en principio, se recomienda no ejecutarlo automaticamente. Pero puede cpnfigurarse para que asi sea.

Para ejecutar el pipeline se han de averiguar primero los sujetos que no se han hecho (opcional pero recomendable), y lanzarlo solo sobre estos.

[osotolongo@brick03 mri_face]$ for x in `xnatapic list_experiments --project_id mriface --modality MRI`; do if [[ -z `curl -f -X GET -b "JSESSIONID=FE98FA6B6AB5520A3A329757C3BC1127" "http://detritus.fundacioace.com:8088/data/experiments/${x}/files?format=json" 2>/dev/null| jq '.ResultSet.Result[].Name' | grep wmh.json` ]]; then echo ${x}; fi; done > falta_wmh.list
[osotolongo@brick03 mri_face]$ for x in `cat falta_wmh.list`; do xnatapic run_pipeline --project_id mriface --pipeline pgsWMH --experiment_id ${x}; done

Los resultados quedan almacenados en el resource WMH/wmh.json.

FSQC en XNAT

This part is obsolete and will be remove in further releases

TL;DR: $ review_vqcfs.pl -x <project>

La idea es que la segmentacion de Freesurfer la estoy haciendo en XNAT pero quiero hacer un quality control decente. Hemos escogido VisualQC pues proporciona informacion visual detallada de la segmentacion y tiene una interfaz amigable

enlace original: Realizar el QC a la segmentacion de Freesurfer de los sujetos analizados en XNAT
script en github: review_vqcfs.pl

Haciendo el QC con visualQC

El QC de la segmentacion lo hacemos con VisualQC. Para esto hay que entrar en el cluster con interfaz grafica y llamar el programa con algo parecido a,

[osotolongo@brick03 test_xnat_fsqc]$ vqcfs -i f5cehbi_experiments.list -f f5cehbi_fsresults -o visualqc_output

donde f5cehbi_experiments.list es un archvo que contiene una lista de las segmentaciones a revisar, f5cehbi_fsresults es el directorio que contiene los datos de Freesurfer y visualqc_output es donde generaremos los resultados. La base de este software esta ya descrita pero la ayuda online es muy buena.

Guardando los resultados

Después que hayamos revisado las imagenes segmentadas tendremos que guardar en XNAT los resultados del QC. Esto se hace con la opcion upload_fsqc, indicando la ubicacion del direcotrio output de visualqc_freesurfer,

[osotolongo@brick03 test_xnat_fsqc]$ xnatapic upload_fsqc --qcdir visualqc_output

Los resultados del QC quedan almacenados dentro de la informacion de cada MRI en XNAT.

En particular, el archivo rating.json almacena los resultados del rating.

{
"ResultSet":
  {
  "Result":
    [
      {
        "rating":"gOOD",
        "notes":"None"
      }
    ]
  }
}

De esta manera los resultados quedan guardados y pueden consultarse en cualquier momento, ya sea visualmente mediante la interfaz web de XNAT, o programaticamente.

un poco de logica (review FSQC)

Ahora lo que quiero hacer es algo para lo que visualQC no ha sido pensado. Quiero lanzar una herramienta que mire si los experimentos estan ya bajados y pasarle un lista de lo que quiero revisar. Si no paso la lista los revisa todos, y si la paso toma el directorio, escoge los que quiero, los reviso y suvbe los resultados.

Esto es contraintuitivo, en relacion a los scripts de xnatapic que bajan y procesan el QC.

Voy a intentar revisar el FSQC despues incluso que se haya hecho. Tengo que hacerle un frontend a esto.

Importante: el frontend hay que lanzarlo en un entorno grafico pues ademas de conectarse a XNAT y hacer cosas por debajo, va a lanzar la herramienta de evaluacion visual. He intentado regular esto añadiendo,

die "Should run this under X!\n" unless $ENV{'DISPLAY'};

pero todavia asi hay escenarios en que puede que este fijada la variable $DISPLAY y el script no encuentre el entorno grafico. O puedo fijar el variable $DISPLAY a un entorno grafico que ya exista desde una terminal que no lo tenga y me aparecera alli.

Para el frontend tengo que contemplar la posibilidad de empezar de cero. Vamos, lo basico que es lanzarlo e intentar procesar todo. Es decir que si hago solo,

$ review_vqcfs.pl -x f5cehbi

el script intentara bajar TODO lo que exista en el proyecto f5cehbi, lanzará vqcfs para todos los que existan y despues subirá los resultados.

Luego puedo definir el directorio donde estan los sujetos y entonces no tendre que bajarme nada, sino leer los datos que esten alli.

Consultando el QC

Se puede obtener un CSV con los resultados del QC haciendo,

[osotolongo@brick03 test_xnat_fsqc]$ xnatapic get_fsqc --project f5cehbi

El archivo de output se puede cambiar con las opciones (ver la ayuda) y obtenemos un archivo CSV de este tipo,

[osotolongo@brick03 test_xnat_fsqc]$ head f5cehbi_fsqc.csv 
F005, "gOOD", "None"
F006, "gOOD", " Pues parece que est aben"
F007, "rEVIEW LATER", ""
F009, "MInOR error", " not realy an error"
F010, "gOOD", ""
F014, "rEVIEW LATER", ""
F015, "gOOD", ""
F001, "gOOD", "None"
F100, "gOOD", "None"
F145, "I'M TIREd", "None"

Read Integrando ENIGMA Freesurfer QC instead.

Obteniendo resultados

Bajando resultados de Freesurfer

El pipeline de Freesurfer, deja los resultados en forma de archivos stats en el recurso FS. Uilizando la interfaz XNATACE es muy sencillo construir un script que baje y ordene estos recursos.

El script xnat_pullfs.pl es capaz de bajar estos archivos.

[osotolongo@brick03 mix_all]$ xnat_pullfs.pl -s aseg -d -x mriface -o mriface_aseg.csv
[osotolongo@brick03 mix_all]$ xnat_pullfs.pl -s aparc -x mriface -o mriface_aparc.csv
[osotolongo@brick03 mix_all]$ join -t, mriface_aseg.csv mriface_aparc.csv > mriface_mri.csv

Bajando WMH

Los resultados de WMH quedan almacenados en su propio recurso como un archivo json. Es asi muy sencillo obtenerlos directamente. El script xnat_get_wmh.pl baja y organiza estos valores.

[osotolongo@brick03 mri_face]$ xnat_get_wmh.pl -x mriface -o mriface_wmh.csv

PET-FBB en XNAT

Este pipeline tampoco se recomienda ejecutarlo automaticamente ya que depende de las imagenes PET y MRI y estas se entregan en fechas no relacionadas por dos proveedores distintos.

Ejecucion

Primero encontramos los experimentos correctos,

[osotolongo@brick03 facehbi]$ xnatapic list_experiments --project_id facehbi --type | grep petSessionData | awk -F"," {'print $1'} > xnat_pet_experiments.list
[osotolongo@brick03 facehbi]$ head xnat_pet_experiments.list
XNAT5_E00276
XNAT5_E00278
XNAT5_E00283
XNAT5_E00296
XNAT5_E00299
XNAT5_E00308
XNAT5_E00316
XNAT5_E00318
XNAT5_E00319
XNAT5_E00323

Y ahora mandamos todo,

[osotolongo@brick03 facehbi]$ for x in `cat xnat_pet_experiments.list`; do xnatapic run_pipeline --project_id facehbi --pipeline RegisterPETwithMRImatch --experiment_id ${x}; done

Extraer resultados

TL;DR: $ xnat_pullcl.pl -p <project> -d -o <output_file.csv>

La idea fundamental aqui es obtener los valores de SUVR y centiloide calculados en XNAT. Estos pueden obtenerse con,

$ xnatapic get_registration_report --project_id f5cehbi --output xnat_pet_results.csv

y esto vamos a limpiarlo un poco,

$ awk -F"," '{print $2","$6","$7}' xnat_pet_results.csv | sed 's/"//g' | tail -n +2 | sort -t, -k 1 > pet.results

La lista de sujetos del proyecto se obtiene con,

$  xnatapic list_subjects --project_id f5cehbi --label > xnat_subjects.list

Esto se puede ordenar rapidamente haciendo,

$ sort -t, -k 1 xnat_subjects.list > xnat_subjects_sorted.list

Y entonces podemos integrar los ID de proyecto con los resultaos, a traves de los ID de XNAT, y seleccionar los datos que queremos guardar,

$ join -t, xnat_subjects_sorted.list pet.results | sort -t, -k 2 | awk -F"," '{if ($3) print $2";"$3";"$4}' | sed '1iSubject;SUVR;Centiloid'"

Todo esto se resume en el script xnat_pullcl.pl. Los archivos intermedios se borran una vez que se obtienen los resultados. Las opciones son,

-p : nombre del proyecto local
-x : nombre del proyecto en XNAT
-o : archivo de salida (default STDOUT)

Troubleshooting

Sacando los PET a ejecutar

¿Que sujetos tengo en el proyecto?

[osotolongo@brick03 f5cehbi]$ xnatapic list_subjects --project_id f5cehbi --label > xnat_subjects.list

¿De aqui, que sujetos tienen MRI?

[osotolongo@brick03 f5cehbi]$ for x in `awk -F"," {'print $1'} xnat_subjects.list`; do e=$(xnatapic list_experiments --project_id f5cehbi --subject_id ${x} --modality MRI); if [[ ${e} ]]; then echo "${x},${e}"; fi; done > xnat_subject_mri.list

¿De estos, cuales tiene PET?

[osotolongo@brick03 f5cehbi]$ for x in `awk -F"," {'print $1'} xnat_subject_mri.list`; do e=$(xnatapic list_experiments --project_id f5cehbi --subject_id ${x} --modality PET); if [[ ${e} ]]; then echo "${x},${e}"; fi; done > xnat_subjects_mri_pet.list

¿De estos PET, cuantos ha sido calculados?

[osotolongo@brick03 f5cehbi]$ awk -F',' '{if($6!="") print $3}' xnat_pet_results.csv | sed 's/"//g' | tail -n +2 > already_done.txt

Y por ultimo, ¿Cuales me falta por calcular?

[osotolongo@brick03 f5cehbi]$ grep -v "`cat already_done.txt`" xnat_subjects_mri_pet.list | awk -F',' '{print $2}' > pet2do.list

Vamos a hacerlos entonces!

[osotolongo@brick03 f5cehbi]$ for x in `cat pet2do.list`; do xnatapic run_pipeline --project_id f5cehbi --pipeline RegisterPETwithMRImatch --experiment_id ${x}; done
[osotolongo@brick03 f5cehbi]$ queue
   JOBID PARTITION                                 NAME     USER ST       TIME  NODES NODELIST(REASON)
   13027      fast           RunFreesurfer.XNAT5_E00001     xnat  R    4:35:39      1 brick01
   13028      fast RegisterPETwithMRImatch.XNAT5_E00710     xnat  R       0:15      1 brick01
   13029      fast RegisterPETwithMRImatch.XNAT5_E00729     xnat  R       0:15      1 brick01
   13030      fast RegisterPETwithMRImatch.XNAT5_E00733     xnat  R       0:15      1 brick01
   13031      fast RegisterPETwithMRImatch.XNAT5_E00735     xnat  R       0:12      1 brick01
   13032      fast RegisterPETwithMRImatch.XNAT5_E00737     xnat  R       0:12      1 brick01
   13033      fast RegisterPETwithMRImatch.XNAT5_E00742     xnat  R       0:12      1 brick01
   13034      fast RegisterPETwithMRImatch.XNAT5_E00743     xnat  R       0:12      1 brick02
   13035      fast RegisterPETwithMRImatch.XNAT5_E00746     xnat  R       0:12      1 brick02
   13036      fast RegisterPETwithMRImatch.XNAT5_E00751     xnat  R       0:12      1 brick02
   13037      fast RegisterPETwithMRImatch.XNAT5_E00752     xnat  R       0:12      1 brick02
   13038      fast RegisterPETwithMRImatch.XNAT5_E00753     xnat  R       0:12      1 brick02

PET con tag erroneo

Si el tag con que estan etiquetadas las imagenes es incorrecto, el registro de PET dara error. Tipicamente recibes un email con un subject como:

XNAT update: Processing failed for F082A

Habria que revisar el sujeto, mirar como se ha etiquetado la imagen,

y relanzar algo como,

[osotolongo@brick03 f5cehbi]$ xnatapic run_pipeline --project_id f5cehbi --pipeline RegisterPETwithMRImatch --experiment_id XNAT5_E00752 --dcmFBBtag PETNAC_Standard
[osotolongo@brick03 f5cehbi]$ queue
   JOBID PARTITION                                 NAME     USER ST       TIME  NODES NODELIST(REASON)
   13039      fast RegisterPETwithMRImatch.XNAT5_E00752     xnat  R       0:03      1 brick01

failsafe PET

Un escenario que puede ocurrir es cuando el DICOM PET no se convierte correctamente con las herramientas automaticas. Es posible entonces, subir el PET en formato NIfTI-1 y ejecutar un pipeline para que registre este ultimo archivo.

[osotolongo@brick03 f2cehbi]$ xnatapic upload_nifti --project_id f2cehbi --subject_id F043 --experiment_id F043F --scan_id 5 /nas/clinic/nii/v2/F043/PET.nii PET.json

El archivo nifti debe ser el correspondiente al total de tiempo integrado,

[osotolongo@brick03 f2cehbi]$ fslinfo /nas/clinic/nii/v2/F043/PET.nii
data_type      FLOAT32
dim1           400
dim2           400
dim3           109
dim4           1
datatype       16
pixdim1        1.018210
pixdim2        1.018210
pixdim3        2.027008
pixdim4        0.000000
cal_max        0.0000
cal_min        0.0000
file_type      NIFTI-1+

El archivo json debe contener la info minima para identificar el protocolo correcto,

[osotolongo@brick03 f2cehbi]$ cat PET.json 
{
	"SeriesDescription": "FACEHBI_Florbetaben_20min",
	"ProtocolName": "FACEHBI_Florbetaben_20min"
}

ambos archivos deben tener el mismo nombre.

Tras esto ya es posible lanzar el pipeline RegisterPETwithMRI,

[osotolongo@brick03 f2cehbi]$ xnatapic run_pipeline --project_id f2cehbi --pipeline RegisterPETwithMRI --experiment_id F043F --dcmT1tag tfl3d1_16 --dcmFBBtag _20min

PET FBB QC interactuando con XNAT

Making of

El procesamiento de las FBB en XNAT genera una imagen para evaluar la calidad del registro de FBB a espacio MNI.

lo que necesitamos es una herramienta sencilla capaz de leer la información de XNAT, mostrar las imagenes, permitir evaluar el QC y volcar de regreso la info en XNAT.

TIMTOWTDI, pero voy a hacerlo a la old fashion con CGI.pm.

Necesito leer tres cosas por cada sujeto, experiment_ID, imagen de QA y estado del QA. Primero la lista de sujetos del proyecto, con la REST API,

curl -X GET -u "$USER:$PASSWD" "$HOST/data/projects/$PROJECT/subjects?format=json" | jq ".ResultSet.Result[].ID" | sed \'s/"/,/g\'

Ahora por cada sujeto saco el experimento PET que le corresponde,

curl -X GET -u "$USER:$PASSWD" "$HOST/data/projects/$PROJECT/subjects/$SUBJECT/experiments?xsiType=xnat:petSessionData" 2>/dev/null | jq ".ResultSet.Result[].ID" | sed \'s/"//g\'

Para este experimento obtengo el QA,

curl -X GET -u "$USER:$PASSWD" "$HOST/data/experiments/$EXPERIMENT/resources/MRI/files/mriSessionMatch.json" 2>/dev/null | jq ".ResultSet.Result[].qa"

y me bajo la imagen a un directorio temporal,

curl -X GET -u "$USER:$PASSWD" "$HOST/data/experiments/$EXPERIMENT/resources/MRI/files/$SUBJECT_fbb_mni.gif" -o $TMP_IMAGE

Todo esto lo voy metiendo en un hash y lo pongo en un formulario con la CGI. Cuando se envia el formulario tengo que crear el json con la info y luego subirlo. Pero lo que voy a hacer es bajar el json existente, cambiar el valor de QA y volverlo a subir. Algo asi,

curl -X GET -u "$USER:$PASSWD" "$HOST/data/experiments/$EXPERIMENT/resources/MRI/files/mriSessionMatch.json" 2>/dev/null

esto lo meto en una variable y le hago,

$resp =~ s/,"qa":\d,/"qa":$qcs{$xpx},/;

la escribo a un archivo temporal y la subo con,

curl -X PUT -u "$USER:$PASSWD" "$HOST/data/experiments/$EXPERIMENT/resources/MRI/files/mriSessionMatch.json?overwrite=true" -F file="@$TMP_FILE"

y ya esta!!!!! ponemos un smiley o algo. el codigo en github.

Usando esto

El output es feo y lento, pero facil de usar. Primero hay que suministrar el nombre del proyecto en XNAt y los datos de login para XNAT.

La aplicacion se conectara entonces usando la REST API y mostrara el resultado del QC que tenemos. Aqui se modifica los necesario, se le da a submit y ya esta.

Nota: no hay que ser muy exigente con el resutado del registro. Hay mucho smooth y ROIs grandes para compensar los pequeños desplazamiento. on que la cabeza este mas o menos en su sitio ya nos vale.

Bien	Aceptable	Mal

Trabajando en local

Manipular las rutas

La configuración del proyecto puede ser modificada editando el archivo de configuración directamente o utilizando la herramienta pipeconf.pl, que puede mostrar, editar una variable o incluso crear un nuevo archivo con los valores por defecto. Ejemplos:,

Mostrar configuración

$ pipeconf.pl -p f5cehbi -r
BIDS = /old_nas/f5cehbi/bids
DATA = /old_nas/f5cehbi
PET = /nas/data/f5cehbi/raw_tau/PI2620
SRC = /nas/data/f5cehbi/raw
WORKING = /old_nas/f5cehbi/working
XNAME = f5cehbi

Editar configuración

$ pipeconf.pl -p f5cehbi -a XNAME=f2cehbi
$ pipeconf.pl -p f5cehbi -r
BIDS = /old_nas/f5cehbi/bids
DATA = /old_nas/f5cehbi
PET = /nas/data/f5cehbi/raw_tau/PI2620
SRC = /nas/data/f5cehbi/raw
WORKING = /old_nas/f5cehbi/working
XNAME = f2cehbi

Crear nuevo archivo

$ pipeconf.pl -p test -c /nas/data/test
[osotolongo@brick03 acenip]$ pipeconf.pl -p test -r
BIDS = /nas/data/test/bids
DATA = /nas/data/test
SRC = /nas/data/test/raw
WORKING = /nas/data/test/working
XNAME = test

Con este conjunto de operaciones basta para manipular la configuración de cualquier proyecto.

de XNAT a local

Aunque se supone que tenemos todos los DICOM en el server, lo cierto es que lo que nos suben puede estar muy desordenado. Para llevar las imagenes a XNAT ya nos hemos tomado el trabajo de ordenar todo y en XNAT deberia estar el proyecto muy ordenado.

Luego, la idea es que una vez que el proyecto en XNAT este cerrado, podemos recrear la colección de DICOMs en el directorio SRC del proyecto local. El script xnat2raw.pl recorre los sujetos del projecto XNAT y baja los DICOM de cada experimento (MRI o PET) y los guarda en un directorio diferenciado utilizando la funcion xget_dicom de la biblioteca de interacción con la API (Ver markdown docs).

Por ejemplo, para el proyecto local unidad, el archivo de configuración,

[osotolongo@brick03 mri_face]$ cat ~/.config/neuro/unidad.cfg
BIDS = /old_nas/mri_face/bids
DATA = /old_nas/mri_face
SRC = /old_nas/mri_face/raw
WORKING = /old_nas/mri_face/working
XNAME = tester

nos dice que debe enlazarse con el proyecto XNAT tester. Asi que al ejecutar,

[osotolongo@brick03 mri_face]$ xnat2raw.pl -p unidad

se crea un directorio por sujeto del proyecto dentro del SRC,

[osotolongo@brick03 mri_face]$ ls raw/
20050456  20081210  20100678

y el archivo CSV, unidad_MRI.csv,

[osotolongo@brick03 mri_face]$ cat unidad_MRI.csv
0001,20081210
0002,20100678
0003,20050456

que podemos utilizar para crear facilmente la DB del proyecto y ejecutar la conversion a BIDS.

DCM2BIDS

El proceso de conversion de los DICOM a formato BIDS ha de hacerse manualmente. Esto es, aunque con ayuda de herramientas (Convertir DICOMs a BIDS), la conversion necesita que se realicen ciertos pasos no automatizados. Dado la disimilitud de los proyectos, los archivos han de editarse y el proceso de conversion de al menos uno de los sujetos ha de revisarse para automatizar el resto.

El primer paso para la conversión ha de ser ejecutar dcm2bids_scaffold dentro del directorio bids del proyecto. Esto crea la estructura necesaria para la ejecucion.

[tester@detritus bids]$ dcm2bids_scaffold
[tester@detritus bids]$ ls
CHANGES  code  dataset_description.json  derivatives  participants.json  participants.tsv  README  sourcedata

El archivo dataset_description.json ha de editarse manualmente, pero es bastante sencillo,

[tester@detritus bids]$ cat dataset_description.json
{
  "Name": "PTEST",
  "BIDSVersion": "1.2.0",
  "License": "CC0"
}

Ahora debemos escoger un sujeto del proyecto para examinar la estructura de los datos,

[tester@detritus bids]$ dcm2bids_helper -d /nas/tester/testproj/032-00001-SCR-GLOBAL/
Example in:
/nas/data/ptest/bids/tmp_dcm2bids/helper

La estructura dejada por este comando, junto al protocolo de maquina permite llenar el archivo de configuracion necesario para la conversion de los sujetos. Esto ha de hacerse con cuidado, pero una vez rellenado correctamente la conversion es sencilla.

Ejemplo basado en el protocolo de EPAD

{
        "descriptions": [
                {        
                        "dataType": "anat",
                        "modalityLabel": "T1w",
                        "criteria": {
                                "SeriesDescription": "MPRAGE_GRAPPA2",
                                "ImageType": ["ORIGINAL", "PRIMARY", "M", "NORM", "DIS3D", "DIS2D", "MFSPLIT"]
                        }
                },
                {        
                        "dataType": "anat",
                        "modalityLabel": "T2w",
                        "criteria": {
                                "SeriesDescription": "3D-T2-FLAIR_SAG*",
                                "ImageType": ["ORIGINAL", "PRIMARY", "M", "NORM", "DIS3D", "DIS2D", "MFSPLIT"]
                        }
                },
                {
                        "dataType": "dwi",
                        "modalityLabel": "dwi",
                        "criteria": {
                                "SeriesDescription": "EPAD-SingleShell-DTI48",
                                "ImageType": ["ORIGINAL", "PRIMARY", "DIFFUSION", "NONE", "ND", "NORM", "MFSPLIT"]
                        }
                },
                {
                        "dataType": "dwi",
                        "modalityLabel": "dwi",
                        "customLabels": "dir-PA_sbref",
                        "criteria": {
                                "SeriesDescription": "EPAD-B0-RevPE"
                        }
 
                },
                {
                        "dataType": "func",
                        "modalityLabel": "bold",
                        "customLabels": "task-rest",
                        "criteria": {
                                "SeriesDescription": "EPAD-rsfMRI*",
                                "ImageType": ["ORIGINAL", "PRIMARY", "M", "ND", "NORM", "MFSPLIT"]
                        }
                },
                {
                        "dataType": "func",
                        "modalityLabel": "bold",
                        "customLabels": "task-rest_dir-PA_sbref",
                        "criteria": {
                                "SeriesDescription": "EPAD-SE-fMRI"
                        }
                },
                {
                        "dataType": "func",
                        "modalityLabel": "bold",
                        "customLabels": "task-rest_dir-AP_sbref",
                        "criteria": {
                                "SeriesDescription": "EPAD-SE-fMRI-RevPE"
                        }
                }
 
        ]
}

En este punto ya podemos ejecutar la conversion de cualquier sujeto usando dcm2bids,

$ dcm2bids -d /src_path/ -p subj_id -c config_file.json -o proj_path/bids/

Ejemplo

[tester@detritus bids]$ dcm2bids -d /nas/tester/testproj/032-00001-SCR-GLOBAL/ -p 0001 -c /nas/data/ptest/bids/conversion.json -o /nas/data/ptest/bids/
INFO:dcm2bids.dcm2bids:--- dcm2bids start ---
INFO:dcm2bids.dcm2bids:OS:version: Linux-2.6.32-431.17.1.el6.centos.plus.x86_64-x86_64-with-centos-6.10-Final
INFO:dcm2bids.dcm2bids:python:version: 3.4.8 (default, Apr  9 2018, 11:43:18) [GCC 4.4.7 20120313 (Red Hat 4.4.7-18)]
INFO:dcm2bids.dcm2bids:dcm2bids:version: 2.1.4
INFO:dcm2bids.dcm2bids:dcm2niix:version: v1.0.20191031
INFO:dcm2bids.dcm2bids:participant: sub-0001
INFO:dcm2bids.dcm2bids:session:
INFO:dcm2bids.dcm2bids:config: /nas/data/ptest/bids/conversion.json
INFO:dcm2bids.dcm2bids:BIDS directory: /nas/data/ptest/bids
INFO:dcm2bids.utils:Running dcm2niix -b y -ba y -z y -f '%3s_%f_%p_%t' -o /nas/data/ptest/bids/tmp_dcm2bids/sub-0001 /nas/tester/testproj/032-00001-SCR-GLOBAL/
INFO:dcm2bids.dcm2niix:Check log file for dcm2niix output
INFO:dcm2bids.sidecar:Sidecars pairing:
............
 
[tester@detritus bids]$ tree sub-0001/
sub-0001/
├── anat
│   ├── sub-0001_T1w.json
│   ├── sub-0001_T1w.nii.gz
│   ├── sub-0001_T2w.json
│   └── sub-0001_T2w.nii.gz
├── dwi
│   ├── sub-0001_dir-PA_sbref_dwi.json
│   ├── sub-0001_dir-PA_sbref_dwi.nii.gz
│   ├── sub-0001_dwi.bval
│   ├── sub-0001_dwi.bvec
│   ├── sub-0001_dwi.json
│   └── sub-0001_dwi.nii.gz
└── func
    ├── sub-0001_task-rest_bold.json
    ├── sub-0001_task-rest_bold.nii.gz
    ├── sub-0001_task-rest_dir-AP_sbref_bold.json
    ├── sub-0001_task-rest_dir-AP_sbref_bold.nii.gz
    ├── sub-0001_task-rest_dir-PA_sbref_bold.json
    └── sub-0001_task-rest_dir-PA_sbref_bold.nii.gz
 
3 directories, 16 files

Lo que es deseable es ejecutar la conversion de todos los sujetos del proyecto usando el cluster. ahora bien,una vez que el archivo de configuracion se ha editado correctamente ya se puede utilizar una herramienta generica.

Lo programamos entonces en el pipeline

[tester@detritus ptest]$ bulk2bids.pl ptest
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114181     devel fs_recon   tester PD       0:00      1 (Dependency)
            114191     devel dcm2bids   tester PD       0:00      1 (Dependency)
            114184     devel dcm2bids   tester  R       0:42      1 brick01
            114185     devel dcm2bids   tester  R       0:42      1 brick01
            114186     devel dcm2bids   tester  R       0:42      1 brick01
            114187     devel dcm2bids   tester  R       0:42      1 brick01
            114188     devel dcm2bids   tester  R       0:42      1 brick01
            114189     devel dcm2bids   tester  R       0:42      1 brick01
[tester@detritus ptest]$ tail -f slurm/dcm2bids0002-114184
INFO:dcm2bids.dcm2bids:--- dcm2bids start ---
INFO:dcm2bids.dcm2bids:OS:version: Linux-3.10.0-862.11.6.el7.x86_64-x86_64-with-centos-7.5.1804-Core
INFO:dcm2bids.dcm2bids:python:version: 3.6.8 (default, Aug  7 2019, 17:28:10) [GCC 4.8.5 20150623 (Red Hat 4.8.5-39)]
INFO:dcm2bids.dcm2bids:dcm2bids:version: 2.1.4
INFO:dcm2bids.dcm2bids:dcm2niix:version: v1.0.20191031
INFO:dcm2bids.dcm2bids:participant: sub-0002
INFO:dcm2bids.dcm2bids:session:
INFO:dcm2bids.dcm2bids:config: /home/data/ptest/bids/conversion.json
INFO:dcm2bids.dcm2bids:BIDS directory: /home/data/ptest/bids
INFO:dcm2bids.utils:Running dcm2niix -b y -ba y -z y -f '%3s_%f_%p_%t' -o /nas/data/ptest/bids/tmp_dcm2bids/sub-0002 /nas/tester/testproj/032-00002-SCR-GLOBAL/
.........

Issues

A veces dmc2niix falla y se interrumpe la ejecucion de dcm2bids. Workaround: Separar los procesos de dcm2niix y dcm2bids. Primero convierto y despues organizo.De paso tengo mas control sobre como convierto (Ver respuesta de Chris Rorden)

        my $order = "mkdir -p $std{'DATA'}/bids/tmp_dcm2bids/sub-$subject; dcm2niix -i y -b y -ba y -z y -f '%3s_%f_%p_%t' -o $std{'DATA'}/bids/tmp_dcm2bids/sub-$subject $std{'SRC'}/$guys{$subject}/; dcm2bids -d $std{'SRC'}/$guys{$subject}/ -p $subject -c $std{'DATA'}/$cfile -o $std{'DATA'}/bids/";

Preparando el proyecto local

Para ejecutar tareas complejas necesitamos crear un proyecto local. Este puede corresponder a un proyecto completo de XNAT o a un subconjunto.

Lo primero es crear la configuracion y el esqueleto del proyecto con algunos valores por defecto,

[osotolongo@brick03 ~]$ make_proj.pl example /old_nas/osotolongo/example

Ahora podemos utilizar la utilidad de edicion de proyectos, ya sea para leer la configuracion,

[osotolongo@brick03 ~]$ pipeconf.pl -r -p example
DATA = /ruby/neuroimagen/example
SRC = /old_nas/osotolongo/example
WORKING = /ruby/neuroimagen/example/working
BIDS = /ruby/neuroimagen/example/bids

o, por ejemplo, para añadir el proyecto de XNAT que debemos enlazar,

[osotolongo@brick03 ~]$ pipeconf.pl -a XNAME=unidad -p example
[osotolongo@brick03 ~]$ pipeconf.pl -r -p example
BIDS = /ruby/neuroimagen/example/bids
DATA = /ruby/neuroimagen/example
SRC = /old_nas/osotolongo/example
WORKING = /ruby/neuroimagen/example/working
XNAME = unidad

Ahora debemos preparar el esqueleto de la conversion a BIDS,

[osotolongo@brick03 ~]$ cd /ruby/neuroimagen/example/bids
[osotolongo@brick03 bids]$ dcm2bids_scaffold 
[osotolongo@brick03 bids]$ ls
CHANGES  code  dataset_description.json  derivatives  participants.json  participants.tsv  README  sourcedata

y deberiamos crear el archivo conversion.json con los datos de los scans que querramos convertir a BIDS. Esto es bastante simple mirando la informacion de la MRI en XNAT,

Por ejemplo, si de estos experimentos queremos los T1w y T2w, creamos el archivo,

conversion.json

{
	"descriptions" :[
		{
			"dataType": "anat",
			"modalityLabel": "T1w",
		       "criteria" : 
		       {
			       "SeriesDescription" : "t1_mprage_sag_p2_iso"
		       }	       
		},
		{
			"dataType": "anat",
			"modalityLabel": "T2w",
			"criteria" :
			{
				"SeriesDescription": "t2_space_FLAIR_sag_p2_ns-t2prep"
			}
		}
	]
}

Entonces, de aqui puedo leer directamente los scans que tengo que bajar por cada experimento haciendo algo como,

[osotolongo@brick03 bids]$ cat conversion.json | jq '.descriptions[].criteria.SeriesDescription'
"t1_mprage_sag_p2_iso"
"t2_space_FLAIR_sag_p2_ns-t2prep"

XNAT to BIDS

Para llevar un grupo de sujetos de XNAT a un proyecto local vamos a establecer un procedimiento operativo sencillo. Primero, vamos a asumir que todos los sujetos pertenecen al mismo proyecto de XNAT y que hemos ya configurado correctamente el proyecto local (lee mas arriba, anda).

[osotolongo@brick03 example]$ pipeconf.pl -r -p example
BIDS = /ruby/neuroimagen/example/bids
DATA = /ruby/neuroimagen/example
SRC = /old_nas/osotolongo/example
WORKING = /ruby/neuroimagen/example/working
XNAME = unidad

ahora voy a hacer una lista de los sujetos del proyecto XNAT que quiero incluir en el proyecto local.

[osotolongo@brick03 example]$ cat sbj.list
20070391
20081210
20100678

En principio esto es todo lo que necesitaria para empezar. El esquema de trabajo es mas o menos asi,

Para hacer esto basta con enviar,

[osotolongo@brick03 example]$ xnat2bids.pl -p example -i sbj.list -j bids/conversion.json

Y, al filtrar los campos de Series Description y pasarlos a la funcion xget_dicom() se bajan solamente los archivos correspondientes a T1w y T2w (en este caso).

[osotolongo@brick03 example]$ ls -l /old_nas/osotolongo/example/*/*/scans
/old_nas/osotolongo/example/20070391/1D22133126/scans:
total 8
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 5001-t1_mprage_sag_p2_iso
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 6001-t2_space_FLAIR_sag_p2_ns_t2prep
 
/old_nas/osotolongo/example/20081210/1D21146858/scans:
total 8
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 5001-t1_mprage_sag_p2_iso
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 6001-t2_space_FLAIR_sag_p2_ns_t2prep
 
/old_nas/osotolongo/example/20100678/1D22003736/scans:
total 16
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 10-t1_mprage_sag_p2_iso
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 11-t2_space_FLAIR_sag_p2_ns_t2prep
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 19-t1_mprage_sag_p2_iso
drwx------ 3 osotolongo osotolongo 4096 Dec 23 09:10 5-t1_mprage_sag_p2_iso

Tras esto se lanza la conversion a BIDS en el schedule manager.

[osotolongo@brick03 example]$ ls -l bids/sub-000*/anat
bids/sub-0001/anat:
total 20072
-rw-rw-r-- 1 osotolongo osotolongo     1572 Dec 23 09:10 sub-0001_T1w.json
-rw-rw-r-- 1 osotolongo osotolongo 13451620 Dec 23 09:10 sub-0001_T1w.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo     1674 Dec 23 09:10 sub-0001_T2w.json
-rw-rw-r-- 1 osotolongo osotolongo  7089107 Dec 23 09:10 sub-0001_T2w.nii.gz
 
bids/sub-0002/anat:
total 48400
-rw-rw-r-- 1 osotolongo osotolongo     1586 Dec 23 09:10 sub-0002_run-01_T1w.json
-rw-rw-r-- 1 osotolongo osotolongo 13485861 Dec 23 09:10 sub-0002_run-01_T1w.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo     1570 Dec 23 09:10 sub-0002_run-02_T1w.json
-rw-rw-r-- 1 osotolongo osotolongo 14063513 Dec 23 09:10 sub-0002_run-02_T1w.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo     1570 Dec 23 09:10 sub-0002_run-03_T1w.json
-rw-rw-r-- 1 osotolongo osotolongo 14433653 Dec 23 09:10 sub-0002_run-03_T1w.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo     1673 Dec 23 09:10 sub-0002_T2w.json
-rw-rw-r-- 1 osotolongo osotolongo  7556939 Dec 23 09:10 sub-0002_T2w.nii.gz
 
bids/sub-0003/anat:
total 19220
-rw-rw-r-- 1 osotolongo osotolongo     1582 Dec 23 09:10 sub-0003_T1w.json
-rw-rw-r-- 1 osotolongo osotolongo 13091825 Dec 23 09:10 sub-0003_T1w.nii.gz
-rw-rw-r-- 1 osotolongo osotolongo     1685 Dec 23 09:10 sub-0003_T2w.json
-rw-rw-r-- 1 osotolongo osotolongo  6576547 Dec 23 09:10 sub-0003_T2w.nii.gz

Ademas, el archivo que enlaza los ID de procesamiento con los ID del proyecto en XNAT queda correctamente configurado,

[osotolongo@brick03 example]$ cat example_mri.csv
0001;20081210
0002;20100678
0003;20070391

Además, los archivos bids/participants.tsv y bids/participants.json se rellenan con toda la informacion necesaria para conectar la estructura BIDS con los datos de XNAT,

[osotolongo@brick03 example]$ cat bids/participants.tsv
0001	20081210
0002	20100678
0003	20070391
[osotolongo@brick03 example]$ cat bids/participants.json  | jq '.'
{
  "0001": {
    "XNAT_experiment_ID": "XNAT05_E00046",
    "XNAT_ID": "XNAT05_S00033",
    "XNAT_label": "20081210"
  },
  "0002": {
    "XNAT_experiment_ID": "XNAT05_E00094",
    "XNAT_label": "20100678",
    "XNAT_ID": "XNAT05_S00081"
  },
  "0003": {
    "XNAT_ID": "XNAT_S00722",
    "XNAT_label": "20070391",
    "XNAT_experiment_ID": "XNAT_E00839"
  }
}

Nota: Al reducir tanto el download como la conversion a los archivos que realmente necesitamos, el tiempo de procesamiento a BIDS disminuye dramaticamente.

GET fake FS subject's pool

Mucho del procesamiento avanzado, por ejemplo el procesamiento DTI, necesita no solo la estructura BIDS sino que toma la MRI postprocesada de Freesurfer. No hemos de ejecutar de nuevo la segmentacion ya que todos los archivos del procesamiento estan almacenados en XNAT. Asi que nos bajamos todo y lo organizamos.

Basicamente lo unico que habria que hacer es leer el archivo de enlaces del proyecto (example_mri.csv) y para cada linea, tomar el ID del sujeto de XNAT, buscar el MRI correspondiente, y bajar el tgxz eu esta dentro del resource FS a descomprimir en ${SUBJECTS_DIR}. Bueno y algun tratamiento de strings por el camino.

El problema seria el tiempo que tarda en descomprimir esto cuando son muchos sujetos. Asi que vamos a ver como enviamos las tareas a SLURM. El uso de SLURM aqui no solo cambia el como se va a usar la interaccion con XNAT sino el uso de directorios temporal.

No es util usar xget_res_file() desde Perl y enviar a SLURM el untar porque es casi igualmente lento. Tengo que ver como paso la orden a SLURM (o sea a bash)
No puedo hacer un directorio temporal desde un nodo y usarlo desde otro ya que no sera accesible (esto sin importar el path!)

Puedo escribir directamente la interaccion con la API pero en lugar de eso voy a usar la funcion xget_res_file pero pasandole un parametro extra para que no ejecute sino que exporte la orden a ejecutar –> undocumented features

ahora, al hacer esto,

[osotolongo@brick03 example]$ xnat_getfs.pl -p example -ovw
Getting XNAT subject list
OK, now I'm going on
Downloading and extracting
D21582920.tar.gz -> /old_nas/subjects/example_0001
D22013217.tar.gz -> /old_nas/subjects/example_0002
D22488864.tar.gz -> /old_nas/subjects/example_0003

se crean los archivos sbatch como,

[osotolongo@brick03 example]$ cat slurm/0001.sh
#!/bin/bash
#SBATCH -J fake_FS_example
#SBATCH -o /ruby/neuroimagen/example/slurm/0001.out-%j
#SBATCH --mail-user=osotolongo
#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT
mkdir /ruby/user_data/osotolongo/tmp/example_0001
curl -f -X GET -b JSESSIONID=DCC76CADC18614B01DD7DFB16CDF31D7 "http://detritus.fundacioace.com:8088/data/experiments/XNAT05_E00046/resources/FS/files/D21582920.tar.gz" -o /ruby/user_data/osotolongo/tmp/example_0001/XNAT05_S00033.tar.gz
tar xzf /ruby/user_data/osotolongo/tmp/example_0001/XNAT05_S00033.tar.gz -C /old_nas/subjects/example_0001/ --transform='s/XNAT05_S00033//' --exclude='fsaverage'
rm -rf /ruby/user_data/osotolongo/tmp/example_0001

que son los encargados de,

crear el directorio temporal
bajar el archivo comprimido de FS
descompactarlo en la ruta final
borrar el directorio temporal

y despues ponemos el warning usual para que avise cuando termine. La ejecucion queda como,

[osotolongo@brick03 example]$ queue
   JOBID PARTITION                                 NAME     USER ST       TIME  NODES NODELIST(REASON)
  326038      fast                      fake_FS_example osotolon PD       0:00      1 (Dependency)
  326035      fast                      fake_FS_example osotolon  R       0:03      1 brick01
  326036      fast                      fake_FS_example osotolon  R       0:03      1 brick01
  326037      fast                      fake_FS_example osotolon  R       0:03      1 brick01

y nos debe dejar todo el procesamiento de FS en su lugar,

[osotolongo@brick03 example]$ ls /old_nas/subjects/example_*
/old_nas/subjects/example_0001:
label  mri  scripts  stats  surf  tmp  touch  trash
 
/old_nas/subjects/example_0002:
label  mri  scripts  stats  surf  tmp  touch  trash
 
/old_nas/subjects/example_0003:
label  mri  scripts  stats  surf  tmp  touch  trash

Hippocampal Subfields

(detalles en Segmentacion del hipocampo)

A esta version se le ha añadido la capacidad de extraer la segmentacion del hipocampo y la amigdala. Primero debemos tener el proyecto creado y la segmentacion de FS en local.

Ahora solo hacemos,

[osotolongo@brick03 mri_face]$ hsfrecon.pl mriface 
[osotolongo@brick03 acenip]$ queue
   JOBID PARTITION                                 NAME     USER ST       TIME  NODES NODELIST(REASON)
  218945      fast                    hsf_recon_mriface osotolon  R       2:16      1 brick02
  218946      fast                    hsf_recon_mriface osotolon  R       2:06      1 brick02
  218947      fast                    hsf_recon_mriface osotolon  R       1:56      1 brick02
  218948      fast                    hsf_recon_mriface osotolon  R       1:46      1 brick02
  218949      fast                    hsf_recon_mriface osotolon  R       1:36      1 brick02
  218950      fast                    hsf_recon_mriface osotolon  R       1:26      1 brick02
  218951      fast                    hsf_recon_mriface osotolon  R       1:16      1 brick02
  218952      fast                    hsf_recon_mriface osotolon  R       1:06      1 brick02
  218953      fast                    hsf_recon_mriface osotolon  R       0:56      1 brick02
  218954      fast                    hsf_recon_mriface osotolon  R       0:46      1 brick02
  218955      fast                    hsf_recon_mriface osotolon  R       0:36      1 brick02
  218956      fast                    hsf_recon_mriface osotolon  R       0:26      1 brick02
  218957      fast                    hsf_recon_mriface osotolon  R       0:16      1 brick02
  218958      fast                    hsf_recon_mriface osotolon  R       0:06      1 brick02
...
...
...

y, tras recibir el aviso de finalización,

se sacan las métricas,

[osotolongo@brick03 mri_face]$ hsf_metrics.pl mriface

con lo que los resultados quedan en mriface_hsf_metrics.csv.

[osotolongo@brick03 mri_face]$ head mriface_hsf_metrics.csv
Subject,lh.CA1-body,lh.CA1-head,lh.CA3-body,lh.CA3-head,lh.CA4-body,lh.CA4-head,lh.GC-ML-DG-body,lh.GC-ML-DG-head,lh.HATA,lh.Hippocampal_tail,lh.Whole_hippocampal_body,lh.Whole_hippocampal_head,lh.Whole_hippocampus,lh.fimbria,lh.hippocampal-fissure,lh.molecular_layer_HP-body,lh.molecular_layer_HP-head,lh.parasubiculum,lh.presubiculum-body,lh.presubiculum-head,lh.subiculum-body,lh.subiculum-head,rh.CA1-body,rh.CA1-head,rh.CA3-body,rh.CA3-head,rh.CA4-body,rh.CA4-head,rh.GC-ML-DG-body,rh.GC-ML-DG-head,rh.HATA,rh.Hippocampal_tail,rh.Whole_hippocampal_body,rh.Whole_hippocampal_head,rh.Whole_hippocampus,rh.fimbria,rh.hippocampal-fissure,rh.molecular_layer_HP-body,rh.molecular_layer_HP-head,rh.parasubiculum,rh.presubiculum-body,rh.presubiculum-head,rh.subiculum-body,rh.subiculum-head,whole.CA1-body,whole.CA1-head,whole.CA3-body,whole.CA3-head,whole.CA4-body,whole.CA4-head,whole.GC-ML-DG-body,whole.GC-ML-DG-head,whole.HATA,whole.Hippocampal_tail,whole.Whole_hippocampal_body,whole.Whole_hippocampal_head,whole.Whole_hippocampus,whole.fimbria,whole.hippocampal-fissure,whole.molecular_layer_HP-body,whole.molecular_layer_HP-head,whole.parasubiculum,whole.presubiculum-body,whole.presubiculum-head,whole.subiculum-body,whole.subiculum-head
mriface_0001,160.140423,643.111512,121.733859,181.450963,136.089629,176.885647,142.768001,214.595877,75.325365,675.274931,1293.282294,2090.659074,4059.216299,74.752421,188.162283,249.214902,399.634974,65.916415,147.868612,140.671683,260.714447,193.066639,160.417880,557.584677,144.919215,150.610933,158.435914,146.932857,167.188920,174.063578,58.895612,646.290298,1294.041853,1785.289454,3725.621605,64.672166,186.310011,244.291453,344.342486,61.677945,112.838806,121.130778,241.277501,170.050586,320.558303,1200.696189,266.653074,332.061896,294.525543,323.818504,309.956921,388.659455,134.220977,1321.565229,2587.324147,3875.948528,7784.837904,139.424587,374.472294,493.506355,743.97746,127.59436,260.707418,261.802461,501.991948,363.117225
mriface_0002,103.287906,482.663069,70.404188,136.628715,108.229725,130.864657,123.395396,154.621026,59.824302,483.804995,1049.338826,1653.070021,3186.213842,77.842256,135.308845,188.309856,312.397416,75.537688,143.016055,140.327036,234.853444,160.206112,94.023260,490.761234,91.329668,133.129240,126.997684,132.734252,141.050346,152.248955,60.223939,448.709529,1104.744740,1655.016815,3208.471084,61.685305,166.993028,202.158720,317.133128,73.087146,146.284817,132.921255,241.214940,162.777666,197.311166,973.424303,161.733856,269.757955,235.227409,263.598909,264.445742,306.869981,120.048241,932.514524,2154.083566,3308.086836,6394.684926,139.527561,302.301873,390.468576,629.530544,148.624834,289.300872,273.248291,476.068384,322.983778
mriface_0003,93.472053,292.211426,63.399442,75.138635,83.350543,79.036618,89.423755,86.168844,28.511534,351.012237,744.342056,991.830003,2087.184295,27.602650,125.530329,142.094642,188.495860,53.322443,103.487495,80.480452,141.511475,108.464190,106.828836,342.747996,71.336504,75.118041,90.448081,86.619239,95.673876,95.907348,32.929040,404.474760,846.803326,1187.683959,2438.962045,30.090636,164.192786,160.769831,228.038129,68.750950,119.683304,102.904924,171.972258,154.668293,200.300889,634.959422,134.735946,150.256676,173.798624,165.655857,185.097631,182.076192,61.440574,755.486997,1591.145382,2179.513962,4526.14634,57.693286,289.723115,302.864473,416.533989,122.073393,223.170799,183.385376,313.483733,263.132483
mriface_0004,90.025085,425.035298,61.411596,95.384809,85.516784,100.256754,93.233374,112.652146,55.431434,376.308526,872.067514,1349.532440,2597.908481,41.677259,159.021456,159.268126,242.911788,69.730972,134.091369,110.008671,206.843921,138.120570,110.674897,370.394924,62.026268,95.396101,92.756480,97.983884,102.218936,105.636656,52.204311,385.544462,849.524182,1210.263846,2445.332491,34.370706,176.590676,159.272561,218.078316,58.590188,103.594471,85.746939,184.609864,126.232527,200.699982,795.430222,123.437864,190.78091,178.273264,198.240638,195.45231,218.288802,107.635745,761.852988,1721.591696,2559.796286,5043.240972,76.047965,335.612132,318.540687,460.990104,128.32116,237.68584,195.75561,391.453785,264.353097
mriface_0005,93.859758,324.806652,77.916251,83.807912,89.245013,87.995024,94.582035,96.679923,49.024978,414.336672,810.804979,1157.485099,2382.626750,37.909804,141.689681,154.416912,208.037965,68.748935,98.145165,107.828875,164.730041,130.554836,119.825141,406.823919,85.275754,95.169639,94.333964,107.044569,102.725566,121.103714,48.089261,506.876257,921.452891,1405.736164,2834.065313,37.298695,155.083729,182.466672,265.513688,81.957020,102.535701,121.564807,196.991399,158.469547,213.684899,731.630571,163.192005,178.977551,183.578977,195.039593,197.307601,217.783637,97.114239,921.212929,1732.25787,2563.221263,5216.692063,75.208499,296.77341,336.883584,473.551653,150.705955,200.680866,229.393682,361.72144,289.024383

Preprocesamiento SBM

El analisis SBM debe prepararse y hacerse independientemente del flujo de proyecto ya que requiere preparar una DB independiente con datos externos a neuroimagen. No obstante, pueden preprocesarse las MRI dentro de la estructura general del proyecto.

Para esto basta con ejecutar un recon_all -qcache -subjid project_subjectid. Asi que es tan simple como recorrer la DB y ejecutar esto en paralelo con,

$ pqcache.pl proyecto

y ahi quedan todas las medias preparadas para distintos valores de FWHM,

[osotolongo@brick03 ~]$ ls /old_nas/subjects/bioface_0001/surf/*fwhm*
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm0.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/lh.volume.fwhm25.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm10.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.volume.fwhm5.fsaverage.mgh           /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm15.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm0.fsaverage.mgh      /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm20.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm10.fsaverage.mgh     /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm25.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm15.fsaverage.mgh     /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.fwhm5.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm20.fsaverage.mgh     /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm0.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm25.fsaverage.mgh     /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm10.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/lh.w-g.pct.mgh.fwhm5.fsaverage.mgh      /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm15.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm0.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm20.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm10.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.sulc.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm25.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm15.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.area.pial.fwhm5.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm20.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm0.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm25.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm10.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.white.H.fwhm5.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm15.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm0.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm20.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm10.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.thickness.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm25.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm15.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.curv.fwhm5.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm20.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm0.fsaverage.mgh   /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm25.fsaverage.mgh         /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm10.fsaverage.mgh  /old_nas/subjects/bioface_0001/surf/lh.white.K.fwhm5.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm15.fsaverage.mgh  /old_nas/subjects/bioface_0001/surf/rh.area.fwhm0.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm20.fsaverage.mgh  /old_nas/subjects/bioface_0001/surf/rh.area.fwhm10.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.volume.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm25.fsaverage.mgh  /old_nas/subjects/bioface_0001/surf/rh.area.fwhm15.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.jacobian_white.fwhm5.fsaverage.mgh   /old_nas/subjects/bioface_0001/surf/rh.area.fwhm20.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm0.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.area.fwhm25.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm10.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.area.fwhm5.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm15.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm0.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm20.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm10.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.w-g.pct.mgh.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm25.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm15.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.sulc.fwhm5.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm20.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm0.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm25.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm10.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.area.pial.fwhm5.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm15.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm0.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm20.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm10.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.white.H.fwhm5.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm25.fsaverage.mgh       /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm15.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm0.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.thickness.fwhm5.fsaverage.mgh        /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm20.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm10.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.volume.fwhm0.fsaverage.mgh           /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm25.fsaverage.mgh            /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm15.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.volume.fwhm10.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.curv.fwhm5.fsaverage.mgh             /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm20.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.volume.fwhm15.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm0.fsaverage.mgh   /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm25.fsaverage.mgh
/old_nas/subjects/bioface_0001/surf/lh.volume.fwhm20.fsaverage.mgh          /old_nas/subjects/bioface_0001/surf/rh.jacobian_white.fwhm10.fsaverage.mgh  /old_nas/subjects/bioface_0001/surf/rh.white.K.fwhm5.fsaverage.mgh

Ver https://github.com/asqwerty666/acenip/blob/main/bin/pqcache.pl

DTI preproc

El procesamiento no cambia mucho del que se realiza en la version anterior. No obstante, al tomar los archivos del formato BIDS, la estructura en que se guarda es por completo abstracta del software y la comprobaicon de los archivos (y posterior toma de decisiones) se puede individualizar en una subrutina externa.

Nota: deben existirlos archivos acqparams.txt y dti_index.txt. Ver como hacer el index.txt y el acqparams.txt

El registro es ahora un poco mas sencillo

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 NEURO4 qw(get_subjects check_subj load_project print_help get_list check_or_make cut_shit);
 
my $study;
my $cfile;
my $debug=0;
my $old=0;
my $chop=0;
my $timeout='8:0:0';
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    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_project($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_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;
 
my $w_dir=$std{'WORKING'};
my $bids_dir=$std{'BIDS'};
my $db = $std{'DATA'}.'/'.$study.'_mri.csv';
my $data_dir=$std{'DATA'};
 
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
 
print "Collecting needed files\n";
 
my @dtis = cut_shit($db, $data_dir."/".$cfile);
 
foreach my $subject (sort @dtis){
	if($subject){
		my %nifti = check_subj($std{'DATA'},$subject);
		if($nifti{'T1w'} && $nifti{'dwi'}){
			my $order;
			if($old){
				if($chop){
					$order = $pipe_dir."/bin/dti_proc_x.sh ".$study." ".$subject." ".$nifti{'dwi'}." ".$nifti{'T1w'}." ".$w_dir;
				}else{
					$order = $pipe_dir."/bin/dti_proc_deprecated.sh ".$subject." ".$nifti{'dwi'}." ".$w_dir;
				}
			}else{
				if($nifti{'dwi_sbref'}){
       		                	$order = $pipe_dir."/bin/dti_proc_epi.sh ".$study." ".$subject." ".$nifti{'dwi'}." ".$nifti{'dwi_sbref'}." ".$nifti{'T1w'}." ".$w_dir;
       		                }else{
                                	$order = $pipe_dir."/bin/dti_proc_uncorr.sh ".$study." ".$subject." ".$nifti{'dwi'}." ".$nifti{'T1w'}." ".$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='.$timeout."\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($old){ print ORD '#SBATCH --gres=gpu:1'."\n"; }
			print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
			unless($old){ print ORD '#SBATCH -p cuda'."\n"; }else{print ORD '#SBATCH -p fast'."\n";}
			print ORD "srun $order\n";
			close ORD;
			system("sbatch $orderfile");
			$debug ? print DBG "$order\n" :0;
		}
	}
}
$debug ? close DBG:0;  
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 --dependency=singleton '.$orderfile;
exec($order);

Los procedimientos individuales de registro si han de cambiarse poco, mas que nada el input, lo demas es estrictamente igual

con epi_reg

dt_proc_epi.sh

#!/bin/sh 
 
Usage() {  
    echo "" 
    echo "Usage: dti_proc_epi.sh  <study> <subject> <dwi> <dwi_pa> <T1> <output_dir>"  
    echo ""  
    echo "You must have FSL installed in order to run this script" 
    echo "" 
    exit 1 
} 
 
[ "$6" = "" ] && Usage 
debug=1
study=$1
shift
b_in=$1
shift  
in=`${FSLDIR}/bin/remove_ext $1` 
shift
p2a=`${FSLDIR}/bin/remove_ext $1`
shift  
t1=`${FSLDIR}/bin/remove_ext $1` 
shift  
out=$1 
shift  
td=${out}'/.tmp_'${b_in}  
if [ ! -d "$out" ]; then
        mkdir $out
fi
if [ ! -d "$td" ]; then
        mkdir $td
fi
echo "DTI preproccessing begins on ${b_in} ..." 
echo [`date`] 
echo 
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1
echo "Preparing topup"
${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${p2a}
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_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
c3d_affine_tool -ref ${td}/hifi_b0.nii.gz -src ${out}/${b_in}_t1_reoriented.nii.gz ${td}/${b_in}_tmp_std2diff.mat -fsl2ras -oitk ${td}/${b_in}_epi_reg_ANTS.mat
antsApplyTransforms -d 3 -i ${out}/${b_in}_t1_reoriented.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_t1_b0.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_mni_t1_warped.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_mni_to_b0.nii.gz
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
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_labels_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_labels.nii.gz
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
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_tracts_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_tracts.nii.gz
${FSLDIR}/bin/imcp ${td}/hifi_b0 ${out}/${b_in}_hifi_b0
echo "Game over :-P"
echo [`date`]
exit 0

[tester@detritus ptest]$ dti_reg.pl ptest
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114225      cuda dti_reg_   tester PD       0:00      1 (Resources)
            114226      cuda dti_reg_   tester PD       0:00      1 (Priority)
            114227      cuda dti_reg_   tester PD       0:00      1 (Priority)
            114228      cuda dti_reg_   tester PD       0:00      1 (Priority)
            114229      cuda dti_reg_   tester PD       0:00      1 (Priority)
            114230      cuda dti_reg_   tester PD       0:00      1 (Priority)
            114231     devel dti_reg_   tester PD       0:00      1 (Dependency)
            114223      cuda dti_reg_   tester  R       0:05      1 brick01
            114224      cuda dti_reg_   tester  R       0:05      1 brick01
[tester@detritus ptest]$ ls working/0001_*
working/0001_brain.nii.gz           working/0001_dti_L1.nii.gz  working/0001_dti_MO.nii.gz  working/0001_dti_V3.nii.gz      working/0001_mni_to_b0.nii.gz
working/0001_dti_brain_mask.nii.gz  working/0001_dti_L2.nii.gz  working/0001_dti_S0.nii.gz  working/0001_hifi_b0.nii.gz     working/0001_t1_b0.nii.gz
working/0001_dti_data.nii.gz        working/0001_dti_L3.nii.gz  working/0001_dti_V1.nii.gz  working/0001_JHU_labels.nii.gz  working/0001_t1_reoriented.nii.gz
working/0001_dti_FA.nii.gz          working/0001_dti_MD.nii.gz  working/0001_dti_V2.nii.gz  working/0001_JHU_tracts.nii.gz

Probando el -cut (que esta reescrito completo),

[tester@detritus ptest]$ cat repdti.txt 
0004
[tester@detritus ptest]$ dti_reg.pl -cut repdti.txt ptest
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114233     devel dti_reg_   tester PD       0:00      1 (Dependency)
            114232      cuda dti_reg_   tester  R    6:25:12      1 brick01

DTI metrics (Atlas)

Esto tampoco cambia demasiado, pero se aprovechan las subs nuevas y se simplifica la ubicacion de los archivos

dti_metrics.pl

#!/usr/bin/perl
 
use strict; use warnings;
use File::Find::Rule;
use NEURO qw(get_subjects print_help get_pair get_list load_project shit_done);
use Data::Dump qw(dump);
use File::Basename qw(basename);
 
my $with_sd = 0;
my $atlas = 2;
my $cfile;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-a1/) { $atlas = 1;}
    if (/^-a2/) { $atlas = 2;}
    if (/^-cut/) { $cfile = shift; chomp($cfile);} 
    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_project($study);
 
my $w_dir=$std{'WORKING'};
my $data_dir=$std{'DATA'};
my $db = $std{'DATA'}.'/'.$study.'_mri.csv';
my $roi_list;
my $roi_group;
if($atlas==1){
	$roi_list = $ENV{'PIPEDIR'}.'/lib/jhu.labels.list';
	$roi_group = "labels";
}elsif($atlas==2){
	$roi_list = $ENV{'PIPEDIR'}.'/lib/jhu.tracts.list';
	$roi_group = "tracts";
}
my $roi_dir = $ENV{'PIPEDIR'}.'/lib/jhu/';
my $jhu = "JHU_".$roi_group;
my %masks = get_pair $roi_list;
print "Collecting needed files\n";
 
my @dtis = cut_shit($db, $data_dir."/".$cfile);
 
my %csv;
 
my $ofile = $data_dir."/".$study."_dti_".$roi_group.".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",";$masks{$rmask}","_MD_Mean;","$masks{$rmask}","_MD_STD";
	}else{
		print OF ";$masks{$rmask}","_FA",";$masks{$rmask}","_MD";
	}
}
print OF "\n";
 
foreach my $subject (@dtis){
        my $dti_fa = $w_dir.'/'.$subject.'_dti_FA.nii.gz';
        my $dti_md = $w_dir.'/'.$subject.'_dti_MD.nii.gz';	
	if($subject){
		foreach my $rmask (sort keys %masks){
			my ($m_index) = $rmask =~ /.*_(\d+)$/;
			if($m_index){
				my $order = "fslmaths ".$w_dir."/".$subject."_".$jhu." -uthr ".$m_index." -thr ".$m_index." -div ".$m_index." ".$w_dir."/.tmp_".$subject."/JHU_".$rmask."_tmp";
			print "$order\n";
			system($order);
			} else {
			my $order = "fslmaths ".$w_dir."/".$subject."_".$jhu." -uthr ".$m_index." -thr ".$m_index." ".$w_dir."/.tmp_".$subject."/JHU_".$rmask."_tmp";               
                        print "$order\n";
                        system($order);
			}
		}
		print OF "$subject";
		foreach my $rmask (sort keys %masks){
			my $order = "fslstats ".$dti_fa." -k ".$w_dir."/.tmp_".$subject."/JHU_".$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";
			}
                        $order = "fslstats ".$dti_md." -k ".$w_dir."/.tmp_".$subject."/JHU_".$rmask."_tmp -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=$std{DATA}."/".$study."_dti_".$roi_group."_results.tgz";
system("tar czf $zfile $ofile");
shit_done basename($ENV{_}), $study, $zfile;

[tester@detritus ptest]$ dti_metrics.pl -a1 ptest
[tester@detritus ptest]$ dti_metrics.pl -a2 ptest
[tester@detritus ptest]$ ls ptest_dti_*
ptest_dti_labels.csv  ptest_dti_labels_results.tgz  ptest_dti_tracts.csv  ptest_dti_tracts_results.tgz

TRACULA

La implementacion de Tracula en FS7.2 ha cambiado para pemitir la correcion por topup. Asi que hay que hacer un dmri.rc distinto. Entonces, ahora hay que correr 7 scripts,

ctrac_dmri.pl –> Crear el dmri.rc
ctrac_prep.pl –> Preprocesamiento DTI
ctrac_bedp.pl –> bedtrack
ctrac_path.pl –> probtrack
ctrac_stat.pl –> Estadisticas globales
ctrac_metrics.pl –> Metricas medias de los principales tractos
ctrac_report.pl –> QC report

DTI metrics VOIs

Las metricas de DTI pueden extraerse en cualquiera de los tractos definidos, ya sea mediante atlas o a traves de TRACULA. Lo que ocurre es que en la realidad estos tractos son demasiado específicos. Es decir, no nos interesa la FA en el Genu of corpus callosum sino mas bien en el Corpus callosum en su totalidad. Hay dos enfoques validos para conseguir esto ultimo.

ICBM VOIs

Para calcular una mascara de los tractos completos (VOI) basta con tomar las mascaras de los trozos individuales y sumarlos todos juntos en una mscara binaria. Esta mascara es la que luego se aplicara a la imagen FA y se calculara el valor medio dentro.

Para esto necesitamos saber que trozos componen cada VOI. Esto lo voy a almacenar en un archivo de texto ya que no cambiara nunca. Ejemplo, para el atlas ICBM-DTI-81,

Corpus_callosum:3,4,5
Internal_capsule:17,18,19,20,21,22
Corona_radiata:23,24,25,26,27,28
Superior_longitudinal_fasciculus:41,42
Fornix:6,39,40
Cingulum:35,36,37,38
Uncinate_fasciculus:45,46

Este archivo lo he de leer y tomar cada archivo XXXX_JHU_labels individual y calcular ahi la porcion con la que he de fabricar la mascara. Luego se suman y se saca el valor deseado.

foreach my $subject (@dtis){
        my $dti_fa = $w_dir.'/'.$subject.'_dti_FA.nii.gz';
        my $dti_md = $w_dir.'/'.$subject.'_dti_MD.nii.gz';
        if($subject && -e $dti_fa){
                foreach my $rmask (sort keys %masks){
                        my @chunks = split /,/,$masks{$rmask};
                        my @orarg;
                        foreach my $chunk (@chunks){
                                my $order = "fslmaths ".$w_dir."/".$subject."_JHU_".$lort{$atlas}." -uthr ".$chunk." -thr ".$chunk." -div ".$chunk." ".$w_dir."/.tmp_".$subject."/JHU_".$chunk."_tmp";
                                print "$order\n";
                                system($order);
                                my $tmpstr = $w_dir."/.tmp_".$subject."/JHU_".$chunk."_tmp";
                                push @orarg, $tmpstr;
                        }
                        my $strarg = join " -add ", @orarg;
                        my $order = "fslmaths ".$strarg." ".$w_dir."/.tmp_".$subject."/JHU_".$rmask."_tmp";
                        print "$order\n";
                        system($order);
                }
                print OF "$subject";
                foreach my $rmask (sort keys %masks){
                        my $order = "fslstats ".$dti_fa." -k ".$w_dir."/.tmp_".$subject."/JHU_".$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";
                        }
                        $order = "fslstats ".$dti_md." -k ".$w_dir."/.tmp_".$subject."/JHU_".$rmask."_tmp -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;

El codigo completo puede encontrarse en github y

para ejecutarlo basta con hacer,

$ dti_metrics_comp.pl project

o,

$ dti_metrics_comp.pl -a1 project

para hacerlo con el atlas JHU-WM-tractography (on development).

En el primer caso los resultados quedarian en el archivo bioface_dti_icbm_comp.csv,

[osotolongo@brick03 bioface]$ head bioface_dti_icbm_comp.csv
Subject;Cingulum_FA;Cingulum_MD;Corona_radiata_FA;Corona_radiata_MD;Corpus_callosum_FA;Corpus_callosum_MD;Fornix_FA;Fornix_MD;Internal_capsule_FA;Internal_capsule_MD;Superior_longitudinal_fasciculus_FA;Superior_longitudinal_fasciculus_MD;Uncinate_fasciculus_FA;Uncinate_fasciculus_MD
0002;0.477487;0.000786;0.433740;0.000744;0.647378;0.000841;0.554593;0.000944;0.556976;0.000720;0.474224;0.000726;0.566348;0.000788
0003;0.494214;0.000800;0.439710;0.000787;0.616569;0.000850;0.577302;0.000925;0.560079;0.000730;0.425376;0.000787;0.431516;0.000902
0004;0.416937;0.000813;0.398660;0.000779;0.595173;0.000903;0.490137;0.001118;0.559346;0.000735;0.410975;0.000748;0.468427;0.000809
0005;0.427343;0.000794;0.404245;0.000765;0.591447;0.000888;0.460360;0.001119;0.563446;0.000757;0.426342;0.000773;0.510532;0.000755
0006;0.475403;0.000801;0.431013;0.000781;0.596390;0.000889;0.521160;0.001085;0.560720;0.000742;0.474176;0.000756;0.471718;0.000817
0007;0.474952;0.000768;0.460331;0.000757;0.575729;0.000977;0.407189;0.001516;0.581289;0.000731;0.465529;0.000735;0.521556;0.000757
0008;0.437555;0.000801;0.416142;0.000778;0.596577;0.000866;0.497081;0.001035;0.554840;0.000767;0.426218;0.000771;0.518556;0.000786
0009;0.442004;0.000812;0.431805;0.000747;0.630885;0.000814;0.499554;0.001055;0.539465;0.000748;0.470043;0.000736;0.537202;0.000756
0010;0.473434;0.000748;0.447876;0.000740;0.626619;0.000852;0.511838;0.001002;0.601466;0.000712;0.470705;0.000716;0.529629;0.000725

TRACULA VOIs

Los valores de FA y MD de los tractos completos tambien puede derivarse de TRACULA. En este caso no reconstruiremos el tracto como tal sino que tomaremos los valores de FA individuales y asumiremos que,

$$ {FA}_{VOI} \approx \frac{\sum{{FA}_{i}*V_{i}}}{\sum{V_{i}}} $$

donde ${FA}_{i}$ son los valores de FA y $V_{i}$ el volumen, de los tractos individuales, respectivamente.

Asi que primero tengo que tener almacenado en un fichero que variables forman cada VOI,

Fornix:lh.fx,rh.fx
Cingulum:lh.cbd,lh.cbv,rh.cbd,rh.cbv
Uncinate_fasciculus:lh.uf,rh.uf
Superior_longitudinal_fasciculus:lh.slf1,lh.slf2,lh.slf3,rh.slf1,rh.slf2,rh.slf3
Inferior_longitudinal_fasciculus:lh.ilf,rh.ilf
Corpus_callosum:cc.bodyc,cc.bodyp,cc.bodypf,cc.bodypm,cc.bodyt,cc.genu,cc.rostrum,cc.splenium

y esto lo leo directamente de los resultados de tracula para cada individuo. Podria leerlo de los achivos de FS pero lo que voy a hacer es,

Cambiar el script de metricas de TRACULA para incluir el volumen de los tractos
Asumir que he ejecutado las metrcas de TRACULA y leer directamente este archivo para calcular los resultados (mucho mas rapido)

El resultado de esto es el script en github que se ejecuta como,

[osotolongo@brick03 bioface]$ ctrac_translate.pl bioface

y deja los resultados en bioface_dti_tracula_composed.csv,

Subject,Cingulum_FA,Cingulum_MD,Corpus_callosum_FA,Corpus_callosum_MD,Fornix_FA,Fornix_MD,Inferior_longitudinal_fasciculus_FA,Inferior_longitudinal_fasciculus_MD,Superior_longitudinal_fasciculus_FA,Superior_longitudinal_fasciculus_MD,Uncinate_fasciculus_FA,Uncinate_fasciculus_MD
0001,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
0002,0.500092918105356,0.000781457339530766,0.519852793601731,0.000781114487752106,0.466710179411765,0.00105999035294118,0.504764250681199,0.000813264,0.436700603806417,0.000742181515280044,0.468751002475248,0.000772750495874587
0003,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
0004,0.471647589821662,0.000811000234449761,0.512361431638802,0.00081236424127981,0.386831987261146,0.000992012576433121,0.452685824742268,0.000871269313402062,0.417407322635963,0.00077327638510534,0.419668255865922,0.000848094241340782
0005,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
0006,0.486058795640327,0.000798137305177112,0.52619407984211,0.000808463441643492,0.449970944306931,0.000960078674504951,0.496323946902655,0.000859897426548673,0.44573164370485,0.000777395058565531,0.463258846660395,0.000819830761994355
0007,0.476584601431981,0.000773005742243437,0.525463044190302,0.000778800766056725,0.308487,0.00163514,0.517803210059172,0.000801609260946746,0.451237894539822,0.00074905124422357,0.467696353195164,0.000769760035405872
0008,0.470472873369794,0.000793977488851493,0.501802720765106,0.000807642284451529,0.433634495232041,0.00100513844373808,0.461641992471769,0.000858071812630698,0.402474214795654,0.000775522433212623,0.442492631840796,0.00079464460199005
0009,0.513032479578393,0.000781708291172596,0.524635845259849,0.000773339304509904,0.437998623059867,0.00105553261640798,0.479583220908231,0.000831665455061495,0.447428453367876,0.00074106617074894,0.461369545135846,0.00078089697283085

Comparando ICBM vs TRACULA

DTI Tractography

La tractografia puede hacerse utilizando regiones definidas en la reconstruccion de FS o un atlas externo. En este caso el atlas externo escogido es el de la UofM.

Voy a simplificar un poco el procesamiento anterior porque hay muchas cosas que no se usan. En caso de necesitarlas ya se volveran a implementar.

Primero el bedpostx y el probtrackx.

Los dos modos de procesamiento que hay implementados tienen identico resultado. Basta con dejar uno (el mas rapido o sencillo) y nos quitamos una variable de encima.
Pasamos el registro de ANTS a la nueva sintaxis, ya que hemos cambiado el registro y las matrices son otras

Queda bastante sencillo ahora

dti_bedtrack_cuda.sh

#!/bin/sh
 
Usage() {                      
    echo ""                                                                                                                                                              
    echo "Usage: dti_bedtrack.sh  <project> <subject> <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}'/../dti_track.seed'
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_aparc.sh ${prj} ${pollo} ${w_dir}
antsApplyTransforms -d 3 -i ${w_dir}/${pollo}_aseg.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${pollo}_epi_reg_ANTS.mat -o ${td}/${pollo}_aseg_warped.nii.gz
#WarpImageMultiTransform 3 ${w_dir}/${pollo}_aseg.nii.gz ${td}/${pollo}_aseg_warped.nii.gz -R ${td}/hifi_b0.nii.gz ${td}/${pollo}_dti_t1_b0Warp.nii.gz ${td}/${pollo}_dti_t1_b0Affine.txt
###########################################
######## Falta calcular las mascaras ######
###########################################
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;
###########################################		
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;
mv ${td}/probtrack_out ${w_dir}/${pollo}_probtrack_out
echo "Done"
echo [`date`]

Vamos a probar uno,

[tester@brick01 ptest]$ cat dti_track.seed
10
49
[tester@brick01 ptest]$ dti_bedtrack_cuda.sh ptest 0001 /nas/data/ptest/working
Copying files
Making bedpostx
[Tue Dec 10 10:36:45 CET 2019]
---------------------------------------------
------------ BedpostX GPU Version -----------
---------------------------------------------
subjectdir is /nas/data/ptest/working/.tmp_0001/bedpostx
Making bedpostx directory structure
Copying files to bedpost directory
Pre-processing stage
Queuing parallel processing stage
 
----- Bedpostx Monitor -----
1 parts processed out of 4
2 parts processed out of 4
3 parts processed out of 4
.....
.....
.....
Done loading samples.
Volume seeds
volume 1
volume 2
 
time spent tracking: 1708 seconds
 
save results
finished
Done
[Tue Dec 10 11:39:22 CET 2019]
[tester@brick01 ptest]$ ls /nas/data/ptest/working/0001_probtrack_out/
fdt_paths.nii.gz  probtrackx.log  waytotal

Esto es eterno, asi que paciencia.

Ahora con el UofM,

[tester@brick01 ptest]$ dti_bedtrack_nodes.sh ptest 0001 /nas/data/ptest/working/ /nas/software/neuro_4_0/lib/uofm/Nodes/LN/
Copying files
Making bedpostx
[Tue Dec 10 15:53:58 CET 2019]
---------------------------------------------
------------ BedpostX GPU Version -----------
---------------------------------------------
subjectdir is /nas/data/ptest/working/.tmp_0001/bedpostx
/nas/data/ptest/working/.tmp_0001/bedpostx has already been processed: /nas/data/ptest/working/.tmp_0001/bedpostx.bedpostX.
Delete or rename /nas/data/ptest/working/.tmp_0001/bedpostx.bedpostX before repeating the process.
So far, so good
[Tue Dec 10 15:53:58 CET 2019]
Getting nodes and making masks
.....
.....
.....
Done loading samples.
Volume seeds
volume 1
volume 2
volume 3
volume 4
volume 5
volume 6
volume 7
 
time spent tracking: 3714 seconds
 
save results
finished
Done
[Tue Dec 10 17:03:25 CET 2019]
[tester@brick01 ptest]$ ls /nas/data/ptest/working/0001_probtrack_out/
fdt_paths.nii.gz  probtrackx.log  waytotal

Nota: Ya tengo el bedpostx hecho asi que demorara menos, pero igual son varios nodos a registrar al DTI.

Parece que todo funciona asi que nos vamos al 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 NEURO4 qw(load_project print_help cut_shit check_or_make);
 
my $study;
my $cfile;
my $debug=1;
my $atlas=0; 
my $network;
my $t1=0;
my $time = '8:0:0';
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
    $_ = shift;
    last if /^--$/;
    if (/^-cut/) { $cfile = shift; chomp($cfile);}
    if (/^-t1/) {$t1 = 1;}
    if (/^-time/) {$time = shift; chomp($time)}
    if (/^-uofm/) {$atlas="uofm"; $network = shift; chomp($network); $network =~ s/\_/\//g;}
    if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/dti_track.hlp'; exit;}
}
$study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/dti_track.hlp'; exit;}
 
my %std = load_project($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;
 
my $w_dir=$std{'WORKING'};
my $data_dir=$std{'DATA'};
my $db = $data_dir.'/'.$study.'_mri.csv';
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
 
my @dtis = cut_shit($db, $data_dir.'/'.$cfile);
foreach my $subject (@dtis){
	#Compulsory: el DTI debe haber pasado el preproc
	my $dti_fa = $w_dir.'/'.$subject.'_dti_FA.nii.gz';
	#Compulsory: Necesito tambien la segmentacion de freesurfer de este pollo
	my $aseg = $subj_dir.'/'.$study.'_'.$subject.'/mri/aseg.mgz';
	if((-e $dti_fa) && (-e $aseg)){
		my $order;
		unless($atlas){
			if($t1){
				$order = $pipe_dir."/bin/dti_bedtrack_cuda_t1.sh ".$study." ".$subject." ".$w_dir;
			}else{
				$order = $pipe_dir."/bin/dti_bedtrack_cuda.sh ".$study." ".$subject." ".$w_dir;
			}
		}else{
			if($t1){
				$order = $pipe_dir."/bin/dti_bedtrack_nodes_t1.sh ".$study." ".$subject." ".$w_dir." ".$pipe_dir."/lib/".$atlas."/Nodes/".$network;
			}else{
				$order = $pipe_dir."/bin/dti_bedtrack_nodes.sh ".$study." ".$subject." ".$w_dir." ".$pipe_dir."/lib/".$atlas."/Nodes/".$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_'.$study."\n";
		print ORD '#SBATCH --time='.$time."\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 --gres=gpu:1'."\n";
		print ORD '#SBATCH -o '.$outdir.'/dti_track-'.$subject.'-'.$study.'-%j'."\n";
		print ORD '#SBATCH -p cuda'."\n";
		print ORD "srun $order\n";
        	close ORD;
		system("sbatch $orderfile")
		$debug ? print DBG "$order\n" :0;
	}
}
 
my $orderfile = $outdir.'/dti_track_end.sh';
open ORD, ">$orderfile";
print ORD '#!/bin/bash'."\n";
print ORD '#SBATCH -J dti_track_'.$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_track_end-%j'."\n";
print ORD ":\n";
close ORD;
my $order = 'sbatch --dependency=singleton '.$orderfile;
exec($order);

[tester@detritus ptest]$ dti_track.pl ptest
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114236      cuda dti_trac   tester PD       0:00      1 (Resources)
            114237      cuda dti_trac   tester PD       0:00      1 (Priority)
            114238      cuda dti_trac   tester PD       0:00      1 (Priority)
            114239      cuda dti_trac   tester PD       0:00      1 (Priority)
            114240      cuda dti_trac   tester PD       0:00      1 (Priority)
            114241     devel dti_trac   tester PD       0:00      1 (Dependency)
            114234      cuda dti_trac   tester  R       0:03      1 brick01
            114235      cuda dti_trac   tester  R       0:03      1 brick01

Importante: : Al terminar, mover los directorios al path correcto,ejemplo,

for x in `ls -d working/*_probtrack_out`; do mv $x `echo $x | sed 's/out/DMN/'`;done

ICA individual

El procedimiento para ICA de un sujeto fue descrito ampliamente usando el protocolo MOPEAD. Aqui hay que cambiar el input y alguna cosa mas,

Todo muy basico

rs_ica_one.pl

#!/usr/bin/perl
 
# Copyright 2019 O. Sotolongo <asqwerty@gmail.com>
 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
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 NEURO4 qw(check_subj load_project print_help cut_shit 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_project($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
my $fsl = $ENV{'FSLDIR'};
 
my $w_dir=$std{'WORKING'};
my $db = $std{'DATA'}.'/'.$study.'_mri.csv';
my $data_dir=$std{'DATA'};
my $template_file = $pipe_dir.'/lib/fsf/lone_ica_template.fsf';
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
print "Collecting needed files\n";
my @fmris = cut_shit($db, $data_dir.'/'.$cfile);
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 @fmris){
	my %nifti = check_subj($std{'DATA'},$subject);
	if($nifti{'func'}){
		my $output_dir = $w_dir.'/'.$subject.'_rsout';
		check_or_make($output_dir);
		my $data = $w_dir.'/'.$subject.'_rs';
		system("$fsl'/bin/imcp' $nifti{'func'} $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.'/'.$subject.'_design.fsf';
		open TPF,"<$template_file";
		open DSG,">$dsg_file";
		while(<TPF>){
			s/<study>/$study/;
			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 = $subject;
			while(<TPF>){
				s/<study>/$study/;
				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="";
	}
}
my $orderfile = $outdir.'/feat_end.sh';
open ORD, ">$orderfile";
print ORD '#!/bin/bash'."\n";
print ORD '#SBATCH -J feat-'.$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.'/feat_end-%j'."\n";
print ORD ":\n";
close ORD;
my $order = 'sbatch --dependency=singleton '.$orderfile;
exec($order);

Ahora a probarlo,

[tester@detritus ptest]$ rs_ica_one.pl ptest
Collecting needed files  
sbatch /nas/data/ptest/working/0001_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114303 /nas/data/ptest/working/0001_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114304 /nas/data/ptest/working/0001_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0002_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114306 /nas/data/ptest/working/0002_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114307 /nas/data/ptest/working/0002_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0003_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114309 /nas/data/ptest/working/0003_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114310 /nas/data/ptest/working/0003_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0004_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114312 /nas/data/ptest/working/0004_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114313 /nas/data/ptest/working/0004_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0005_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114315 /nas/data/ptest/working/0005_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114316 /nas/data/ptest/working/0005_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0006_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114318 /nas/data/ptest/working/0006_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114319 /nas/data/ptest/working/0006_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0007_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114321 /nas/data/ptest/working/0007_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114322 /nas/data/ptest/working/0007_rs.ica/scripts/feat4.sh
sbatch /nas/data/ptest/working/0008_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114324 /nas/data/ptest/working/0008_rs.ica/scripts/feat2.sh
sbatch --depend=afterok:114325 /nas/data/ptest/working/0008_rs.ica/scripts/feat4.sh
Submitted batch job 114327
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114304     devel feat-pte   tester PD       0:00      1 (Dependency)
            114305     devel feat-pte   tester PD       0:00      1 (Dependency)
            114307     devel feat-pte   tester PD       0:00      1 (Dependency)
            114308     devel feat-pte   tester PD       0:00      1 (Dependency)
            114310     devel feat-pte   tester PD       0:00      1 (Dependency)
            114311     devel feat-pte   tester PD       0:00      1 (Dependency)
            114313     devel feat-pte   tester PD       0:00      1 (Dependency)
            114314     devel feat-pte   tester PD       0:00      1 (Dependency)
            114316     devel feat-pte   tester PD       0:00      1 (Dependency)
            114317     devel feat-pte   tester PD       0:00      1 (Dependency)
            114319     devel feat-pte   tester PD       0:00      1 (Dependency)
            114320     devel feat-pte   tester PD       0:00      1 (Dependency)
            114322     devel feat-pte   tester PD       0:00      1 (Dependency)
            114323     devel feat-pte   tester PD       0:00      1 (Dependency)
            114325     devel feat-pte   tester PD       0:00      1 (Dependency)
            114326     devel feat-pte   tester PD       0:00      1 (Dependency)
            114327     devel feat-pte   tester PD       0:00      1 (Dependency)
            114303     devel feat-pte   tester  R       0:16      1 brick01
            114306     devel feat-pte   tester  R       0:13      1 brick01
            114309     devel feat-pte   tester  R       0:13      1 brick01
            114312     devel feat-pte   tester  R       0:10      1 brick01
            114315     devel feat-pte   tester  R       0:07      1 brick01
            114318     devel feat-pte   tester  R       0:07      1 brick01
            114321     devel feat-pte   tester  R       0:04      1 brick01
            114324     devel feat-pte   tester  R       0:01      1 brick01
 [tester@detritus ptest]$ ls -d working/*rs.ica
working/0001_rs.ica  working/0002_rs.ica  working/0003_rs.ica  working/0004_rs.ica  working/0005_rs.ica  working/0006_rs.ica  working/0007_rs.ica  working/0008_rs.ica

so far so good

ICA grupal

Este procedimiento es mucho mas complicado pero por esa misma razon esta escrito mucho mas abstracto.

Asi que apenas hay que cambiar nada

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 NEURO4 qw(cut_shit check_subj load_project print_help 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_project($study);
 
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $pipe_dir = $ENV{'PIPEDIR'};
my $fsl = $ENV{'FSLDIR'};
 
my $w_dir=$std{'WORKING'};
my $db = $std{'DATA'}.'/'.$study.'_mri.csv';
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');
my $error_data_file = $data_dir.'/fmri_image_errors.txt';
 
print "Collecting needed files\n";
my @fmris = cut_shit($db, $data_dir.'/'.$cfile);
 
my $count = 0;
my %subjects;
print "Counting available subjects\n";
foreach my $subject (sort @fmris){
	my %nifti = check_subj($std{'DATA'},$subject);
	if($nifti{'func'}){
		my $sname = $subject;
		$subjects{$sname} = $nifti{'func'}; 
		$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 "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;
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;
}
#print "Everything is fine. You pass the test\n"; dump %subjects; exit;
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/<study>/$study/;
				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/<study>/$study/;
	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);

[tester@detritus ptest]$ rs_ica_group.pl ptest
Collecting needed files
Counting available subjects
Getting info from images
Checking images and excluding wrong subjects
Copying FSL files and setting directories
Making global .fsf file
Making individual .fsf files and scripts
sbatch /nas/data/ptest/working/0001_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114328 /nas/data/ptest/working/0001_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0002_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114330 /nas/data/ptest/working/0002_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0003_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114332 /nas/data/ptest/working/0003_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0004_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114334 /nas/data/ptest/working/0004_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0005_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114336 /nas/data/ptest/working/0005_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0006_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114338 /nas/data/ptest/working/0006_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0007_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114340 /nas/data/ptest/working/0007_rs.ica/scripts/feat2.sh
sbatch /nas/data/ptest/working/0008_rs.ica/scripts/feat1.sh
sbatch --depend=afterok:114342 /nas/data/ptest/working/0008_rs.ica/scripts/feat2.sh
Making global script
sbatch --depend=afterok:114329,114331,114333,114335,114337,114339,114341,114343 /nas/data/ptest/working/rs.gica/scripts/feat4_ica.sh
Submitted batch job 114344
[tester@detritus ptest]$ squeue
             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
            114344     devel feat-pte   tester PD       0:00      1 (Dependency)
            114329     devel feat-pte   tester  R       0:05      1 brick02
            114331     devel feat-pte   tester  R       0:11      1 brick01
            114333     devel feat-pte   tester  R       0:11      1 brick01
            114335     devel feat-pte   tester  R       0:05      1 brick01
            114337     devel feat-pte   tester  R       0:05      1 brick02
            114339     devel feat-pte   tester  R       0:05      1 brick01
            114341     devel feat-pte   tester  R       0:05      1 brick02
            114343     devel feat-pte   tester  R       0:05      1 brick01

Table of Contents