Table of Contents
WMH pgs
- Se supone que un metodo dedicado a sacar las WMH es mejor
- De esta manera lo que saco es una mascara de las WMH y podria sacar los volumenes en cada tracto.
- Por lo menos no tendre que exlpicar a mas reviewers que es igual uno que otro para sacar relaciones lineales
HowTo
[osotolongo@brick03 f5cehbi]$ singularity run --cleanenv -B /nas:/nas -B /old_nas/f5cehbi/output:/output -B /old_nas/f5cehbi/input:/input /nas/usr/local/opt/singularity/pgs.simg sh /WMHs_segmentation_PGS.sh sub-0020_T1w.nii.gz sub-0020_T2w.nii.gz sub-0020_WMH.nii.gz Step 1 of 9: reading in images...Done! It took roughly 1 seconds Step 2 of 9: registration... 2a) Similarity transform... 2b) Affine transform... 2c) Resampling volume... Registration done! It took roughly 6 seconds Step 3 of 9: rough bias field correction...Done! It took roughly 6 seconds Step 4 of 9: calculating features...Done! It took roughly 28 seconds Step 5 of 9: voxel classification...Done! It took roughly 17 seconds Step 6 of 9: fitting the shape model...Done! It took roughly 6 seconds Step 7 of 9: free deformation...Done! It took roughly 1 seconds Step 8 of 9: building volume from mesh...Done! It took roughly 3 seconds Step 9 of 9: warping back to original space...Done! It took roughly 0 seconds Output ready The whole process took roughly 68 seconds Freeing memory and deleting temporary files... Traceback (most recent call last): File "/WMHs_segmentation_PGS.py", line 89, in <module> FLAIR_bet_data = FLAIR_data*T1_bet_mask_data ValueError: operands could not be broadcast together with shapes (176,252,256) (176,256,256)
[osotolongo@brick03 f5cehbi]$ fslinfo output/sub-0020_T1w.nii.gz data_type INT16 dim1 176 dim2 256 dim3 256 dim4 1 datatype 4 pixdim1 1.000000 pixdim2 0.976562 pixdim3 0.976562 pixdim4 2.200000 cal_max 0.0000 cal_min 0.0000 file_type NIFTI-1+ [osotolongo@brick03 f5cehbi]$ fslinfo output/sub-0020_T2w.nii.gz data_type INT16 dim1 176 dim2 252 dim3 256 dim4 1 datatype 4 pixdim1 1.000000 pixdim2 1.015625 pixdim3 1.015625 pixdim4 5.000000 cal_max 0.0000 cal_min 0.0000 file_type NIFTI-1+
Nota: Registrar T2 a T1 primero!
[osotolongo@brick03 f5cehbi]$ antsRegistrationSyNQuick.sh -d 3 -t a -f input/pre/sub-0020_T1w.nii.gz -m input/pre/sub-0020_T2w.nii.gz -o input/pre/t2tot1_ | tee t2tot1Output.txt ... ... ... [osotolongo@brick03 f5cehbi]$ antsApplyTransforms -d 3 -i input/pre/sub-0020_T2w.nii.gz -r input/pre/sub-0020_T1w.nii.gz -t input/pre/t2tot1_0GenericAffine.mat -o input/pre/sub-0020_T2w_resampled.nii.gz [osotolongo@brick03 f5cehbi]$ ls input/pre/ sub-0020_T1w.nii.gz sub-0020_T2w_resampled.nii.gz t2tot1_InverseWarped.nii.gz sub-0020_T2w.nii.gz t2tot1_0GenericAffine.mat t2tot1_Warped.nii.gz
Ahora si,
[osotolongo@brick03 f5cehbi]$ singularity run --cleanenv -B /nas:/nas -B /old_nas/f5cehbi/output:/output -B /old_nas/f5cehbi/input:/input /nas/usr/local/opt/singularity/pgs.simg sh /WMHs_segmentation_PGS.sh sub-0020_T1w.nii.gz sub-0020_T2w_resampled.nii.gz sub-0020_WMH.nii.gz
This is a highly massive computation, Dont take it lightly!
but results are amazing,
[osotolongo@brick05 f5cehbi]$ fslstats output/sub-0020_WMH.nii.gz -V 2353 2243.995605
Scripting
Primero hay que sacar las iamgenes de WMH para el proyecto,
https://github.com/asqwerty666/acenip/blob/main/bin/wmh.pl
y luego las metricas de cada imagen,
https://github.com/asqwerty666/acenip/blob/main/bin/wmh_metrics.pl
comparando
Para cada proyecto,
[osotolongo@brick03 facehbi]$ xnatapic list_subjects --project_id facehbi --label > xnat_subjects.list [osotolongo@brick03 facehbi]$ for x in `awk -F"," '{print $2}' xnat_subjects.list`; do e=$(xnatapic list_experiments --project_id facehbi --subject_id ${x} --modality MRI --label); echo "${x},${e}"; done > xnat_sub_exp.list [osotolongo@brick03 facehbi]$ for l in `cat xnat_sub_exp.list`; do s=$(echo ${l} | awk -F"," '{print $1}'); e=$(echo ${l} | awk -F"," '{print $2}'); mkdir -p fsresults/${s}; xnatapic get_fsresults --experiment_id ${e} --stats aseg fsresults/${s}/; done [osotolongo@brick03 facehbi]$ for x in `ls fsresults/*/aseg.stats`; do s=$(echo ${x} | awk -F"/" '{print $2}'); v=$(grep " WM-hypointensities" ${x} | awk '{print $4}'); echo "${s},${v}"; done | sed '1iPSubject,WM-hypointensities' > fswmh_results.csv [osotolongo@brick03 facehbi]$ sed 's/;/,/g' facehbi_mri.csv | sed '1iSubject,PSubject' > joiner.csv [osotolongo@brick03 facehbi]$ wmh_metrics.pl facehbi [osotolongo@brick03 facehbi]$ join -t"," joiner.csv facehbi_wmh_metrics.csv | awk -F"," '{print $2","$3}' > pgs_wmh.csv [osotolongo@brick03 facehbi]$ (head -n 1 pgs_wmh.csv; tail -n +2 pgs_wmh.csv | sort -t",") > pgs_wmh_sorted.csv [osotolongo@brick03 facehbi]$ join -t"," fswmh_results.csv pgs_wmh_sorted.csv > compare.csv [osotolongo@brick03 facehbi]$ head compare.csv
y despues se pegan,
[osotolongo@brick03 f5cehbi]$ cat compare.csv /old_nas/facehbi/compare.csv | grep -v PSubject | sed '1iPSubject,WM-hypointensities,WMH' > compare_v05.csv
y el script en R hace el plot,
- myplot.r
w <- read.csv("compare_v05.csv") wl <- lm(w$WMH ~ w$WM.hypointensities) pwl <- summary(wl) postscript("wmh_comparing.ps", width=1024, height=600, bg="white") plot(w$WM.hypointensities, w$WMH, main="WMH comparison pgs VS FS", xlab="Freesurfer WMH", ylab="pgs WMH", pch=15, col="blue") abline(wl$coefficients[1],wl$coefficients[2],col="red",lwd=3) text(2500,25000, bquote(R^2 == .(round(pwl$adj.r.squared, digits=3))), cex=2) dev.off()
Poniendo las tres visitas de FACEHBI ($n \sim 550$) juntas tenemos algo asi,
Guardando WMH en XNAT (deprecated)
Para guardar los valores calculados de WMH en XNAT necesito un archivo con el formato,
PSubject,WMH B014,1240.071167 B021,15298.660156 B030,906.359924 B035,966.427979 B037,3140.890869 B038,2161.114502 ... ... ...
Ahora voy a aprovechar dos funciones de interaccion con XNAT; xput_res y xget_res. (ver XNATACE lib)
- xput_res toma la referencia a un hash y lo sube como un .json. Asi que para cada sujeto descrito en el archivo de datos habria que armar un pequeño hash e invocar la funcion,
foreach my $sbj (sort keys %wmhs){ my $experiment = xget_mri($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $sbj); my %wmh_data; $wmh_data{'WMH'} = $wmhs{$sbj}; xcreate_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'data'); xput_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'data', 'wmh.json', \%wmh_data); }
- xget_res toma a localizacion de un archivo .json y devuelve un hash. Asi que por cada sujeto invoco la funcion y escribo el contenido del hash,
foreach my $sbj (sort keys %subjects){ my $experiment = xget_mri($xconf{'HOST'}, $xconf{'JSESSION'}, $xprj, $sbj); my $label = xget_sbj_data($xconf{'HOST'}, $xconf{'JSESSION'}, $sbj, 'label'); my %wmh_data = xget_res($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'data', 'wmh.json'); print "$label,$wmh_data{'WMH'}\n"; }
Los scripts xnat_up_wmh.pl y xnat_get_wmh.pl son entonces sencillos de usar para subir y bajar el contenido del archivo CSV.
Calculando WMH en XNAT
Hemos conseguido configurar el pipeline pgsWMH para que calcule el volumen de WMH directamente en XNAT. (Lo archivos xml aqui) ¿Porque si ya disponiamos de un metodo valido?
- Al integrar el calculo de WMH en un pipeline de XNAT nos ahorramos lo pasos mas costosos (en esfuerzo humano y tiempo) que son localizar, extraer y organizar los DICOM, así como convertir los T1w y T2w a formato nifti.
- Si el pipeline esta corectamente configurado los datos se guardan silenciosamente en XNAT y estan disponibles automaticamente para consultarse en cualquier momento.
- La localizacion de los sujeto para los que aun no se ha calculado el volumen de WMH es ahora trivial y su calculo por separado aun mas simple.
- Al dismminuir el nivel de interaccion necesario para este calculo, los posibles errores de procesamiento intermedio se disminuyen a un minimo.
Así que, ahora para calcular el WMH de un proyecto completo basta hacer algo como,
[osotolongo@brick03 f5cehbi]$ for x in `xnatapic list_experiments --project_id f5cehbi --modality MRI`; do xnatapic run_pipeline --project_id f5cehbi --pipeline pgsWMH --experiment_id ${x}; done
Dependiendo lo bien escritos que esten los TAGs de las imagenes, puede haber algun error, pero llegara un email anunciandolo.
Ahora, con un pequeño cambio en el script de bajar los datos de WMH, podemos obtener los volumenes en mm3.
my %wmh_data = xget_res_data($xconf{'HOST'}, $xconf{'JSESSION'}, $experiment, 'WMH', 'wmh.json'); if (exists($wmh_data{'WMH mm3'}) and $wmh_data{'WMH mm3'}){ print "$label,$wmh_data{'WMH mm3'}\n"; }else{ print "$label,NA\n"; }
Esto es gracias a que la funcion xget_res_data es capaz de bajar los datos de cualquier JSON y almacenarlos en un hash. Solo necesitamos saber la estructura del JSON y preguntar por la variable que se quiera.
Entonces, el resultado se baja como,
[osotolongo@brick03 f5cehbi]$ xnat_get_wmh.pl -x f5cehbi Subject_ID,WMH F204,5162.239258 F246,663.757324 F230,3814.697266 F227,379.562378 F013,5701.064941 F132,387.191772 F033,1088.142334 F073,4213.333008 F027,7450.103516 F074,794.410706 F071,3265.380859 ... ... ...
o si lo deseamos, especificando un archivo de output.
Ahora tenemso el problema de ir ejecutando por lotes, pero utilizando la API de XNAT esto deberia ser simple. Algo como,
[osotolongo@brick03 mri_face]$ for x in `xnatapic list_experiments --project_id mriface --modality MRI`; do if [[ -z `curl -f -X GET -b "JSESSIONID=542119A78B8CD068BACB5862961601F3" "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
Debria darnos la lista de los experimentos que NO tienen el archivo wmh.json. Asi que bastaria con enviar el pipeline a ejecutarse sobre estos experimentos.
Por ejemplo,
[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 XNAT05_E00040 XNAT05_E00066 XNAT_E00937 XNAT_E00938 XNAT_E00939 XNAT_E00940 [osotolongo@brick03 mri_face]$ for x in $(seq 37 40); do xnatapic run_pipeline --project_id mriface --pipeline pgsWMH --experiment_id XNAT_E009${x}; done [osotolongo@brick03 mri_face]$ queue JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON) 323751 fast pgsWMH.XNAT_E00937 xnat R 2:55 1 brick01 323752 fast pgsWMH.XNAT_E00938 xnat R 2:55 1 brick02 323753 fast pgsWMH.XNAT_E00939 xnat R 2:51 1 brick05 323754 fast pgsWMH.XNAT_E00940 xnat R 2:51 1 brick04
Y si, las dos primeras las he excluido porque ya se que han fallado antes por no tener bien el T2w. Entonces lo que voy a hacer es excluir estos experimentos. Voy a hacer un archivo wmh.json para los fallidos,
[osotolongo@brick03 mri_face]$ cat wmh.json {"ResultSet":{"Result":[{"ID":"NA","WMH voxels":"NA","WMH mm3":"NA"}]}}
y entonces,
[osotolongo@brick03 mri_face]$ xnatapic get_jsession FE98FA6B6AB5520A3A329757C3BC1127 [osotolongo@brick03 mri_face]$ curl -f -X PUT -b "JSESSIONID=FE98FA6B6AB5520A3A329757C3BC1127" "http://detritus.fundacioace.com:8088/data/experiments/XNAT05_E00040/resources/WMH/files/wmh.json?overwrite=true" -F file="@wmh.json" [osotolongo@brick03 mri_face]$ curl -f -X PUT -b "JSESSIONID=FE98FA6B6AB5520A3A329757C3BC1127" "http://detritus.fundacioace.com:8088/data/experiments/XNAT05_E00066/resources/WMH/files/wmh.json?overwrite=true" -F file="@wmh.json"
y cuando haga la proxima busqueda estos resultados no apareceran.