Analisis de PET-tau
Procesando
Basicamente se han de hacer algunas tareas que conducen a dos metricas distintas.
Registro de PET-Tau a espacio nativo MRI
Aplicar mascaras provenientes de la segmentacion FS a imagen PET-Tau
Calcular imagen SUVR segun valor en ROR
Calcular metricas sin corregir
Hacer PVC a imagen registrada
Calcular SUVR con correccion
Registro a T1w
El primer paso es colocar la imagen tau (4D de 6 slices) en espacio nativo T1. Debemos haber procesado el T1 con Freesurfer, asi que traemos primero el struc de FS. Ahora hacemos un split de la imagen tau y registramos cada imagen resultante al T1. La imagen tau final sera el valor medio de cada una de estas imagenes registradas.
Utilizamos las 6 imágenes de 5min para la correccion de movimiento
- tau_reg.sh
#!/bin/sh
study=$1
shift
id=$1
shift
wdir=$1
shift
src=$1
shift
td=${wdir}'/.tmp_'${id}
if [ ! -d "$td" ]; then
mkdir $td
fi
debug=0
${FREESURFER_HOME}/bin/mri_convert --in_type mgz --out_type nii ${SUBJECTS_DIR}/${study}_${id}/mri/rawavg.mgz ${wdir}/${id}_struc.nii.gz
${FSLDIR}/bin/fslsplit ${src} ${td}/${id}_piece_ -t
for x in ${td}/${id}_piece_*; do
${ANTS_PATH}/antsRegistrationSyNQuick.sh -d 3 -f ${wdir}/${id}_struc.nii.gz -m ${x} -t a -o ${x%.nii.gz}_movingToFixed_;
${ANTS_PATH}/antsApplyTransforms -d 3 -r ${wdir}/${id}_struc.nii.gz -i ${x} -t ${x%.nii.gz}_movingToFixed_0GenericAffine.mat -o ${x%.nii.gz}_reg.nii.gz;
done
a=`for i in ${td}/*_reg.nii.gz; do echo " $i"; done`
${FSLDIR}/bin/fslmerge -t ${td}/${id}_corr.nii.gz $a
${FSLDIR}/bin/fslmaths ${td}/${id}_corr.nii.gz -Tmean ${wdir}/${id}_tau.nii.gz
if [ $debug = 0 ] ; then
rm -rf ${td}/*_piece_*
fi
Y para ir preparando el wrapper
el call, usando el modulo SLURM.pm se hace como,
#Registro de PET a T1w
$ptask{'command'} = $ENV{'PIPEDIR'}."/bin/tau_reg.sh ".$study." ".$subject." ".$w_dir." ".$spet{'tau'};
$ptask{'filename'} = $outdir.'/'.$subject.'_tau_reg.sh';
$ptask{'output'} = $outdir.'/tau_reg_'.$subject;
my $job_id = send2slurm(\%ptask);
push @r_jobs, $job_id;
Tras efectuar el registro se puede generar el informe para ver como ha ido. Aqui el procedimiento es muy parecido al de FBB asi que no vale la pena comentar nada
- make_tau_report.pl
#!/usr/bin/perl
# Copyright 2021 O. Sotolongo <asqwerty@gmail.com>
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);
my $full = 0;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
last if /^--$/;
if (/^-f/) {$full = 1}
if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/make_pet_report.hlp'; exit;}
}
my $study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/make_pet_report.hlp'; exit;}
my %std = load_project($study);
my $w_dir=$std{'WORKING'};
my $d_dir=$std{'DATA'};
# Redirect ouput to logfile (do it only when everything is fine)
my $debug = "$d_dir/.debug_report.log";
open STDOUT, ">$debug" or die "Can't redirect stdout";
open STDERR, ">&STDOUT" or die "Can't dup stdout";
my $order = "slicesdir -o ";
my @names = find(file => 'name' => "*_tau.nii.gz", in => $w_dir);
foreach my $name (sort @names){
(my $struc = $name) =~ s/_tau/_struc/;
$order .= $name.' '.$struc." ";
}
chdir $w_dir;
print "$order\n";
system($order);
dirmove('slicesdir', 'taus');
Hacer mascaras
Ahora debemos seleccionar las regiones de interes donde queremos medir el SUVR. Hay varias opciones, como las regiones de Braak donde normalmente se observa la progresion del AD, un combinacion de estas regiones o una meta-region que definamos.
Las definiciones de cada region se han tomado de los metodos de ADNI para el AV1451.
[osotolongo@brick03 atau]$ ls /nas/software/neuro.dev/lib/tau/*.roi
/nas/software/neuro.dev/lib/tau/braak_12.roi /nas/software/neuro.dev/lib/tau/braak_4.roi
/nas/software/neuro.dev/lib/tau/braak_1.roi /nas/software/neuro.dev/lib/tau/braak_56.roi
/nas/software/neuro.dev/lib/tau/braak_2.roi /nas/software/neuro.dev/lib/tau/braak_5.roi
/nas/software/neuro.dev/lib/tau/braak_34.roi /nas/software/neuro.dev/lib/tau/braak_6.roi
/nas/software/neuro.dev/lib/tau/braak_3.roi /nas/software/neuro.dev/lib/tau/meta_temporal.roi
[osotolongo@brick03 atau]$ cat /nas/software/neuro.dev/lib/tau/braak_3.roi
1016,L_parahippocampal
1007,L_fusiform
1013,L_lingual
18,L_amygdala
2016,R_parahippocampal
2007,R_fusiform
2013,R_lingual
54,R_amygdala
Ahora tomamos la segmentacion de FS (aparc+aseg), la llevamos a espacio nativo y seleccionamos, segun los LUT, los elementos de cada ROI. Se hace una mascara con cada ROI y simplemente se suman para tener la mascara de la ROI deseada.
Es bastante rapido de hacer
- get_troi.sh
#!/bin/sh
subject=$1
shift
tmp_dir=$1
shift
roi=$1
shift
if [ ! -f ${tmp_dir}/rois/register.dat ]; then
mkdir -p ${tmp_dir}/rois/;
tkregister2 --mov $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --noedit --s ${subject} --regheader --reg ${tmp_dir}/rois/register.dat;
fi
if [ ! -f ${tmp_dir}/all_aseg.nii.gz ]; then
mri_label2vol --seg $SUBJECTS_DIR/${subject}/mri/aparc+aseg.mgz --temp $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --o ${tmp_dir}/all_aseg.nii.gz --reg ${tmp_dir}/rois/register.dat;
fi
mkdir ${tmp_dir}/rois/${roi};
for x in `cat ${PIPEDIR}/lib/tau/${roi}.roi`; do
sleep 5;
rlabel=$(echo ${x} | awk -F"," '{print $1}');
nlabel=$(echo ${x} | awk -F"," '{print $2}');
${FSLDIR}/bin/fslmaths ${tmp_dir}/all_aseg.nii.gz -uthr ${rlabel} -thr ${rlabel} -div ${rlabel} ${tmp_dir}/rois/${roi}/${nlabel};
done
a=$(for x in ${tmp_dir}/rois/${roi}/*.nii.gz; do echo "${x} -add "; done)
a=$(echo ${a} | sed 's/\(.*\)-add$/\1/')
#echo ${a}
${FSLDIR}/bin/fslmaths ${a} ${tmp_dir}/rois/${roi}
Esto en el wrapper ha de hacerse para cada roi que deseemos
#Hacer mascara para cada ROI
$ptask{'job_name'} = 'tau_rois_'.$subject.'_'.$study;
$ptask{'dependency'} = 'afterok:'.$job_id;
foreach my $roi (@rois){
$ptask{'output'} = $outdir.'/tau_roi_'.$roi.'_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_roi_'.$roi.'.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_troi.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject.' '.$roi;
$mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/'.$roi.'.nii.gz ';
send2slurm(\%ptask);
}
Normalizacion
Mascara para normalizar
Vamos a normalizar por la sustancia gris del cerebelo. Para garantizar que no haya contaminacion en las mediciones de la materia gris del cerebelo tomamos solo la zona inferior (inferior cerebellum gray matter). Para ello se toma el atlas del cerebelo de SUIT en espacio MNI que proporciona AFNI. Se calcula el registro del T1w a espacio MNI y con la inversa se lleva este atlas a espacio nativo. De aqui se seleccionan las regiones inferiores del cerebelo y se hace una mascara del cerebelo calculado por FS. Ahora bien, dado que el cerebelo de SUIT no encaja exactamente con el segmentado de FS, tomamos la parte inferior y la parte superior (de SUIT) y a cada una se le hace un smooth con un kernel gaussiano de 8mm y se comparan. Si la inferior es mayor que la superior y pertenece a la segmentacion de FS, entonces se toma este punto como valido. La mascara resultante es la que se utilizara para normalizar. Esto parece muy complicado pero son dos operaciones matematicas, $B(I-S)$, donde
\[
B(x) = \{
\begin{array}{cc}
1, & x>0 \\
0, & \textrm{otherwise}
\end{array}
\].
Un poco mas complicado pero similar a la extraccion de ROIs
- get_tref_cgm.sh
#!/bin/sh
subject=$1
shift
tmp_dir=$1
shift
roi=cerebgm.ref
#Primero sacar el cerebelo de FS
if [ ! -f ${tmp_dir}/rois/register.dat ]; then
mkdir -p ${tmp_dir}/rois/;
tkregister2 --mov $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --noedit --s ${subject} --regheader --reg ${tmp_dir}/rois/register.dat;
#tkregister2 --mov $SUBJECTS_DIR/${subject}/mri/nu.mgz --noedit --s ${subject} --regheader --reg ${tmp_dir}/rois/register.dat;
fi
if [ ! -f ${tmp_dir}/all_aseg.nii.gz ]; then
mri_label2vol --seg $SUBJECTS_DIR/${subject}/mri/aparc+aseg.mgz --temp $SUBJECTS_DIR/${subject}/mri/rawavg.mgz --o ${tmp_dir}/all_aseg.nii.gz --reg ${tmp_dir}/rois/register.dat;
#mri_label2vol --seg $SUBJECTS_DIR/${subject}/mri/aparc+aseg.mgz --temp $SUBJECTS_DIR/${subject}/mri/nu.mgz
fi
mkdir ${tmp_dir}/rois/${roi%.ref};
for x in `cat ${PIPEDIR}/lib/tau/${roi}`; do
sleep 5;
rlabel=$(echo ${x} | awk -F"," '{print $1}');
nlabel=$(echo ${x} | awk -F"," '{print $2}');
${FSLDIR}/bin/fslmaths ${tmp_dir}/all_aseg.nii.gz -uthr ${rlabel} -thr ${rlabel} -div ${rlabel} ${tmp_dir}/rois/${roi%.ref}/${nlabel}
done
a=$(for x in ${tmp_dir}/rois/${roi%.ref}/*.nii.gz; do echo "${x} -add "; done)
a=$(echo ${a} | sed 's/\(.*\)-add$/\1/')
#echo ${a}
${FSLDIR}/bin/fslmaths ${a} ${tmp_dir}/rois/${roi%.ref}.nii.gz
#Ahora sacar el mapa de cerebelo de SUIT
${FREESURFER_HOME}/bin/mri_convert --in_type mgz --out_type nii ${SUBJECTS_DIR}/${subject}/mri/rawavg.mgz ${tmp_dir}/${subject}_tmp.nii.gz
${FSLDIR}/bin/fslreorient2std ${tmp_dir}/${subject}_tmp ${tmp_dir}/${subject}_struc
${ANTSPATH}/antsRegistrationSyNQuick.sh -d 3 -f ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -m ${tmp_dir}/${subject}_struc.nii.gz -o ${tmp_dir}/structToMNI_ -t s
${ANTSPATH}/antsApplyTransforms -d 3 -i ${PIPEDIR}/lib/tau/Cerebellum-MNIsegment.nii.gz -r ${tmp_dir}/${subject}_struc.nii.gz -t [${tmp_dir}/structToMNI_0GenericAffine.mat, 1] -t ${tmp_dir}/structToMNI_1InverseWarp.nii.gz -n GenericLabel -o ${tmp_dir}/CerebinNS.nii.gz
#Dejar solo a parte de abajo?
mkdir ${tmp_dir}/rois/incsuit
for x in `cat ${PIPEDIR}/lib/tau/incsuit.ref`; do
${FSLDIR}/bin/fslmaths ${tmp_dir}/CerebinNS.nii.gz -uthr ${x} -thr ${x} -div ${x} ${tmp_dir}/rois/incsuit/${x};
done
a=$(for x in ${tmp_dir}/rois/incsuit/*.nii.gz; do echo "${x} -add "; done)
a=$(echo ${a} | sed 's/\(.*\)-add$/\1/')
${FSLDIR}/bin/fslmaths ${a} ${tmp_dir}/rois/incsuit.nii.gz
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/incsuit.nii.gz -kernel gauss 3.4 -fmean ${tmp_dir}/rois/incsuit_smooth.nii.gz
#Excluir la parte de arriba correctamente
mkdir ${tmp_dir}/rois/excsuit
for x in `cat ${PIPEDIR}/lib/tau/excsuit.ref`; do
${FSLDIR}/bin/fslmaths ${tmp_dir}/CerebinNS.nii.gz -uthr ${x} -thr ${x} -div ${x} ${tmp_dir}/rois/excsuit/${x};
done
a=$(for x in ${tmp_dir}/rois/excsuit/*.nii.gz; do echo "${x} -add "; done)
a=$(echo ${a} | sed 's/\(.*\)-add$/\1/')
${FSLDIR}/bin/fslmaths ${a} ${tmp_dir}/rois/excsuit.nii.gz
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/excsuit.nii.gz -kernel gauss 3.4 -fmean ${tmp_dir}/rois/excsuit_smooth.nii.gz
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/incsuit_smooth.nii.gz -sub ${tmp_dir}/rois/excsuit_smooth.nii.gz -bin ${tmp_dir}/rois/icereb_mask.nii.gz
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/${roi%.ref} -mas ${tmp_dir}/rois/icereb_mask.nii.gz ${tmp_dir}/rois/icgm.nii.gz
Y en el wrapper se lanza como un proceso independiente en SLURM
#Hacer mascara de cerebelo
$ptask{'output'} = $outdir.'/tau_cgm_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_cgm.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_tref_cgm.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject;
$mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/icgm.nii.gz';
send2slurm(\%ptask);
Merge Masks
Calculo de valores e imagen SUVR
PVC correction
Usando la imagen registrada y las mascaras 4D que hemos construido, calcularemos la correccion de volumen parcial (PVC) correspondiente a cada mascara. Vamos a usar PETPVC, desarrollada en University College London para corregir usando el metodo GTM.
Esto es una sola orden simple
- pvc.sh
#!/bin/sh
subject=$1
shift
wdir=$1
shift
petpvc -i ${wdir}/${subject}_tau.nii.gz -m ${wdir}/${subject}_masks.nii.gz -o ${wdir}/${subject}_pvc.csv -p GTM -x 6.0 -y 6.0 -z 6.0
y los lanzamos desde el wrapper cuando tengamos el SUVR
#PVC
$ptask{'output'} = $outdir.'/tau_pvc_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_pvc.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/pvc.sh '.$subject.' '.$w_dir;
$ptask{'dependency'} = 'afterok:'.$sjob_id;
my $pjob_id = send2slurm(\%ptask);
push @p_jobs, $pjob_id;
PVC correction (metodos)
Por defecto hemos elegido el GMT como metodo de correccion pero es posible elegir otro metodo de manera relativamente sencilla.
Esta, por ejemplo, es la implementación del Multi-target correction (MTC),
- petmtc.pl
#!/usr/bin/perl
# Copyright 2021 O. Sotolongo <asqwerty@gmail.com>
#
# This script run the partial volumen correction method of MTC
# over a PET image using a supplied set of masks. The masks file
# is a 4D image where the last 3D should be the reference region.
# The outputs are a txt file with mean values for each mask and a SUVR
# PET image
use strict; use warnings;
use File::Temp qw(tempdir);
use Data::Dump qw(dump);
my $ifile;
my $ofile;
my $mask;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
last if /^--$/;
if (/^-i/) {$ifile = shift; chomp($ifile);}
if (/^-m/) {$mask = shift; chomp($mask);}
if (/^-o/) {$ofile = shift; chomp($ofile);}
}
my $tdir = tempdir( CLEANUP => 1);
(my $mtc = $ifile) =~ s/_tau.nii.gz/_mtc.nii.gz/;
my $gomtc = 'petpvc -i '.$ifile.' -m '.$mask.' -o '.$mtc.'.nii.gz -p MTC -x 6.0 -y 6.0 -z 6.0';
system($gomtc);
my $splord = "$ENV{'FSLDIR'}/bin/fslsplit $mask $tdir/masks -t";
system($splord);
opendir my $dir, "$tdir" or die "Cannot find temp dir\n";
my @files = grep !/^\./, readdir $dir;
close $dir;
my $mean;
open OF, ">$ofile" or die "Could not open output file\n";
print OF "REGION\tMEAN\n";
foreach my $imask (@files){
(my $tag)= $imask =~ /masks(\d+)\.nii\.gz/;
$tag += 1;
my $ord = "fslstats ".$mtc. " -k ".$tdir.'/'.$imask." -M";
$mean = qx/$ord/;
chomp $mean;
print OF "$tag\t$mean\n";
}
close OF;
(my $suvr = $ifile) =~ s/_tau.nii.gz/_mtc_suvr.nii.gz/;
my $ord = "fslmaths $ifile -div $mean $suvr";
system($ord);
que puede añadirse al wrapper con solo algunas lineas
#PVC MTC
$ptask{'output'} = $outdir.'/tau_mtc_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_mtc.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/petmtc.pl -i '.$w_dir.'/'.$subject.'_tau.nii.gz -m '.$w_dir.'/'.$subject.'_masks.nii.gz -o '.$w_dir.'/'.$subject.'_mtc.csv';
$ptask{'dependency'} = 'afterok:'.$sjob_id;
$pjob_id = send2slurm(\%ptask);
push @p_jobs, $pjob_id;
Pero que tambien habria que añadir como opcion en el archivo de las metricas (Ver mas abajo)
Metricas
Cuando tenemos el output del PVC (y los SUVR sin corregir), hemos de escribir toda la informacion en una tabla que agrupe las metricas y los sujetos. Ha de ser una tabla para los SUVR corregidos y la misma tablas para los SUVR sin corregir.
Pero como tenemos los resultados en el mismo formato, solo hay que escribir un algoritmo
- tau_metrics.pl
#!/usr/bin/perl
# Copyright 2021 O. Sotolongo <asqwerty@gmail.com>
use strict; use warnings;
use FSMetrics qw(tau_rois);
use NEURO4 qw(print_help load_project cut_shit);
use Data::Dump qw(dump);
my $cfile="";
my @rois = tau_rois();
my $style = "";
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
last if /^--$/;
if (/^-cut/) { $cfile = shift; chomp($cfile);}
if (/^-r/) {$style = shift; chomp($style);}
}
my $study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/pet_metrics.hlp'; exit;}
my %std = load_project($study);
my $w_dir=$std{'WORKING'};
my $data_dir=$std{'DATA'};
my $db = $data_dir.'/'.$study.'_pet.csv';
our @subjects = cut_shit($db, $data_dir."/".$cfile);
my @subs = ("pvc", "unc");
if($style){@rois = tau_rois($style);}
my $norm = @rois;
my %measures;
foreach my $subject (@subjects){
foreach my $msub (@subs){
my $ifile = $w_dir.'/'.$subject.'_'.$msub.'.csv';
open IDF, "<$ifile" or next;
while (<IDF>){
if(/\d\t\d+\.*\d*/){
my ($index, $tau) = /(\d)\t(\d+\.*\d*)/;
$measures{$subject}{$msub}[$index-1] = $tau;
}
}
close IDF;
}
}
#dump %measures;
foreach my $msub (@subs){
my $ofile = $data_dir.'/'.$study."_tau_suvr_".$msub.".csv";
print "Writing $ofile\n";
open ODF, ">$ofile";
print ODF "Subject";
foreach my $roi (@rois){
print ODF ", $roi";
}
print ODF "\n";
foreach my $subject (@subjects){
if(exists($measures{$subject}) && exists($measures{$subject}{$msub}) && exists($measures{$subject}{$msub}[$norm]) && $measures{$subject}{$msub}[$norm]){
print ODF "$subject";
for (my $i=0; $i<$norm; $i++){
my $mean = 0;
if(exists($measures{$subject}{$msub}[$i]) && $measures{$subject}{$msub}[$i]){
$mean = $measures{$subject}{$msub}[$i]/$measures{$subject}{$msub}[$norm];
}
print ODF ", $mean";
}
print ODF "\n";
}
}
close ODF;
}
Esto ha de lanzarse tras hacer los PVC
# Calculating Metrics
$warn{'command'} = $ENV{'PIPEDIR'}."/bin/tau_metrics.pl ".$study;
$warn{'filename'} = $outdir.'/tau_metrics.sh';
$warn{'job_name'} = 'tau_metrics_'.$study;
$warn{'time'} = '2:0:0'; #si no ha terminado en X horas matalo
$warn{'mailtype'} = 'FAIL,END'; #email cuando termine o falle
$warn{'output'} = $outdir.'/tau_metrics';
$warn{'dependency'} = 'afterok:'.join(',',@p_jobs);
send2slurm(\%warn);
Mas sobre metodos de PVC
Para extraer las metricas siguiendo un metodo extra de correccion habra que cambiar la linea,
my @subs = ("pvc", "unc");
y añadir la cabecera del metodo. Por ejemplo, para el MTC,
my @subs = ("pvc", "unc","mtc");
Wrapper
Al final nos ha quedado un algoritmo extremadamente paralelizado y mucho mas rapido (aunque mas complejo) que, por ejemplo el procesamiento PET-FBB. Se hacen muchisimos procesos cortos. Dependiendo de cuantas ROI se procesen puede llegar a dividirse un sujeto en mas de 10 tareas relativamente rapidas. Por ejemplo, el procesamiento por defecto de 57 sujetos se divide en 675 tareas independientes con su arbol de dependencias bien definido.
Todo el procesamiento Tau, del NIfTI a las metricas finales
- tau_proc.pl
#!/usr/bin/perl
# Copyright 2021 O. Sotolongo <asqwerty@gmail.com>
use strict; use warnings;
use Data::Dump qw(dump);
use File::Find::Rule;
use File::Copy::Recursive qw(dirmove);
use NEURO4 qw(check_pet check_subj load_project print_help check_or_make cut_shit);
use SLURM qw(send2slurm);
use FSMetrics qw(tau_rois);
my $cfile="";
my $time = '2:0:0';
my $style = "";
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
last if /^--$/;
if (/^-cut/) { $cfile = shift; chomp($cfile);}
if (/^-time/) {$time = shift; chomp($time);}
if (/^-r/) {$style = shift; chomp($style);}
if (/^-h/) { print_help $ENV{'PIPEDIR'}.'/doc/tau_reg.hlp'; exit;}
}
my $study = shift;
unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/tau_reg.hlp'; exit;}
my %std = load_project($study);
my $subj_dir = $ENV{'SUBJECTS_DIR'};
my $w_dir=$std{'WORKING'};
my $db = $std{'DATA'}.'/'.$study.'_pet.csv';
my $data_dir=$std{'DATA'};
my $outdir = "$std{'DATA'}/slurm";
check_or_make($outdir);
print "Collecting needed files\n";
my @pets = cut_shit($db, $data_dir.'/'.$cfile);
#defino aqui las propiedades comunes de ejecucion
my %ptask;
$ptask{'job_name'} = 'tau_reg_'.$study;
$ptask{'cpus'} = 4;
$ptask{'time'} = $time;
my @ok_pets;
my @rois = tau_rois($style);
my @r_jobs;
my @p_jobs;
print "Running shit\n";
foreach my $subject (@pets){
my %spet = check_pet($std{'DATA'},$subject);
my %smri = check_subj($std{'DATA'},$subject);
if($spet{'tau'} && $smri{'T1w'}){
push @ok_pets, $subject;
if(exists($ptask{'mailtype'})){ delete($ptask{'mailtype'}) };
if(exists($ptask{'dependency'})){ delete($ptask{'dependency'}) };
#Registro de PET a T1w
$ptask{'command'} = $ENV{'PIPEDIR'}."/bin/tau_reg.sh ".$study." ".$subject." ".$w_dir." ".$spet{'tau'};
$ptask{'filename'} = $outdir.'/'.$subject.'_tau_reg.sh';
$ptask{'output'} = $outdir.'/tau_reg_'.$subject;
my $job_id = send2slurm(\%ptask);
push @r_jobs, $job_id;
my $mask_chain = "";
my $reg_img = $w_dir.'/'.$subject.'_tau.nii.gz';
#Hacer mascara para cada ROI
$ptask{'job_name'} = 'tau_rois_'.$subject.'_'.$study;
$ptask{'dependency'} = 'afterok:'.$job_id;
foreach my $roi (@rois){
$ptask{'output'} = $outdir.'/tau_roi_'.$roi.'_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_roi_'.$roi.'.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_troi.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject.' '.$roi;
$mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/'.$roi.'.nii.gz ';
send2slurm(\%ptask);
}
#Hacer mascara de cerebelo
$ptask{'output'} = $outdir.'/tau_cgm_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_cgm.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_tref_cgm.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject;
$mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/icgm.nii.gz';
send2slurm(\%ptask);
#Juntar todas las mascaras en un 4D
$ptask{'output'} = $outdir.'/tau_merge_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_merge.sh';
$ptask{'command'} = $ENV{'FSLDIR'}.'/bin/fslmerge -t '.$w_dir.'/'.$subject.'_masks '.$mask_chain;
#$ptask{'mailtype'} = 'FAIL,END';
$ptask{'dependency'} = 'singleton';
my $mjob_id = send2slurm(\%ptask);
#Calculo SUVR
$ptask{'output'} = $outdir.'/tau_suvr_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_suvr.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/petunc.pl -i '.$w_dir.'/'.$subject.'_tau.nii.gz -m '.$w_dir.'/'.$subject.'_masks.nii.gz -o '.$w_dir.'/'.$subject.'_unc.csv';
$ptask{'dependency'} = 'afterok:'.$mjob_id;
my $sjob_id = send2slurm(\%ptask);
#PVC
$ptask{'output'} = $outdir.'/tau_pvc_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_pvc.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/pvc.sh '.$subject.' '.$w_dir;
$ptask{'dependency'} = 'afterok:'.$sjob_id;
my $pjob_id = send2slurm(\%ptask);
push @p_jobs, $pjob_id;
}
}
# Generating report
my %warn;
$warn{'command'} = $ENV{'PIPEDIR'}."/bin/make_tau_report.pl ".$study;
$warn{'filename'} = $outdir.'/tau_report.sh';
$warn{'job_name'} = 'tau_reg_'.$study;
$warn{'time'} = '2:0:0'; #si no ha terminado en X horas matalo
#$warn{'mailtype'} = 'FAIL,END'; #email cuando termine o falle
$warn{'output'} = $outdir.'/tau_report';
$warn{'dependency'} = 'afterok:'.join(',',@r_jobs);
send2slurm(\%warn);
# Calculating Metrics
$warn{'command'} = $ENV{'PIPEDIR'}."/bin/tau_metrics.pl ".($style?" -r $style ":"").$study;
$warn{'filename'} = $outdir.'/tau_metrics.sh';
$warn{'job_name'} = 'tau_metrics_'.$study;
$warn{'time'} = '2:0:0'; #si no ha terminado en X horas matalo
$warn{'mailtype'} = 'FAIL,END'; #email cuando termine o falle
$warn{'output'} = $outdir.'/tau_metrics';
$warn{'dependency'} = 'afterok:'.join(',',@p_jobs);
send2slurm(\%warn);
Testing with ADNI
Procesamiento
Ahora tengo que probar que todo funciona correctamente y comprar los resultados contra los reportados en ADNI. el primer paso debe ser escoger los sujetos que voy a procesar.
Lo que hago es ir a R, instalar el paquete ADNIMERGE y encontrar la correspondencia entre los RID y los PTID ya que los resultados del AV1451 se reportan con el RID pero las imagenes se identifican por el PTID.
Cuando tenga las lista de los PTID debo ir a la pagina de ADNI y buscar, en es ta lista, aquellos que tienen MRI T1 y PET-Tau al mismo tiempo. Entonces simplemente escogo algunos sujetos y bajo las dos imagenes (MRI y PET) que matcheen bien (uy que palabreja). Los tags son muy distintos aqui asi que debo revisar bien el proceso de conversion a BIDS. Pero con unos 50 MRI que se convierten automaticamente tengo para empezar. Despues tocara revisar que los PETs esten OK.
Entonces sigo el procedimiento usual de cualquier proyecto, hago la DB, convierto a BIDS, paso el recon a todas las MRI. Primero vamos a construir el proyecto de la manera usual y a ejecutar la segmentacion de Freesurfer y sacamos al final los datos que necesitemos para comparar.
[osotolongo@brick03 ttau]$ make_proj.pl ttau /old_nas/ttau/raw/ADNI/
[osotolongo@brick03 ttau]$ update_mri_db.pl ttau
[osotolongo@brick03 ttau]$ ls
bids ids.csv raw ttau_mri.csv working
[osotolongo@brick03 ttau]$ cd bids/
[osotolongo@brick03 bids]$ dcm2bids_scaffold
[osotolongo@brick03 bids]$ dcm2bids_helper -d /old_nas/ttau/raw/ADNI/002_S_0413/
Example in:
/old_nas/ttau/bids/tmp_dcm2bids/helper
[osotolongo@brick03 bids]$ cat /old_nas/ttau/bids/tmp_dcm2bids/helper/002_002_S_0413_Accelerated_Sagittal_MPRAGE_20170621132338.json | jq ".SeriesDescription"
"Accelerated_Sagittal_MPRAGE"
[osotolongo@brick03 bids]$ cat /old_nas/ttau/bids/tmp_dcm2bids/helper/324840_002_S_0413_ADNI_Brain_PET__Raw_20170621152221.json | jq ".SeriesDescription"
"ADNI_Brain_PET:_Raw"
...
...
...
[osotolongo@brick03 ttau]$ bulk2bids.pl ttau
...
...
...
[osotolongo@brick03 ttau]$ precon.pl ttau
Bueno, no es tan simple porque hay la conversion a BIDS hay que hacerla a mano en algunos casos, revisar tres veces todo, etc, pero basicamente es lo que hay que hacer para cualquier proyecto.
Y entonces es que lanzo el tau_proc.pl.
[osotolongo@brick03 ttau]$ tau_proc.pl ttau
Collecting needed files
Running shit
[osotolongo@brick03 ttau]$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
24897 fast tau_rois osotolon PD 0:00 1 (Dependency)
24898 fast tau_rois osotolon PD 0:00 1 (Dependency)
24899 fast tau_rois osotolon PD 0:00 1 (Dependency)
...
25468 fast tau_rois osotolon R 0:01 1 brick04
25479 fast tau_rois osotolon R 0:01 1 brick04
25490 fast tau_rois osotolon R 0:01 1 brick04
25501 fast tau_rois osotolon R 0:01 1 brick04
[osotolongo@brick03 ttau]$ squeue | wc -l
619
Puede que falle alguno de los sujetos en cuyo caso se repite para este y las metricas salen para todos
[osotolongo@brick03 ttau]$ ls ttau_tau_suvr_*
ttau_tau_suvr_ewm_pvc.csv ttau_tau_suvr_ewm_unc.csv ttau_tau_suvr_icgm_pvc.csv ttau_tau_suvr_icgm_unc.csv
Voy a coger uno cualquiera,
[osotolongo@brick03 ttau]$ head ttau_tau_suvr_ewm_pvc.csv
Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6
0001, 1.44483375159353, 1.542037154334, 1.68067637885554, 1.64863738043504, 1.7229901461639, 1.67999486862888
0002, 1.48708315130333, 1.33939614581816, 1.61751371292656, 1.61595553613902, 1.73767778263191, 1.71424688121936
0003, 1.9650963436256, 1.8743393102064, 1.79994898471214, 1.64887682735088, 1.66973455323989, 1.7693272672532
0004, 1.86972155598144, 1.47172525122697, 1.77916508872797, 1.7092461469606, 1.88224178228654, 1.83245160031674
0005, 1.37786202498854, 1.27174541179284, 1.52766873243927, 1.43915270111392, 1.50966462746348, 1.61955741984337
0006, 1.4514188475063, 1.17718144895938, 1.53763984073279, 1.55371280053709, 1.55088108801574, 1.52661876242198
0007, 2.56283422459893, 1.39576528400058, 2.57994652406417, 2.66047839283133, 2.43327793033675, 1.79630726983668
0008, 2.677682894831, 1.52617894711733, 2.26320873983255, 2.95362330940057, 2.61091763471126, 1.93361042864299
0009, 1.87077072644842, 1.75846335993224, 1.82703173579494, 1.77024997834227, 1.79149380588898, 1.69506877532751
[osotolongo@brick03 ttau]$ sed 's/;/,/;1iSubject,PTID' ttau_mri.csv > codes.csv
[osotolongo@brick03 ttau]$ join -t, codes.csv ttau_tau_suvr_ewm_pvc.csv > ewm_pvc.csv
[osotolongo@brick03 ttau]$ head ewm_pvc.csv
Subject,PTID, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6
0001,002_S_0413, 1.44483375159353, 1.542037154334, 1.68067637885554, 1.64863738043504, 1.7229901461639, 1.67999486862888
0002,002_S_1155, 1.48708315130333, 1.33939614581816, 1.61751371292656, 1.61595553613902, 1.73767778263191, 1.71424688121936
0003,002_S_1261, 1.9650963436256, 1.8743393102064, 1.79994898471214, 1.64887682735088, 1.66973455323989, 1.7693272672532
0004,002_S_1280, 1.86972155598144, 1.47172525122697, 1.77916508872797, 1.7092461469606, 1.88224178228654, 1.83245160031674
0005,002_S_4213, 1.37786202498854, 1.27174541179284, 1.52766873243927, 1.43915270111392, 1.50966462746348, 1.61955741984337
0006,002_S_4229, 1.4514188475063, 1.17718144895938, 1.53763984073279, 1.55371280053709, 1.55088108801574, 1.52661876242198
0007,002_S_4262, 2.56283422459893, 1.39576528400058, 2.57994652406417, 2.66047839283133, 2.43327793033675, 1.79630726983668
0008,002_S_4521, 2.677682894831, 1.52617894711733, 2.26320873983255, 2.95362330940057, 2.61091763471126, 1.93361042864299
0009,002_S_4654, 1.87077072644842, 1.75846335993224, 1.82703173579494, 1.77024997834227, 1.79149380588898, 1.69506877532751
y voy a buscar las fechas de los PETs, pues las necesitare mas tarde. Primero voy a copiar cualquier DICOM PET para mirar el AcquisitionDate, en un sitio aparte para que no se me confundan con los MRI.
[osotolongo@brick03 raw]$ for x in `ls ADNI/`; do mkdir CHECK/${x}; done
[osotolongo@brick03 raw]$ for x in `ls ADNI/`; do y=$(ls -d ADNI/${x}/* | grep PET); if [[ ${y} ]]; then f=$(find ${y} -type f | head -n 1); cp $f CHECK/${x}/; fi; done
[osotolongo@brick03 raw]$ tree CHECK/
CHECK/
├── 002_S_0413
│ └── ADNI_002_S_0413_PT_ADNI_Brain_PET__Raw_br_raw_20170622122126095_198_S574062_I863053.dcm
├── 002_S_1155
│ └── ADNI_002_S_1155_PT_ADNI_Brain_PET__Raw_br_raw_20170425132632053_105_S558344_I843576.dcm
├── 002_S_1261
│ └── ADNI_002_S_1261_PT_ADNI_Brain_PET__Raw_br_raw_20170320134605163_460_S547104_I831093.dcm
├── 002_S_1280
│ └── ADNI_002_S_1280_PT_ADNI_Brain_PET__Raw_br_raw_20170315120003807_259_S545652_I829443.dcm
├── 002_S_4213
│ └── ADNI_002_S_4213_PT_ADNI_Brain_PET__Raw_br_raw_20170818105234393_5_S598553_I892242.dcm
├── 002_S_4229
│ └── ADNI_002_S_4229_PT_ADNI_Brain_PET__Raw_br_raw_20171006103626160_184_S618186_I915198.dcm
├── 002_S_4262
│ └── ADNI_002_S_4262_PT_ADNI_Brain_PET__Raw_br_raw_20160226102445166_420_S383619_I641202.dcm
├── 002_S_4521
│ └── ADNI_002_S_4521_PT_ADNI_Brain_PET__Raw_br_raw_20160406143337045_210_S411667_I672710.dcm
├── 002_S_4654
│ └── ADNI_002_S_4654_PT_ADNI_Brain_PET__Raw_br_raw_20160422123234789_222_S426421_I689622.dcm
├── 002_S_4799
│ └── ADNI_002_S_4799_PT_ADNI_Brain_PET__Raw_br_raw_20180615124556282_63_S696406_I1010905.dcm
├── 002_S_5178
│ └── ADNI_002_S_5178_PT_ADNI_Brain_PET__Raw_br_raw_20170601120820611_179_S569644_I857832.dcm
├── 002_S_5230
│ └── ADNI_002_S_5230_PT_ADNI_Brain_PET__Raw_br_raw_20170814143412364_155_S594270_I886723.dcm
├── 002_S_6007
│ └── ADNI_002_S_6007_PT_ADNI_Brain_PET__Raw_br_raw_20170427121132154_62_S558985_I844301.dcm
├── 002_S_6009
│ └── ADNI_002_S_6009_PT_ADNI_Brain_PET__Raw_br_raw_20170516133754120_124_S564181_I851391.dcm
├── 002_S_6030
│ └── ADNI_002_S_6030_PT_ADNI_Brain_PET__Raw_br_raw_20170726144651878_339_S588018_I879426.dcm
├── 002_S_6053
│ └── ADNI_002_S_6053_PT_ADNI_Brain_PET__Raw_br_raw_20170828120721764_236_S601531_I895803.dcm
├── 002_S_6103
│ └── ADNI_002_S_6103_PT_ADNI_Brain_PET__Raw_br_raw_20180322120342372_261_S669417_I976654.dcm
├── 002_S_6456
│ └── ADNI_002_S_6456_PT_ADNI_Brain_PET__Raw_br_raw_20180807132825528_274_S713482_I1032082.dcm
├── 002_S_6652
│ └── ADNI_002_S_6652_PT_ADNI_Brain_PET__Raw_br_raw_20190425134639516_506_S818823_I1158635.dcm
├── 002_S_6695
│ └── ADNI_002_S_6695_PT_ADNI_Brain_PET__Raw_br_raw_20190430114756478_164_S819870_I1159991.dcm
├── 003_S_1122
├── 003_S_4119
├── 003_S_4288
├── 003_S_4441
├── 003_S_4644
├── 003_S_5154
├── 003_S_6014
├── 003_S_6067
├── 003_S_6256
├── 003_S_6257
├── 003_S_6258
├── 003_S_6259
├── 003_S_6260
├── 003_S_6264
├── 003_S_6268
├── 003_S_6307
├── 003_S_6432
├── 003_S_6606
├── 003_S_6644
├── 011_S_0021
├── 011_S_6303
├── 011_S_6367
├── 011_S_6465
├── 011_S_6714
├── 012_S_4094
│ └── ADNI_012_S_4094_PT_PET_AC_3D_BRAIN__br_raw_20191217141112685_123_S906083_I1267876.dcm
├── 012_S_4188
│ └── ADNI_012_S_4188_PT_PET_AC_3D_BRAIN__br_raw_20191120135313594_198_S898709_I1258492.dcm
├── 012_S_4643
│ └── ADNI_012_S_4643_PT_PET_AC_3D_BRAIN__br_raw_20160609062740359_108_S459983_I729585.dcm
├── 012_S_6073
│ └── ADNI_012_S_6073_PT_PET_AC_3D_BRAIN__br_raw_20180111122501862_185_S649546_I952534.dcm
├── 012_S_6760
│ └── ADNI_012_S_6760_PT_PET_AC_3D_BRAIN__br_raw_20190812082418032_76_S851921_I1200704.dcm
├── 013_S_2389
│ └── ADNI_013_S_2389_PT_PET_Brain_DYNAMIC_6X5--__br_raw_20180308153859437_613_S665196_I971459.dcm
├── 013_S_6206
│ └── ADNI_013_S_6206_PT_PET_Brain_DYNAMIC_6X5__br_raw_20180726141550509_565_S709086_I1026668.dcm
├── 013_S_6768
│ └── ADNI_013_S_6768_PT_PET_Brain_AV-1451_DYNAMIC__br_raw_20191016130947649_483_S884501_I1241169.dcm
├── 013_S_6780
│ └── ADNI_013_S_6780_PT_PET_Brain_AV-1451_DYNAMIC__br_raw_20191010165106316_85_S883192_I1239588.dcm
├── 014_S_4401
│ └── ADNI_014_S_4401_PT_ADNI_Brain_PET__Raw_br_raw_20171121162212250_113_S635424_I935447.dcm
├── 023_S_0031
│ └── ADNI_023_S_0031_PT_PET_Brain_TAU__2___AC__br_raw_20190423173518567_12_S817876_I1157437.dcm
├── 067_S_0056
└── 067_S_0059
Me faltan,
[osotolongo@brick03 raw]$ for x in `ls ADNI/`; do y=$(ls -d ADNI/${x}/* | grep 1451); if [[ ${y} ]]; then f=$(find ${y} -type f | head -n 1); cp $f CHECK/${x}/; fi; done
[osotolongo@brick03 raw]$ tree CHECK/
CHECK/
├── 002_S_0413
│ └── ADNI_002_S_0413_PT_ADNI_Brain_PET__Raw_br_raw_20170622122126095_198_S574062_I863053.dcm
├── 002_S_1155
│ └── ADNI_002_S_1155_PT_ADNI_Brain_PET__Raw_br_raw_20170425132632053_105_S558344_I843576.dcm
├── 002_S_1261
│ └── ADNI_002_S_1261_PT_ADNI_Brain_PET__Raw_br_raw_20170320134605163_460_S547104_I831093.dcm
├── 002_S_1280
│ └── ADNI_002_S_1280_PT_ADNI_Brain_PET__Raw_br_raw_20170315120003807_259_S545652_I829443.dcm
├── 002_S_4213
│ └── ADNI_002_S_4213_PT_ADNI_Brain_PET__Raw_br_raw_20170818105234393_5_S598553_I892242.dcm
├── 002_S_4229
│ └── ADNI_002_S_4229_PT_ADNI_Brain_PET__Raw_br_raw_20171006103626160_184_S618186_I915198.dcm
├── 002_S_4262
│ └── ADNI_002_S_4262_PT_ADNI_Brain_PET__Raw_br_raw_20160226102445166_420_S383619_I641202.dcm
├── 002_S_4521
│ └── ADNI_002_S_4521_PT_ADNI_Brain_PET__Raw_br_raw_20160406143337045_210_S411667_I672710.dcm
├── 002_S_4654
│ └── ADNI_002_S_4654_PT_ADNI_Brain_PET__Raw_br_raw_20160422123234789_222_S426421_I689622.dcm
├── 002_S_4799
│ └── ADNI_002_S_4799_PT_ADNI_Brain_PET__Raw_br_raw_20180615124556282_63_S696406_I1010905.dcm
├── 002_S_5178
│ └── ADNI_002_S_5178_PT_ADNI_Brain_PET__Raw_br_raw_20170601120820611_179_S569644_I857832.dcm
├── 002_S_5230
│ └── ADNI_002_S_5230_PT_ADNI_Brain_PET__Raw_br_raw_20170814143412364_155_S594270_I886723.dcm
├── 002_S_6007
│ └── ADNI_002_S_6007_PT_ADNI_Brain_PET__Raw_br_raw_20170427121132154_62_S558985_I844301.dcm
├── 002_S_6009
│ └── ADNI_002_S_6009_PT_ADNI_Brain_PET__Raw_br_raw_20170516133754120_124_S564181_I851391.dcm
├── 002_S_6030
│ └── ADNI_002_S_6030_PT_ADNI_Brain_PET__Raw_br_raw_20170726144651878_339_S588018_I879426.dcm
├── 002_S_6053
│ └── ADNI_002_S_6053_PT_ADNI_Brain_PET__Raw_br_raw_20170828120721764_236_S601531_I895803.dcm
├── 002_S_6103
│ └── ADNI_002_S_6103_PT_ADNI_Brain_PET__Raw_br_raw_20180322120342372_261_S669417_I976654.dcm
├── 002_S_6456
│ └── ADNI_002_S_6456_PT_ADNI_Brain_PET__Raw_br_raw_20180807132825528_274_S713482_I1032082.dcm
├── 002_S_6652
│ └── ADNI_002_S_6652_PT_ADNI_Brain_PET__Raw_br_raw_20190425134639516_506_S818823_I1158635.dcm
├── 002_S_6695
│ └── ADNI_002_S_6695_PT_ADNI_Brain_PET__Raw_br_raw_20190430114756478_164_S819870_I1159991.dcm
├── 003_S_1122
│ └── ADNI_003_S_1122_PT_ADNI3_av-1451__AC__br_raw_20171019121957192_634_S622410_I920110.dcm
├── 003_S_4119
│ └── ADNI_003_S_4119_PT_ADNI3_av-1451__AC__br_raw_20191024101644788_318_S887685_I1245124.dcm
├── 003_S_4288
│ └── ADNI_003_S_4288_PT_ADNI3_av-1451__AC__br_raw_20180301152746113_94_S663344_I969266.dcm
├── 003_S_4441
│ └── ADNI_003_S_4441_PT_ADNI3_av-1451__AC__br_raw_20190425141845749_602_S818825_I1158637.dcm
├── 003_S_4644
│ └── ADNI_003_S_4644_PT_ADNI3_av-1451__AC__br_raw_20180418155808610_128_S677845_I987120.dcm
├── 003_S_5154
│ └── ADNI_003_S_5154_PT_ADNI3_av-1451__AC__br_raw_20180913144127683_168_S726965_I1048272.dcm
├── 003_S_6014
│ └── ADNI_003_S_6014_PT_ADNI3_av-1451__AC__br_raw_20180411104621424_608_S675426_I984154.dcm
├── 003_S_6067
│ └── ADNI_003_S_6067_PT_ADNI3_av-1451__AC__br_raw_20180314150648891_332_S667020_I973677.dcm
├── 003_S_6256
│ └── ADNI_003_S_6256_PT_ADNI3_AV-1451__AC__br_raw_20180710101539021_58_S702929_I1019022.dcm
├── 003_S_6257
│ └── ADNI_003_S_6257_PT_ADNI3_av-1451__AC__br_raw_20180330133302322_350_S672236_I980135.dcm
├── 003_S_6258
│ └── ADNI_003_S_6258_PT_ADNI3_av-1451__AC__br_raw_20180502165119864_598_S681920_I992116.dcm
├── 003_S_6259
│ └── ADNI_003_S_6259_PT_ADNI3_av-1451__AC__br_raw_20180518161343520_30_S687345_I999771.dcm
├── 003_S_6260
│ └── ADNI_003_S_6260_PT_ADNI3_AV-1451__AC__br_raw_20180730160933950_194_S710691_I1028655.dcm
├── 003_S_6264
│ └── ADNI_003_S_6264_PT_ADNI3_av-1451__AC__br_raw_20181206135747230_262_S755211_I1082828.dcm
├── 003_S_6268
│ └── ADNI_003_S_6268_PT_ADNI3_av-1451__AC__br_raw_20180822110628592_620_S720092_I1039770.dcm
├── 003_S_6307
│ └── ADNI_003_S_6307_PT_ADNI3_av-1451__AC__br_raw_20190814170919421_474_S858406_I1208729.dcm
├── 003_S_6432
│ └── ADNI_003_S_6432_PT_ADNI3_av-1451__AC__br_raw_20190726154159037_552_S845026_I1191956.dcm
├── 003_S_6606
│ └── ADNI_003_S_6606_PT_ADNI3_av-1451__AC__br_raw_20190304155918547_66_S803257_I1138703.dcm
├── 003_S_6644
│ └── ADNI_003_S_6644_PT_ADNI3_av-1451__AC__br_raw_20190219125639038_134_S798424_I1132421.dcm
├── 011_S_0021
│ └── ADNI_011_S_0021_PT_ADNI3_AV1451__AC___br_raw_20180209092946315_455_S657348_I961935.dcm
├── 011_S_6303
│ └── ADNI_011_S_6303_PT_ADNI3_AV1451__AC___br_raw_20180508091001800_37_S683222_I994488.dcm
├── 011_S_6367
│ └── ADNI_011_S_6367_PT_ADNI3_AV1451__AC___br_raw_20180605090017512_268_S691948_I1005530.dcm
├── 011_S_6465
│ └── ADNI_011_S_6465_PT_ADNI3_AV1451__AC___br_raw_20180816124040325_474_S717388_I1036653.dcm
├── 011_S_6714
│ └── ADNI_011_S_6714_PT_ADNI3_AV1451__AC___br_raw_20190725105029277_209_S844033_I1190652.dcm
├── 012_S_4094
│ └── ADNI_012_S_4094_PT_PET_AC_3D_BRAIN__br_raw_20191217141112685_123_S906083_I1267876.dcm
├── 012_S_4188
│ └── ADNI_012_S_4188_PT_PET_AC_3D_BRAIN__br_raw_20191120135313594_198_S898709_I1258492.dcm
├── 012_S_4643
│ └── ADNI_012_S_4643_PT_PET_AC_3D_BRAIN__br_raw_20160609062740359_108_S459983_I729585.dcm
├── 012_S_6073
│ └── ADNI_012_S_6073_PT_PET_AC_3D_BRAIN__br_raw_20180111122501862_185_S649546_I952534.dcm
├── 012_S_6760
│ └── ADNI_012_S_6760_PT_PET_AC_3D_BRAIN__br_raw_20190812082418032_76_S851921_I1200704.dcm
├── 013_S_2389
│ └── ADNI_013_S_2389_PT_PET_Brain_DYNAMIC_6X5--__br_raw_20180308153859437_613_S665196_I971459.dcm
├── 013_S_6206
│ └── ADNI_013_S_6206_PT_PET_Brain_DYNAMIC_6X5__br_raw_20180726141550509_565_S709086_I1026668.dcm
├── 013_S_6768
│ └── ADNI_013_S_6768_PT_PET_Brain_AV-1451_DYNAMIC__br_raw_20191016130947649_483_S884501_I1241169.dcm
├── 013_S_6780
│ └── ADNI_013_S_6780_PT_PET_Brain_AV-1451_DYNAMIC__br_raw_20191010165106316_85_S883192_I1239588.dcm
├── 014_S_4401
│ └── ADNI_014_S_4401_PT_ADNI_Brain_PET__Raw_br_raw_20171121162212250_113_S635424_I935447.dcm
├── 023_S_0031
│ └── ADNI_023_S_0031_PT_PET_Brain_TAU__2___AC__br_raw_20190423173518567_12_S817876_I1157437.dcm
├── 067_S_0056
│ └── ADNI_067_S_0056_PT_AV1451.TAU_Dyn_6x5min_2Di_336_2z_AllPass__AC___br_raw_20190115195420753_396_S785330_I1116637.dcm
└── 067_S_0059
└── ADNI_067_S_0059_PT_AV1451.TAU_Dyn_6x5min_2Di_336_2z_AllPass__AC___br_raw_20171214155739398_493_S644301_I946067.dcm
57 directories, 57 files
yahoo! Ahora si que puedo sacar las fechas.
[osotolongo@brick03 ttau]$ for x in /nas/data/ttau/raw/CHECK/*; do f=$(find ${x} -type f | head -n 1); d=$(dckey -k "AcquisitionDate" ${f} 2>&1 ); n=$(echo ${x} | awk -F"/" '{print $7}'); echo "${n},${d}"; done | sed '1iPTID,Date'| sort -t, -n -k 1 > dates.csv
[osotolongo@brick03 ttau]$ head dates.csv
PTID,Date
002_S_0413,20170621
002_S_1155,20170424
002_S_1261,20170315
002_S_1280,20170313
002_S_4213,20170817
002_S_4229,20171003
002_S_4262,20160225
002_S_4521,20160405
002_S_4654,20160421
[osotolongo@brick03 ttau]$ sort -t, -n -k 2 ewm_pvc.csv > ewm_pvc_sorted.csv
[osotolongo@brick03 ttau]$ join -t, -1 1 -2 2 dates.csv ewm_pvc_sorted.csv > ewm_pvc_ok.csv
[osotolongo@brick03 ttau]$ head ewm_pvc_ok.csv
PTID,Date,Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6
002_S_0413,20170621,0001, 1.44483375159353, 1.542037154334, 1.68067637885554, 1.64863738043504, 1.7229901461639, 1.67999486862888
002_S_1155,20170424,0002, 1.48708315130333, 1.33939614581816, 1.61751371292656, 1.61595553613902, 1.73767778263191, 1.71424688121936
002_S_1261,20170315,0003, 1.9650963436256, 1.8743393102064, 1.79994898471214, 1.64887682735088, 1.66973455323989, 1.7693272672532
002_S_1280,20170313,0004, 1.86972155598144, 1.47172525122697, 1.77916508872797, 1.7092461469606, 1.88224178228654, 1.83245160031674
002_S_4213,20170817,0005, 1.37786202498854, 1.27174541179284, 1.52766873243927, 1.43915270111392, 1.50966462746348, 1.61955741984337
002_S_4229,20171003,0006, 1.4514188475063, 1.17718144895938, 1.53763984073279, 1.55371280053709, 1.55088108801574, 1.52661876242198
002_S_4262,20160225,0007, 2.56283422459893, 1.39576528400058, 2.57994652406417, 2.66047839283133, 2.43327793033675, 1.79630726983668
002_S_4521,20160405,0008, 2.677682894831, 1.52617894711733, 2.26320873983255, 2.95362330940057, 2.61091763471126, 1.93361042864299
002_S_4654,20160421,0009, 1.87077072644842, 1.75846335993224, 1.82703173579494, 1.77024997834227, 1.79149380588898, 1.69506877532751
Data ADNI
A ver que hay dentro,
[osotolongo@brick03 tau]$ head -n 1 UCBERKELEYAV1451_PVC_01_15_21.csv | sed 's/,/\n/g' | cat -n
1 "RID"
2 "VISCODE"
3 "VISCODE2"
4 "EXAMDATE"
5 "INFERIOR_CEREBGM_SUVR"
6 "INFERIOR_CEREBGM_VOLUME"
7 "HEMIWM_SUVR"
8 "HEMIWM_VOLUME"
9 "BRAAK1_SUVR"
10 "BRAAK1_VOLUME"
11 "BRAAK34_SUVR"
12 "BRAAK34_VOLUME"
13 "META_TEMPORAL_SUVR"
14 "META_TEMPORAL_VOLUME"
15 "BRAAK56_SUVR"
16 "BRAAK56_VOLUME"
17 "BRAAK2_SUVR"
18 "BRAAK2_VOLUME"
19 "AMYGDALA_SUVR"
20 "AMYGDALA_VOLUME"
21 "BRAIN_STEM_SUVR"
22 "BRAIN_STEM_VOLUME"
23 "LEFT_MIDDLEFR_SUVR"
24 "LEFT_MIDDLEFR_VOLUME"
25 "LEFT_ORBITOFR_SUVR"
26 "LEFT_ORBITOFR_VOLUME"
27 "LEFT_PARSFR_SUVR"
28 "LEFT_PARSFR_VOLUME"
29 "LEFT_ACCUMBENS_AREA_SUVR"
30 "LEFT_ACCUMBENS_AREA_VOLUME"
31 "LEFT_AMYGDALA_SUVR"
32 "LEFT_AMYGDALA_VOLUME"
33 "LEFT_CAUDATE_SUVR"
34 "LEFT_CAUDATE_VOLUME"
35 "LEFT_HIPPOCAMPUS_SUVR"
36 "LEFT_HIPPOCAMPUS_VOLUME"
37 "LEFT_PALLIDUM_SUVR"
38 "LEFT_PALLIDUM_VOLUME"
39 "LEFT_PUTAMEN_SUVR"
40 "LEFT_PUTAMEN_VOLUME"
41 "LEFT_THALAMUS_PROPER_SUVR"
42 "LEFT_THALAMUS_PROPER_VOLUME"
43 "RIGHT_MIDDLEFR_SUVR"
44 "RIGHT_MIDDLEFR_VOLUME"
45 "RIGHT_ORBITOFR_SUVR"
46 "RIGHT_ORBITOFR_VOLUME"
47 "RIGHT_PARSFR_SUVR"
48 "RIGHT_PARSFR_VOLUME"
49 "RIGHT_ACCUMBENS_AREA_SUVR"
50 "RIGHT_ACCUMBENS_AREA_VOLUME"
51 "RIGHT_AMYGDALA_SUVR"
52 "RIGHT_AMYGDALA_VOLUME"
53 "RIGHT_CAUDATE_SUVR"
54 "RIGHT_CAUDATE_VOLUME"
55 "RIGHT_HIPPOCAMPUS_SUVR"
56 "RIGHT_HIPPOCAMPUS_VOLUME"
57 "RIGHT_PALLIDUM_SUVR"
58 "RIGHT_PALLIDUM_VOLUME"
59 "RIGHT_PUTAMEN_SUVR"
60 "RIGHT_PUTAMEN_VOLUME"
61 "RIGHT_THALAMUS_PROPER_SUVR"
62 "RIGHT_THALAMUS_PROPER_VOLUME"
63 "CHOROID_SUVR"
64 "CHOROID_VOLUME"
65 "CTX_LH_BANKSSTS_SUVR"
66 "CTX_LH_BANKSSTS_VOLUME"
67 "CTX_LH_CAUDALANTERIORCINGULATE_SUVR"
68 "CTX_LH_CAUDALANTERIORCINGULATE_VOLUME"
69 "CTX_LH_CUNEUS_SUVR"
70 "CTX_LH_CUNEUS_VOLUME"
71 "CTX_LH_ENTORHINAL_SUVR"
72 "CTX_LH_ENTORHINAL_VOLUME"
73 "CTX_LH_FUSIFORM_SUVR"
74 "CTX_LH_FUSIFORM_VOLUME"
75 "CTX_LH_INFERIORPARIETAL_SUVR"
76 "CTX_LH_INFERIORPARIETAL_VOLUME"
77 "CTX_LH_INFERIORTEMPORAL_SUVR"
78 "CTX_LH_INFERIORTEMPORAL_VOLUME"
79 "CTX_LH_INSULA_SUVR"
80 "CTX_LH_INSULA_VOLUME"
81 "CTX_LH_ISTHMUSCINGULATE_SUVR"
82 "CTX_LH_ISTHMUSCINGULATE_VOLUME"
83 "CTX_LH_LATERALOCCIPITAL_SUVR"
84 "CTX_LH_LATERALOCCIPITAL_VOLUME"
85 "CTX_LH_LINGUAL_SUVR"
86 "CTX_LH_LINGUAL_VOLUME"
87 "CTX_LH_MIDDLETEMPORAL_SUVR"
88 "CTX_LH_MIDDLETEMPORAL_VOLUME"
89 "CTX_LH_PARACENTRAL_SUVR"
90 "CTX_LH_PARACENTRAL_VOLUME"
91 "CTX_LH_PARAHIPPOCAMPAL_SUVR"
92 "CTX_LH_PARAHIPPOCAMPAL_VOLUME"
93 "CTX_LH_PERICALCARINE_SUVR"
94 "CTX_LH_PERICALCARINE_VOLUME"
95 "CTX_LH_POSTCENTRAL_SUVR"
96 "CTX_LH_POSTCENTRAL_VOLUME"
97 "CTX_LH_POSTERIORCINGULATE_SUVR"
98 "CTX_LH_POSTERIORCINGULATE_VOLUME"
99 "CTX_LH_PRECENTRAL_SUVR"
100 "CTX_LH_PRECENTRAL_VOLUME"
101 "CTX_LH_PRECUNEUS_SUVR"
102 "CTX_LH_PRECUNEUS_VOLUME"
103 "CTX_LH_ROSTRALANTERIORCINGULATE_SUVR"
104 "CTX_LH_ROSTRALANTERIORCINGULATE_VOLUME"
105 "CTX_LH_SUPERIORFRONTAL_SUVR"
106 "CTX_LH_SUPERIORFRONTAL_VOLUME"
107 "CTX_LH_SUPERIORPARIETAL_SUVR"
108 "CTX_LH_SUPERIORPARIETAL_VOLUME"
109 "CTX_LH_SUPERIORTEMPORAL_SUVR"
110 "CTX_LH_SUPERIORTEMPORAL_VOLUME"
111 "CTX_LH_SUPRAMARGINAL_SUVR"
112 "CTX_LH_SUPRAMARGINAL_VOLUME"
113 "CTX_LH_TEMPORALPOLE_SUVR"
114 "CTX_LH_TEMPORALPOLE_VOLUME"
115 "CTX_LH_TRANSVERSETEMPORAL_SUVR"
116 "CTX_LH_TRANSVERSETEMPORAL_VOLUME"
117 "CTX_RH_BANKSSTS_SUVR"
118 "CTX_RH_BANKSSTS_VOLUME"
119 "CTX_RH_CAUDALANTERIORCINGULATE_SUVR"
120 "CTX_RH_CAUDALANTERIORCINGULATE_VOLUME"
121 "CTX_RH_CUNEUS_SUVR"
122 "CTX_RH_CUNEUS_VOLUME"
123 "CTX_RH_ENTORHINAL_SUVR"
124 "CTX_RH_ENTORHINAL_VOLUME"
125 "CTX_RH_FUSIFORM_SUVR"
126 "CTX_RH_FUSIFORM_VOLUME"
127 "CTX_RH_INFERIORPARIETAL_SUVR"
128 "CTX_RH_INFERIORPARIETAL_VOLUME"
129 "CTX_RH_INFERIORTEMPORAL_SUVR"
130 "CTX_RH_INFERIORTEMPORAL_VOLUME"
131 "CTX_RH_INSULA_SUVR"
132 "CTX_RH_INSULA_VOLUME"
133 "CTX_RH_ISTHMUSCINGULATE_SUVR"
134 "CTX_RH_ISTHMUSCINGULATE_VOLUME"
135 "CTX_RH_LATERALOCCIPITAL_SUVR"
136 "CTX_RH_LATERALOCCIPITAL_VOLUME"
137 "CTX_RH_LINGUAL_SUVR"
138 "CTX_RH_LINGUAL_VOLUME"
139 "CTX_RH_MIDDLETEMPORAL_SUVR"
140 "CTX_RH_MIDDLETEMPORAL_VOLUME"
141 "CTX_RH_PARACENTRAL_SUVR"
142 "CTX_RH_PARACENTRAL_VOLUME"
143 "CTX_RH_PARAHIPPOCAMPAL_SUVR"
144 "CTX_RH_PARAHIPPOCAMPAL_VOLUME"
145 "CTX_RH_PERICALCARINE_SUVR"
146 "CTX_RH_PERICALCARINE_VOLUME"
147 "CTX_RH_POSTCENTRAL_SUVR"
148 "CTX_RH_POSTCENTRAL_VOLUME"
149 "CTX_RH_POSTERIORCINGULATE_SUVR"
150 "CTX_RH_POSTERIORCINGULATE_VOLUME"
151 "CTX_RH_PRECENTRAL_SUVR"
152 "CTX_RH_PRECENTRAL_VOLUME"
153 "CTX_RH_PRECUNEUS_SUVR"
154 "CTX_RH_PRECUNEUS_VOLUME"
155 "CTX_RH_ROSTRALANTERIORCINGULATE_SUVR"
156 "CTX_RH_ROSTRALANTERIORCINGULATE_VOLUME"
157 "CTX_RH_SUPERIORFRONTAL_SUVR"
158 "CTX_RH_SUPERIORFRONTAL_VOLUME"
159 "CTX_RH_SUPERIORPARIETAL_SUVR"
160 "CTX_RH_SUPERIORPARIETAL_VOLUME"
161 "CTX_RH_SUPERIORTEMPORAL_SUVR"
162 "CTX_RH_SUPERIORTEMPORAL_VOLUME"
163 "CTX_RH_SUPRAMARGINAL_SUVR"
164 "CTX_RH_SUPRAMARGINAL_VOLUME"
165 "CTX_RH_TEMPORALPOLE_SUVR"
166 "CTX_RH_TEMPORALPOLE_VOLUME"
167 "CTX_RH_TRANSVERSETEMPORAL_SUVR"
168 "CTX_RH_TRANSVERSETEMPORAL_VOLUME"
169 "OTHER_SUVR"
170 "OTHER_VOLUME"
171 "update_stamp"
Ah, voy a sacar una medida de aqui a ver,
[osotolongo@brick03 tau]$ awk -F"," '{print $1","$4","$9","$7}' UCBERKELEYAV1451_PVC_01_15_21.csv | sed 's/"//g;s/-//g;s/EXAMDATE/Date/' > adni_pvc.csv
[osotolongo@brick03 tau]$ head adni_pvc.csv
RID,Date,BRAAK1_SUVR,HEMIWM_SUVR
21,20180202,2.2462,1.0244
31,20180424,3.1467,0.8688
31,20190423,3.0552,0.8856
56,20180220,2.4762,1.1414
56,20190110,2.2607,1.2007
56,20191203,2.3225,1.1212
59,20171212,2.2866,1.0762
69,20180403,2.343,1.1652
69,20190123,2.2939,1.0749
y a preprocesar un poco
saco del ADNIMERGE la correspondencia de los RID,
[osotolongo@brick03 tau]$ head idcorr.csv
"RID","PTID"
21,"011_S_0021"
31,"023_S_0031"
56,"067_S_0056"
59,"067_S_0059"
69,"100_S_0069"
89,"073_S_0089"
96,"022_S_0096"
112,"127_S_0112"
120,"027_S_0120"
y a mezclar un poco,
[osotolongo@brick03 tau]$ join -t, idcorr_ok.csv adni_pvc.csv > adni_pvc_ok.csv
[osotolongo@brick03 tau]$ head adni_pvc_ok.csv
RID,PTID,Date,BRAAK1_SUVR,HEMIWM_SUVR
21,011_S_0021,20180202,2.2462,1.0244
31,023_S_0031,20180424,3.1467,0.8688
31,023_S_0031,20190423,3.0552,0.8856
56,067_S_0056,20180220,2.4762,1.1414
56,067_S_0056,20190110,2.2607,1.2007
56,067_S_0056,20191203,2.3225,1.1212
59,067_S_0059,20171212,2.2866,1.0762
69,100_S_0069,20180403,2.343,1.1652
69,100_S_0069,20190123,2.2939,1.0749
y vamos a intentar unir los registros correctos,
[osotolongo@brick03 ttau]$ awk -F"," '{print $1"_"$2","$0}' ewm_pvc_ok.csv > ewm_pvc_refact.csv
[osotolongo@brick03 tau]$ awk -F"," '{print $2"_"$3","$0}' adni_pvc_ok.csv > adni_pvc_refact.csv
[osotolongo@brick03 tau]$ head adni_pvc_refact.csv
PTID_Date,RID,PTID,Date,BRAAK1_SUVR,HEMIWM_SUVR
011_S_0021_20180202,21,011_S_0021,20180202,2.2462,1.0244
023_S_0031_20180424,31,023_S_0031,20180424,3.1467,0.8688
023_S_0031_20190423,31,023_S_0031,20190423,3.0552,0.8856
067_S_0056_20180220,56,067_S_0056,20180220,2.4762,1.1414
067_S_0056_20190110,56,067_S_0056,20190110,2.2607,1.2007
067_S_0056_20191203,56,067_S_0056,20191203,2.3225,1.1212
067_S_0059_20171212,59,067_S_0059,20171212,2.2866,1.0762
100_S_0069_20180403,69,100_S_0069,20180403,2.343,1.1652
100_S_0069_20190123,69,100_S_0069,20190123,2.2939,1.0749
[osotolongo@brick03 check]$ pwd
/old_nas/ttau/check
[osotolongo@brick03 check]$ cp /old_nas/ttau/ewm_pvc_refact.csv ./
[osotolongo@brick03 check]$ cp /old_nas/tau/adni_pvc_refact.csv ./
[osotolongo@brick03 check]$ sort -t, -n -k 1 ewm_pvc_refact.csv > ewm_pvc_refact_sorted.csv
[osotolongo@brick03 check]$ sort -t, -n -k 1 adni_pvc_refact.csv > adni_pvc_refact_sorted.csv
Comparando con ADNI
Pegando en bash
[osotolongo@brick03 ttau]$ cat adni_sort.sh
ifile=$1
tfile=$(mktemp adni1.XXXXXXX)
tfile2=$(mktemp adni2.XXXXX)
join -t, idcorr_ok.csv ${ifile} > ${tfile}
awk -F"," '{print $2"_"$3","$0}' ${tfile} > ${tfile2}
sort -t, -n -k 1 ${tfile2}
[osotolongo@brick03 ttau]$ cat ace_sort.sh
ifile=$1
tfile=$(mktemp tau1.XXXXXXX)
tfile2=$(mktemp tau2.XXXXXXX)
join -t, codes.csv ${ifile} > ${tfile}
sort -t, -n -k 2 ${tfile} > ${tfile2}
join -t, -1 1 -2 2 dates.csv ${tfile2} > ${tfile}
awk -F"," '{print $1"_"$2","$0}' ${tfile} > ${tfile2}
sort -t, -n -k 1 ${tfile2}
Ejemplo con el PVCed
[osotolongo@brick03 ttau]$ awk -F"," '{print $1","$4","$13","$5}' /old_nas/tau/UCBERKELEYAV1451_PVC_01_15_21.csv | sed 's/"//g;s/-//g;s/EXAMDATE/Date/' > adni_mpvc.csv
[osotolongo@brick03 ttau]$ ./adni_sort.sh adni_mpvc.csv > check/adni_mpvc_refact_sorted.csv
[osotolongo@brick03 ttau]$ ./ace_sort.sh ttau_tau_suvr_pvc.csv > check/icgm_mpvc_refact_sorted.csv
[osotolongo@brick03 ttau]$ join -t, check/icgm_mpvc_refact_sorted.csv check/adni_mpvc_refact_sorted.csv > check/check_mpvc.csv
Esto lo puedo mirar con R
> read.csv("check_mpvc.csv", header=TRUE) -> cp
> cp$adni_suvr = cp$META_TEMPORAL_SUVR/cp$INFERIOR_CEREBGM_SUVR
> mp <- lm(cp$meta_temporal ~ cp$adni_suvr)
> summary(mp)
Call:
lm(formula = cp$meta_temporal ~ cp$adni_suvr)
Residuals:
Min 1Q Median 3Q Max
-0.219275 -0.045731 0.003388 0.056858 0.193383
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.31842 0.02868 11.10 1.49e-15 ***
cp$adni_suvr 0.96495 0.01455 66.31 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07806 on 54 degrees of freedom
Multiple R-squared: 0.9879, Adjusted R-squared: 0.9876
F-statistic: 4397 on 1 and 54 DF, p-value: < 2.2e-16
> plot(cp$adni_suvr,cp$meta_temporal,col="darkgreen",pch=15,cex=1.5, main = "Corrected Meta-Temporal", xlab = "ADNI SUVR", ylab = "ACE SUVR")
> abline(mp$coefficients[1],mp$coefficients[2],col="red",lwd=3)
> text(1.6,4.3,expression(ACE == 0.31842 + 0.96495*ADNI))
> text(1.4,4.0, expression(R^2 == 0.9876))
Después de haber revisado todos los procedimientos, comparado de varias maneras con los datos de ADNI y usado varias maneras de normalización hemos de tomar decisiones.
La normalización se hará por inferior cerebellum gray matter (ICGM) que parece ser la mas estable.
Aqui las comparaciones,
| Uncorrected Braak 1 | PVCed Braak 1 |
 |  |
| Uncorrected Braak 3+4 | PVCed Braak 3+4 |
 |  |
| Uncorrected Braak 5+6 | PVCed Braak 5+6 |
 |  |
| Uncorrected Meta Temporal | PVCed Meta Temporal |
 |  |
$R^2$
| | Uncorrected | Corrected |
| Braak 1 | 0.97 | 0.96 |
| Braak 3+4 | 0.99 | 0.98 |
| Braak 5+6 | 0.99 | 0.95 |
| Meta temporal | 0.99 | 0.99 |
$slope$
| | Uncorrected | Corrected |
| Braak 1 | 1.11 | 0.90 |
| Braak 3+4 | 1.05 | 0.98 |
| Braak 5+6 | 1.02 | 0.92 |
| Meta temporal | 1.07 | 0.98 |
Comparando con otro metodo de correccion (MTC)
Comparando Geometric Matrix Transfer (GMT) com Multi-target Correction (MTC)
| Braak 1 | Meta-temporal |
 |  |
y el ajuste,
| | $R^2$ | $slope$ |
| Braak 1 | 0.99 | 1.02 |
| Meta-temporal | 1.0 | 0.98 |
How to cite
PET-Tau scans have been processing following ADNI method, freely available at ADNI website for registered researchers [1, 2]. First, the 4D 30min PET image was separated into 6 5min images for movement correction and registered into the T1w subject native space. The regions of interest for calculating SUVR, as well as for normalization, where built using the Freesurfer 7.2 segmentation of the T1w. Partial volume correction was calculated using GTM method [3] with PETPVC [4]. Corrected and uncorrected SUVR where calculated using inferior cerebellum gray matter as region of reference.
[1] http://loni.adni.usc.edu. Flortaucipir (AV-1451) processing methods. (https://adni.bitbucket.io/reference/docs/UCBERKELEYAV1451/UCBERKELEY_AV1451_Methods_2021-01-14.pdf)
[2] Maass, A., et al., Comparison of multiple tau-PET measures as biomarkers in aging and Alzheimer's disease. Neuroimage, 2017. 157: p. 448-463.
[3] Rousset, O.G., Y. Ma, and A.C. Evans, Correction for partial volume effects in PET:
principle and validation. J Nucl Med, 1998. 39(5): p. 904-11.
[4] BA Thomas, V Cuplov, A Bousse, A Mendes, K Thielemans, BF Hutton, K Erlandsson. PETPVC: a toolbox for performing partial volume correction techniques in positron emission tomography. Physics in Medicine and Biology 61 (22), 7975.
Plotting with ggseg
Para mostrar los valores de $\tau$-SUVR voy a usar ggseg. Esto es un script python que voy a generar a partir del output del analisis.
Para empezar, necesito saber: proyecto, radiotracer, sujeto y modalidad que voy a graficar. Asi que estos valores los voy a meter como inputs de un script Perl y a partir d ahi decido que hacer,
#!/usr/bin/perl
use strict;
use warnings;
use NEURO4 qw(check_pet check_subj load_project);
use Data::Dump qw(dump);
my $tracer;
my $subject;
my $mod = 'unc';
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
if (/^-tracer/) {$tracer = shift; chomp($tracer);}
if (/^-id/) {$subject = shift; chomp($subject);}
if (/^-mod/) {$mod = shift; chomp($mod);}
}
my $study = shift;
die "Should supply project name\n" unless $study;
die "Should supply subject ID\n" unless $subject;
Ahora cargo las propiedades del proyecto y leo los datos en un hash,
my %std = load_project($study);
my $ifile = $std{'DATA'}.'/'.$study.'_tau_suvr_'.$tracer.'_'.$mod.'.csv';
open IDF, "<$ifile";
my $seeker;
my %taud;
my @levels;
while(<IDF>){
unless ($seeker){
@levels = split ", ", $_;
chomp @levels;
}else{
my @values = split ", ", $_;
chomp @values;
if ($values[0] eq $subject){
for (my $i = 1; $i < scalar @values; $i++){
$taud{$levels[$i]} = $values[$i] unless $levels[$i] eq 'extra';
}
}
}
$seeker++;
}
close IDF;
Ahora, segun estos datos toca escribir el script python. En principio no sabemos que ROIs han sido definidas para calcular el SUVR pero si que estan definidas dentro del pipeline. Asi que lo que hago es aislar cada una y traducirlas a las ROI de FS.
my $odata = "data = {\n";
foreach my $tag (sort keys %taud){
my $libfile = $ENV{'PIPEDIR'}.'/lib/tau/'.$tag.'.roi';
open ILF, "<$libfile";
while (<ILF>){
my ($hand, $roi) = /\d+,(\w)_(\w+)/;
unless ($hand eq 'R'){
$odata.="'".$roi."_left': ".$taud{$tag}.",\n";
$odata.="'".$roi."_right': ".$taud{$tag}.",\n";
}
}
close ILF;
}
$odata.="}\n";
Eso lo vuelco dentro del script python,
my $pfile = $std{'DATA'}.'/'.$study.'_'.$tracer.'_'.$subject.'_'.$mod.'.py';
open OPF, ">$pfile";
print OPF '#!/usr/bin/python3'."\n";
print OPF "import ggseg\n";
print OPF "$odata\n";
print OPF "ggseg.plot_dk(data, cmap='Spectral', figsize=(15,15),";
print OPF "background='k', edgecolor='w', bordercolor='gray',";
print OPF "ylabel='PET-".$tracer." SUVR', title='".$subject."')\n";
close OPF;
El call es algo como,
./generate_tau_pyplot.pl -tracer ro948 -id 0075 -mod unc f5cehbi
que nos deja un script python como,
#!/usr/bin/python3
import ggseg
data = {
'entorhinal_left': 0.880505023804403,
'entorhinal_right': 0.880505023804403,
....
'paracentral_left': 1.31710059291554,
'paracentral_right': 1.31710059291554,
}
ggseg.plot_dk(data, cmap='Spectral', figsize=(15,15),background='k', edgecolor='w', bordercolor='gray',ylabel='PET-ro948 SUVR', title='0075')
y cuando se ejecuta podemos salvar una imagen con los distintos valores de $\tau$-SUVR por las ROIs de FS,
PET-Tau FACEHBI project
La particularidad de este proyecto es que hay dos marcadores PI2620 y RO948. En principio esta diferencia se indica en los tags. Los json deberían indicar el marcador,
"SeriesDescription": "[DY_NAC] Dynamic Brain",
"ProtocolName": "Head/FACEHBI_PI2620_",
Así que añadimos un par de reglas al archivo de conversion,
{
dataType: pet,
modalityLabel: tau,
customLabels: pi2620,
criteria: {
SeriesDescription: *DY_NAC*,
ProtocolName : *PI2620*
}
},
{
dataType: pet,
modalityLabel: tau,
customLabels: ro948,
criteria: {
SeriesDescription: *DY_NAC*,
ProtocolName : *RO948*
}
}
Ahora, la configuracion del proyecto es la siguiente,
[osotolongo@brick03 f5cehbi]$ cat ~/.config/neuro/f5cehbi.cfg
DATA = /old_nas/f5cehbi
SRC = /nas/data/f5cehbi/raw
PET = /nas/data/f5cehbi/raw_tau/PI2620
WORKING = /old_nas/f5cehbi/working
BIDS = /old_nas/f5cehbi/bids
y dependiendo del marcador que vayamos a procesar cambiamos el path de PET al correcto,
/nas/data/f5cehbi/raw_tau/PI2620
/nas/data/f5cehbi/raw_tau/RO948
y ejecutamos,
$ pet2bids.pl f5cehbi
Pero, hay alguna imagen mal identificada, y para encontrarla hacemos,
[osotolongo@brick03 f5cehbi]$ ls /nas/data/f5cehbi/raw_tau/PI2620 > upload_pi2620.list
[osotolongo@brick03 f5cehbi]$ find bids/ -name "*pi2620_*tau.json" | awk -F"-" '{print $2}' | awk -F"/" '{print $1}' > onbids.list
[osotolongo@brick03 f5cehbi]$ grep -f upload_pi2620.list f5cehbi_pet.csv | awk -F";" '{print $1}' > subjects_uploaded.list
[osotolongo@brick03 f5cehbi]$ grep -v "`cat onbids.list`" subjects_uploaded.list
0071
0078
y ahora estos los he de copiar manualmente,
[osotolongo@brick03 f5cehbi]$ grep "DY_NAC" bids/tmp_dcm2bids/sub-0071/*.json
bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.json: "SeriesDescription": "[DY_NAC] Dynamic Brain",
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.json bids/sub-0071/pet/sub-0071_pi2620_tau.json
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.nii.gz bids/sub-0071/pet/sub-0071_pi2620_tau.nii.gz
[osotolongo@brick03 f5cehbi]$ grep "DY_NAC" bids/tmp_dcm2bids/sub-0078/*.json
bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.json: "SeriesDescription": "[DY_NAC] Dynamic Brain",
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.json bids/sub-0078/pet/sub-0078_pi2620_tau.json
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.nii.gz bids/sub-0078/pet/sub-0078_pi2620_tau.nii.gz
Y ya estamos. Ahora a procesar, pero antes hay que hacer limpieza,
[osotolongo@brick03 f5cehbi]$ rm -rf working/taus working/.tmp_* working/*_tau.nii.gz working/*_suvr* working/*_masks* working/*_pvc* working/*_mtc* working/*_unc*
Ahora si,
[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer pi2620 f5cehbi
y al terminar tendremos los resultados en los archivos,
[osotolongo@brick03 f5cehbi]$ ls f5cehbi_tau_suvr_pi2620_*
f5cehbi_tau_suvr_pi2620_mtc.csv f5cehbi_tau_suvr_pi2620_pvc.csv f5cehbi_tau_suvr_pi2620_unc.csv
[osotolongo@brick03 f5cehbi]$ head f5cehbi_tau_suvr_pi2620_unc.csv
Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6, extra
0071, 0.642265450913569, 0.445293067511959, 0.579153982059088, 0.747920179741428, 0.72149021446756, 0.610043587665485, 0.456761589685176
0075, 0.712122676066806, 0.38863874053704, 0.572840374589894, 0.800212087107015, 0.729010277863309, 0.662967042884597, 0.394067798859667
0076, 0.83589891768702, 0.763072917986403, 0.98202394592715, 1.1291513524068, 1.07060320066878, 0.942356703562845, 0.809746693748801
0078, 0.862191018953226, 0.627670886199243, 0.862667901143999, 1.10188426829118, 1.04589664481027, 0.884161149229649, 0.63109502246773
0087, 0.830211995236038, 0.456955305688661, 0.528028370642198, 1.02278147223258, 0.706801020407484, 0.549391337405959, 0.484003816604253
0092, 0.745907131688215, 0.688052328589026, 0.991095394357976, 1.35968730905804, 1.30310065337532, 1.1900620384451, 0.807880424703717
0098, 0.702882998845167, 0.711627059244375, 0.856521933131609, 0.990771318090343, 0.936773228442926, 0.842062696773501, 0.764783413964935
0099, 0.823298330211384, 0.705876423260448, 0.875380807683978, 1.1361288650671, 0.961009113607904, 0.86318920811104, 0.92071927514608
0100, 0.857882757785999, 0.656176690905404, 1.12023931254886, 1.17050947142544, 1.30016270374768, 1.19007631440357, 0.857547667880239
Y ahora voy a hacer el otra marcador. Edito el config
[osotolongo@brick03 f5cehbi]$ cat ~/.config/neuro/f5cehbi.cfg
DATA = /old_nas/f5cehbi
SRC = /nas/data/f5cehbi/raw
PET = /nas/data/f5cehbi/raw_tau/RO948
WORKING = /old_nas/f5cehbi/working
BIDS = /old_nas/f5cehbi/bids
Hago los bids,
[osotolongo@brick03 f5cehbi]$ pet2bids.pl f5cehbi
Limpio el directorio de trabajo,
[osotolongo@brick03 f5cehbi]$ cat clean_tau.sh
#!/bin/bash
rm -rf working/taus working/.tmp_* working/*_tau.nii.gz working/*_suvr* working/*_masks* working/*_pvc* working/*_mtc* working/*_unc*
[osotolongo@brick03 f5cehbi]$ ./clean_tau.sh
Lanzo los procesos,
[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer ro948 f5cehbi
y espero los resultados,
[osotolongo@brick03 f5cehbi]$ head f5cehbi_tau_suvr_ro948_pvc.csv
Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6, extra
0075, 1.04088050314465, 0.64527944539737, 1.15851915380217, 1.52799456832476, 2.21913236134934, 2.39967124070898, 1.30373785020011
0115, 0.908310815175864, 0.492099523944246, 1.10060882563261, 1.43477011430239, 1.83175427770736, 1.77909186294406, 0.796637090468415
0120, 1.46854534539977, 1.07698127739582, 1.8255003454564, 1.66630044512906, 2.06223198812989, 2.256736399778, 1.52880200251447
Integrando el proyecto en XNAT
Subir un PET a XNAT es trivial. El problema aqui es que hay dos PETs y los DICOM no contienen suficiente información 
para que XNAT distinga ambos. Asi que tiende a agruparlos en el mismo PET.
Lo que voy a hacer es fijar el experiment_id y asignar a cada uno el PET que le corresponde. Ejemplo,
$ xnatapic upload_dicom --project_id tau0 --subject_id ${sbj} --experiment_id ${sbj}_PI2620 /nas/data/f5cehbi/raw_tau/PI2620/${pet}
O en contexto,
$ while read -r line; do pet=$(echo ${line} | awk -F"," '{print $3}'); sbj=$(echo ${line} | awk -F"," '{print $1}'); if [ -d /nas/data/f5cehbi/raw_tau/PI2620/${pet} ]; then xnatapic upload_dicom --project_id tau0 --subject_id ${sbj} --experiment_id ${sbj}_PI2620 /nas/data/f5cehbi/raw_tau/PI2620/${pet}; else echo "No ${pet}"; fi; done < fpt_matches.csv
Lo que quedaría como,
Pero como llego a tener el archivo necesario para esto (fpt_matches.csv)?
Lo primero es obtener la lista de los PET del archivo compartido y arreglarle los endline,
$ tr -d '\r' < pet_tau_list.csv > pet_tau_list_proper.csv
Ahora tomamos los codigos del proyecto (f5cehbi),
$ sed 's/;/,/' f5cehbi_mri.csv > f5cehbi_codes.csv
$ sort -t, -k 2 f5cehbi_codes.csv > f5cehbi_codes_sorted.csv
Y todo esto lo juntamos,
$ tail -n +2 pet_tau_list_proper.csv | awk -F"," '{print $1","$3}' | sort -t, > fpt_list.csv
$ join -t, -1 2 -2 1 f5cehbi_codes_sorted.csv fpt_list.csv > fpt_matches.csv
Por supuesto que este archivo es necesario solo una vez, para ambos tracers. Pero debe actualizarse cuando tenemos nuevos sujetos.
Ya hemos calculado los SVUR pero ahora necesitamos ponerlo de manera comprensible en un archivo para entregar. Primero sacamos las fechas, aprovechando el formato de almacenamiento en XNAT,
$ xnatapic list_experiments --project_id tau0 --modality PET --label --date > xnat_tau0_experiments.csv
$ grep RO948 xnat_tau0_experiments.csv | awk -F"," '{print $3","$2}' | sed 's/_.*,/,/' | sed '1iPSubject,Date' > ro948_dates.csv
$ grep PI2620 xnat_tau0_experiments.csv | awk -F"," '{print $3","$2}' | sed 's/_.*,/,/' | sed '1iPSubject,Date' > pi2620_dates.csv
Ahora hay que unirlo con los resultados pero necesitamos unificar los codigos,
$ sed 's/;/,/' f5cehbi_mri.csv | sort -t, -k 1 | sed '1iSubject,PSubject' > codes.csv
$ join -t, codes.csv f5cehbi_tau_suvr_pi2620_unc.csv | sed 's/ //g' > pi2620_unc_wcodes.csv
$ (head -n 1 pi2620_unc_wcodes.csv && tail -n +2 pi2620_unc_wcodes.csv | sort -t, -k 2) > pi2620_unc_wcodes_sorted.csv
$ join -t, -1 1 -2 2 pi2620_dates.csv pi2620_unc_wcodes_sorted.csv > tau_tmp/PI2620_UNC.csv
La ultima parte hay que hacerla para cada uno de los sheets que se quieran exportar,
$ ls -l tau_tmp/
total 28
-rw-r--r-- 1 osotolongo osotolongo 625 Apr 19 10:48 INFO.csv
-rw-rw-r-- 1 osotolongo osotolongo 4758 Apr 19 10:49 PI2620_PVC.csv
-rw-rw-r-- 1 osotolongo osotolongo 4870 Apr 19 10:47 PI2620_UNC.csv
-rw-rw-r-- 1 osotolongo osotolongo 1055 Apr 19 10:51 RO948_PVC.csv
-rw-rw-r-- 1 osotolongo osotolongo 1071 Apr 19 10:52 RO948_UNC.csv
y ahora voy a hacer la conversion a xls,
- dir2xls.pl
#!/usr/bin/perl
#
# To convert a bunch of CSV files into a single XLS
#
# Copyleft 2022 O. Sotolongo <asqwerty@gmail.com>
#
use strict;
use warnings;
use Text::CSV qw( csv );
use Spreadsheet::Write;
my $idir;
my $ofile;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
$_ = shift;
last if /^--$/;
if (/^-i/) { $idir = shift; chomp($idir);}
if (/^-o/) { $ofile = shift; chomp($ofile);}
}
die "Should supply results directory\n" unless $idir;
die "Should output filename\n" unless $ofile;
$ofile =~ s/\.(\w*)?$/.xls/;
opendir (DIR, $idir);
my @ifiles = grep(/\.csv/, readdir(DIR));
close DIR;
my $workbook = Spreadsheet::Write->new(file => $ofile, sheet => 'INFO');
my $infof = $idir.'/INFO.csv';
if( -e $infof ){
my $inf_data = csv (in => $infof);
for my $i (0 .. $#{$inf_data}) {
my $row = $inf_data->[$i];
$workbook->addrow($row);
}
}
foreach my $ifile (@ifiles){
unless( $ifile eq 'INFO.csv'){
my $ipath = $idir.'/'.$ifile;
my $idata = csv (in => $ipath); # as array of array
(my $shname = $ifile) =~ s/\.csv$//;
$workbook->addsheet($shname);
for my $i (0 .. $#{$idata}) {
my $row = $idata->[$i];
$workbook->addrow($row);
}
}
}
$workbook->close();
Voy a dejar una copia de esto en github. 
$ ./dir2xls.pl -i tau_tmp/ -o tau_results_20220419.xls
Thinking around
$ sed 's/ //g' f5cehbi_tau_suvr_ro948_unc.csv | awk '{print $0",ro948"}' | sed 's/extra,ro948/extra,tracer/' > f5cehbi_tau_suvr_ro948_unc_tracer.csv
$ sed 's/ //g' f5cehbi_tau_suvr_pi2620_unc.csv | awk '{print $0",pi2620"}' | sed 's/extra,pi2620/extra,tracer/' > f5cehbi_tau_suvr_pi2620_unc_tracer.csv
$ tail -n +2 f5cehbi_tau_suvr_pi2620_unc_tracer.csv > tmp.csv
$ cat f5cehbi_tau_suvr_ro948_unc_tracer.csv tmp.csv > tau_suvr_unc.csv
Y esto voy a compararlo, a ver como se distinguen, en media, los uptakes de los tracers,
- comp_tracers.r
library("ggplot2")
input_file="tau_suvr_unc.csv"
unc <- read.csv(input_file)
rois = c("braak_1","braak_2","braak_3","braak_4","braak_5","braak_6")
for (roi in rois){
output_fig <- paste0(trimws(roi),'.ps')
postscript(output_fig, width=1024, height=600, bg="white")
myp1 <- ggplot(unc, aes_string( x="tracer", y=trimws(roi))) + geom_boxplot()
print(myp1)
dev.off()
}
$ Rscript comp_tracers.r
$ for x in braak_*.ps; do convert "${x}" -rotate 90 ${x%.ps}_mean_suvr.png; done
| braak 1 | braak 2 |
 |  |
| braak 3 | braak 4 |
 |  |
| braak 5 | braak 6 |
 |  |
O todo junto,
$ echo "Subject,SUVR,Tracer,ROI" > header.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$2","$9",braak_1"}' > braak1.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$3","$9",braak_2"}' > braak2.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$4","$9",braak_3"}' > braak3.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$5","$9",braak_4"}' > braak4.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$6","$9",braak_5"}' > braak5.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$7","$9",braak_6"}' > braak6.data
$ cat header.data braak1.data braak2.data braak3.data braak4.data braak5.data braak6.data > tau_suvr_unc_grouped.csv
y
- comp_tau_all.r
library("ggplot2")
theme_set(theme_minimal())
input_file="tau_suvr_unc_grouped.csv"
unc <- read.csv(input_file)
output_fig <- 'output_tau_all_rois.ps'
postscript(output_fig, width=1024, height=600, bg="white")
pp <- ggplot(unc, aes_string( x="ROI", y="SUVR", fill="Tracer")) + geom_boxplot() + geom_hline(yintercept=1, color = "red", size=1)
print(pp)
dev.off()
O puedo comparar los valores de SUVR de tau con los valores de SUVR de FBB. lo mas facil es con los globales,
$ xnat_pullcl.pl -x f5cehbi -o f5cehbi_fbb_cl.csv
$ sed 's/;/,/g;s/Subject/PSubject/' f5cehbi_fbb_cl.csv > f5cehbi_fbb_cl_lone.csv
$ (head -n 1 f5cehbi_fbb_cl_lone.csv && tail -n +2 f5cehbi_fbb_cl_lone.csv | sort -t,) > f5cehbi_fbb_cl_lone_sorted.csv
$ join -t, -1 2 -2 1 codes_sorted.csv f5cehbi_fbb_cl_lone_sorted.csv > f5cehbi_fbb_ok.csv
$ (head -n 1 f5cehbi_fbb_ok.csv && tail -n +2 f5cehbi_fbb_ok.csv | sort -t, -k 2) > f5cehbi_fbb_resorted.csv
$ join -t, -1 1 -2 2 f5cehbi_tau_suvr_pi2620_unc_tracer.csv f5cehbi_fbb_resorted.csv > comp_fbb_pi2620.csv
$ join -t, -1 1 -2 2 f5cehbi_tau_suvr_ro948_unc_tracer.csv f5cehbi_fbb_resorted.csv > comp_fbb_ro948.csv
$ tail -n +2 comp_fbb_ro948.csv > tmp.csv
$ cat comp_fbb_pi2620.csv tmp.csv > comp_fbb_tau.csv
- comp_tracers.r
library("ggplot2")
theme_set(theme_minimal())
input_file="tau_suvr_unc.csv"
unc <- read.csv(input_file)
rois = c("braak_1","braak_2","braak_3","braak_4","braak_5","braak_6")
for (roi in rois){
output_fig <- paste0(trimws(roi),'.ps')
postscript(output_fig, width=1024, height=600, bg="white")
myp1 <- ggplot(unc, aes_string( x="tracer", y=trimws(roi))) + geom_boxplot(aes(fill=factor(tracer)), show.legend=F)+labs(x="Radiotracer", y="SUVR", title=paste0("SUVR for ",roi))
print(myp1)
dev.off()
}
input_file="comp_fbb_tau.csv"
uac <- read.csv(input_file)
for (roi in rois){
output_fig <- paste0('fbb_tau_',roi,'.ps')
postscript(output_fig, width=1024, height=600, bg="white")
myp1 <- ggplot(uac, aes_string(x=roi, y="SUVR", color="factor(tracer)")) + geom_point() + labs(x="tau SUVR", y="FBB SUVR", title=paste0("Comparing FBB Global SUVR and tau ",roi," SUVR"))
myp1 <- myp1 + guides(color = guide_legend(title = "tracer"))
print(myp1)
dev.off
}
$ Rscript comp_tracers.r
$ for x in fbb_tau_braak_*.ps; do convert ${x} -rotate 90 ${x\.ps}.png; done
| braak 1 | braak 2 |
 |  |
| braak 3 | braak 4 |
 |  |
| braak 5 | braak 6 |
 |  |
Con un poco de esfuerzo podemos usar el mismo metodo para medir el SUVR del FBB en las regiones de Braak.
Solo habria que reutilizar el procedimiento de PET-tau y engañarlo para que use los PET-FBB en su lugar. Como voy a suponer que ya he hecho el analisis de FBB, las imagenes puedo bajarmelas de XNAT, asi que me ahorro el registro a espacio nativo del PET.
Solo tengo que sacar la URI del fbb registrado y guardado y bajarme la imagen.
foreach my $subject (@pets){
my $psubject = $pet_data{$subject};
my $fake_tau = $w_dir.'/'.$subject.'_tau.nii.gz';
#print "$fake_tau\n";
my $xrd = xget_pet($xconf_data{'HOST'}, $jsession, $xprj, $psubject);
if ($xrd){
my $xld = 'curl -f -X GET -b "JSESSIONID='.$jsession.'" "'.$xconf_data{'HOST'}.'/data/experiments/'.$xrd.'/files?format=json" 2>/dev/null';
#print "$xld\n";
my $jres = qx/$xld/;
my $xfres = decode_json $jres;
foreach my $xres (@{$xfres->{'ResultSet'}{'Result'}}){
if ($xres->{'file_content'} eq 'PET_reg'){
my $xuri = $xres->{'URI'};
my $grd = 'curl -f -b "JSESSIONID='.$jsession.'" -X GET "'.$xconf_data{'HOST'}.$xuri.'" -o '.$fake_tau.' 2>/dev/null';
system($grd);
}
}
}
....
....
}
El resto del procedimiento es el mimo excepto que la region de referencia debe ser la materia gris del cerebelo. Por lo que debe ser mas rapido todo.
De momento voy a mirar solo los uncorrected.
pi2620 <- read.csv("f5cehbi_tau_suvr_pi2620_unc.csv")
fbb <- read.csv("f5cehbi_fbb_suvr_fake.csv")
ro948 <- read.csv("f5cehbi_tau_suvr_ro948_unc.csv")
merge(fbb, pi2620, by=c("Subject")) -> c1
for (roi in rois){
output_fig <- paste0('fbb_pi2620_local_',roi,'.ps')
postscript(output_fig, width=1024, height=600, bg="white")
myp1 <- ggplot(c1, aes_string( x=paste0(roi,".y"), y=paste0(roi,".x"))) + geom_point() + labs(x="PI2620 SUVR", y="FBB local SUVR", title=paste0("Comparing FBB ",roi, " SUVR and PI2620 ",roi," SUVR"))
print(myp1)
dev.off()
}
merge(fbb, ro948, by=c("Subject")) -> c1
for (roi in rois){
output_fig <- paste0('fbb_ro948_local_',roi,'.ps')
postscript(output_fig, width=1024, height=600, bg="white")
myp1 <- ggplot(c1, aes_string( x=paste0(roi,".y"), y=paste0(roi,".x"))) + geom_point() + labs(x="RO948 SUVR", y="FBB local SUVR", title=paste0("Comparing FBB ",roi, " SUVR and RO948 ",roi," SUVR"))
print(myp1)
dev.off()
}
$ Rscript comp_tracers.r
$ for x in *_local_*.ps; do convert ${x} -rotate 90 ${x%.ps}.png; done
**RO948**
| braak 1 | braak 2 |
 |  |
| braak 3 | braak 4 |
 |  |
| braak 5 | braak 6 |
 |  |
**PI2620**
| braak 1 | braak 2 |
 |  |
| braak 3 | braak 4 |
 |  |
| braak 5 | braak 6 |
 |  |
Si quiero comparar los tracers entre ellos,
$ join -t, -j 1 f5cehbi_tau_suvr_ro948_unc.csv f5cehbi_tau_suvr_pi2620_unc.csv > tau_ropi.csv
$ sed 's/,/ /g' tau_ropi.csv > tau_ropi.dat
y voy a pintar con gnuplot porque con R es demasiado complicado,
- tau_plot.plt
set terminal postscript eps color enhanced size 10,7 "Times-Roman, 32"
set output "braak_1_comp.eps"
m(x) = a*x +b
fit m(x) "tau_ropi.dat" u 9:2 via a,b
unset key
f(x) = x
set xlabel "PI2620 SUVR" font "Times-Roman, 24"
set ylabel "RO948 SUVR" font "Times-Roman, 24"
set title "BRAAK 1 ROI"
plot [0.5:1.5][0.5:1.5] "tau_ropi.dat" u 9:2 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
set title "BRAAK 2 ROI"
set output "braak_2_comp.eps"
fit m(x) "tau_ropi.dat" u 10:3 via a,b
plot [0.2:1.1][0.2:1.1] "tau_ropi.dat" u 10:3 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
set title "BRAAK 3 ROI"
set output "braak_3_comp.eps"
fit m(x) "tau_ropi.dat" u 11:4 via a,b
plot [0.3:1.8][0.3:1.8] "tau_ropi.dat" u 11:4 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
set title "BRAAK 4 ROI"
set output "braak_4_comp.eps"
fit m(x) "tau_ropi.dat" u 12:5 via a,b
plot [0.6:2.0][0.6:2.0] "tau_ropi.dat" u 12:5 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
set title "BRAAK 5 ROI"
set output "braak_5_comp.eps"
fit m(x) "tau_ropi.dat" u 13:6 via a,b
plot [0.5:1.8][0.5:1.8] "tau_ropi.dat" u 13:6 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
set title "BRAAK 6 ROI"
set output "braak_6_comp.eps"
fit m(x) "tau_ropi.dat" u 14:7 via a,b
plot [0.4:1.6][0.4:1.6] "tau_ropi.dat" u 14:7 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
esto puede extenderse a la combinacion de ROIs
Para intentar caracterizar el off-target binding he descomentado la parte en que se calcula la retencion en eroded WM,
#Hacer mascara de Eroded WM
$ptask{'output'} = $outdir.'/tau_ewm_'.$subject;
$ptask{'filename'} = $outdir.'/'.$subject.'_ewm.sh';
$ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_tref_ewm.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject;
$mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/ewm.nii.gz ';
send2slurm(\%ptask);
y entonces,
[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer ro948 -r meta -p f5cehbi
...
...
...
[osotolongo@brick03 f5cehbi]$ for x in `find working/ -name "*_unc.csv"`; do id=$(echo ${x} | sed 's/.*\/\(.*\)_.*/\1/'); ewm=$(grep "^3" ${x} | awk '{print $2}'); icgm=$(grep "^4" ${x} | awk '{print $2}'); echo "${id} ${ewm} ${icgm}"; done | sort -n > ewm_icgm_ro948.csv
[osotolongo@brick03 f5cehbi]$ cat ewm_icgm_ro948.csv
0065 29.982110 20.723816
0074 20.095760 19.340429
0083 16.295862 13.545850
0089 35.377898 42.095356
0090 36.817361 45.491038
0092 39.262153 46.518434
0094 15.278639 22.369623
0097 21.845054 24.733724
0115 25.884322 34.812070
0118 20.809293 16.964745
0121 40.504989 36.016138
0141 30.426140 28.183759
0142 14.925220 15.207533
0149 42.510422 25.144353
Las columnas son: ID, EWM, ICGM
[osotolongo@brick03 f5cehbi]$ sed '1iid,ewm,icgm' ewm_icgm_ro948.csv | sed 's/ /,/g' > ewm_icgm_ro948_fixed.csv
Ahora hacemos lo mismo para pi2620 y,
[osotolongo@brick03 f5cehbi]$ for x in `find working/ -name "*_unc.csv"`; do id=$(echo ${x} | sed 's/.*\/\(.*\)_.*/\1/'); ewm=$(grep "^3" ${x} | awk '{print $2}'); icgm=$(grep "^4" ${x} | awk '{print $2}'); echo "${id} ${ewm} ${icgm}"; done | sort -n > ewm_icgm_pi2620.csv
[osotolongo@brick03 f5cehbi]$ awk '{print $0" r0948"}' ewm_icgm_ro948.csv | sed 's/ /,/g;1iid,ewm,icgm,tracer' > ewm_icgm_1.csv
[osotolongo@brick03 f5cehbi]$ awk '{print $0" pi2620"}' ewm_icgm_pi2620.csv | sed 's/ /,/g' > ewm_icgm_2.csv
[osotolongo@brick03 f5cehbi]$ cat ewm_icgm_1.csv ewm_icgm_2.csv > ewm_icgm.csv
Ahora, aqui hay una cosa rara. Estas son las regiones de referencia y supuestamente la relacion entre ambas deberia ser cercana a 1. Al ser menor que 1 implica que se esta reteniendo el trazador en el cerebelo.
Problemas y Workarounds
Las imagenes PET estan hechas muy ajustadas, practicamente desde la base del craneo (ver imagenes abajo). Las imagenes MRI son ma generosas y tienen mas de la mitad del cuello. El algoritmo de registro (script antsRegistrationSyNQuick.sh) tiene problemas con esto y muchas veces no logra registrar correctamente.
La solucion pasa por recortar el cuello de las MRI. Esto es tan sencillo como convertir el 3D en un numero determinado de imagenes 2D a lo largo del eje Z, eliminar las inferiores y reemsamblar el 3D con lo que queda.
- cutslices.sh
Usage() {
echo ""
echo "Usage: cutslices.sh <input> <output> <n>"
echo ""
echo "Cut <n> slices from <input> and save result as <output>"
echo ""
echo "You must have FSL installed in order to run this script"
echo ""
echo ""
exit 1
}
[ "$3" = "" ] && Usage
debug=0
# }}}
# {{{ parse main arguments
in=`${FSLDIR}/bin/remove_ext $1`
shift
out=`${FSLDIR}/bin/remove_ext $1`
shift
n=$1
shift
${FSLDIR}/bin/fslreorient2std ${in} ${in}_tmp
${FSLDIR}/bin/fslsplit ${in}_tmp ${in}_2crop -z
for ((x=0;x<=$n;x+=1)); do a=`printf "%04d" $x`; ${FSLDIR}/bin/imrm ${in}_2crop$a; done
${FSLDIR}/bin/fslmerge -z $out ${in}_2crop*
if [ $debug = 0 ] ; then
rm ${in}_tmp*
rm ${in}_2crop*
fi
Para esto necesito el numero de slices que tengo que eliminar de la MRI. Esto requiere una inspeccion visual previa y mirar donde esta la base del craneo (aprox) de cada MRI. Si fala el registro se puede ir variando el nuemro de slices hasta que mejor se ajuste.
Lo que he hecho es que desde el wraper se lee un archivo con el numero de slices a recortar,
[osotolongo@brick03 acenip]$ head /old_nas/f5cehbi/cuts_pi2620.csv
0065,80
0070,70
0074,75
0075,80
0077,60
0081,60
0083,70
0085,80
0086,80
0087,75
El archivo se especifica en el input del wrapper,
if (/^-s/) {$slices_file = shift; chomp($slices_file);}
y se lee apropiadamente,
open IDF, "<$cut_file" or die "No file with cutting slices\n";
while (<IDF>) {
my ($id, $slice) = /(\d*),(\d*)/;
$slices{$id} = $slice;
}
close IDF;
entonces se pasa la info de slices a recortar al script de registro,
if(exists $slices{$subject} and $slices{$subject}){
$ptask{'command'} = $ENV{'PIPEDIR'}."/bin/tau_reg.sh ".$study." ".$subject." ".$w_dir." ".$spet{'tau'}." ".$slices{$subject};
}else{
$ptask{'command'} = $ENV{'PIPEDIR'}."/bin/tau_reg.sh ".$study." ".$subject." ".$w_dir." ".$spet{'tau'}." 0";
}
dentro de este script se llama a cutslices.sh y se calcula la info de registro del T1 recortado a la original pues esta info debera utilizarse luego.
if [[ $slices != 0 ]]; then
${PIPEDIR}/bin/cutslices.sh ${wdir}/${id}_struc.nii.gz ${wdir}/${id}_struc_croped.nii.gz ${slices};
${ANTS_PATH}/antsRegistrationSyNQuick.sh -d 3 -j 1 -f ${wdir}/${id}_struc.nii.gz -m ${wdir}/${id}_struc_croped.nii.gz -t t -o ${wdir}/${id}_cropedToFixed_;
fi
Entonces cada imagen PET (6x5min) se registra al T1 recortado y se pasa despues a espacio original,
${ANTS_PATH}/antsRegistrationSyNQuick.sh -d 3 -f ${wdir}/${id}_struc_croped.nii.gz -m ${x%.nii.gz}_mod.nii.gz -t a -o ${x%.nii.gz}_movingToCroped_;
${ANTS_PATH}/antsApplyTransforms -d 3 -r ${wdir}/${id}_struc_croped.nii.gz -i ${x} -t ${x%.nii.gz}_movingToCroped_0GenericAffine.mat -o ${x%.nii.gz}_reg_croped.nii.gz;
${ANTS_PATH}/antsApplyTransforms -d 3 -r ${wdir}/${id}_struc.nii.gz -i ${x%.nii.gz}_reg_croped.nii.gz -t ${wdir}/${id}_cropedToFixed_0GenericAffine.mat -o ${x%.nii.gz}_reg.nii.gz;
Todo el codigo en github, tanto para el wrapper como para el script de registro.
Las imagenes con poca informacion (con poca o ninguna captacion) tienes problemas para registrarse correctamente a espacio nativo MRI. Esto es particularmente cierto para PI2620, donde la señal es muy baja. Para RO948 la señal es mas fuerte y se aprecian mejor los matices, dando mas informacion y facilitando el registro.
La solución temporal para esto puede se resaltar las areas con menos señal mientras las areas con menor señal se mantienen mas o menos inalteradas (enfasis en + o -). Por ejemplo, si elevo la imagen a 3/2,
fslmaths ${x} -sqr -mul ${x} -sqrt ${x%.nii.gz}_mod.nii.gz
${ANTS_PATH}/antsRegistrationSyNQuick.sh -d 3 -f ${wdir}/${id}_struc_croped.nii.gz -m ${x%.nii.gz}_mod.nii.gz -t a -o ${x%.nii.gz}_movingToCroped_;
registro esa imagen modificada y luego uso esta info para registrar la imagen original.
Lo que voy a hacer es aplicar el pet original (ya registrado) como una mascara y excluir las regiones del cerebelo que no esten. Pierdo una parte del cerebelo inferior pero se supone que sea homogeneo, blah, blah, blah.
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/icgm_full.nii.gz -mas ${tmp_dir}/${subject}_tau_mask.nii.gz ${tmp_dir}/rois/icgm_masked.nii.gz
Esto es bastante jodido. Lo que hare es quitar el borde (una especie de eroded icgm)
y esperar que la contaminacion exterior se elimine.
${FSLDIR}/bin/fslmaths ${tmp_dir}/rois/icgm_masked.nii.gz -kernel gauss 3.3973 -fmean -thr 0.7 -bin ${tmp_dir}/rois/icgm.nii.gz
FSGA
Para hacer el archivo principal del analisis SBM, voy a añadir las covariables. La edad y el genero puedo sacarlos de XNAT. Luego he de ir pegando archivos,
$ join -t, f5cehbi_tau_suvr_pi2620_unc_icgm.csv f5cehbi_age_gender_recoded_sorted.csv > f5cehbi_unc_covs.csv
para lograr un archivo con este formato,
[osotolongo@brick03 fsga]$ head f5cehbi_unc_covs.csv
Subject,braak_1,braak_2,braak_3,braak_4,braak_5,braak_6,extra,AGE,Gender
0065,1.89537957199309,1.36073293452894,1.92020844877494,1.88817601591982,2.08289077909769,1.75747791947603,1.4028012031543,74.7,2
0070,0.935908461890794,0.872846276787541,1.04387128228284,1.09844634695409,1.15393097798958,1.01620313973646,0.822044288302219,66.4,1
0074,1.09330958102051,0.621896466311178,0.87135791910626,0.979736826194312,1.00314766091566,0.972597967319522,0.761565234181486,72.7,2
0075,0.920501428843966,0.899281738999484,1.12014603844927,1.29233468764307,1.24634622596148,1.10127253846343,0.911070304464376,71.4,1
0077,1.01453341640456,0.895199397062266,1.07513513726626,1.20334944569985,1.24142400884331,1.10915864591335,0.806642437912619,74.1,1
0081,0.893958373845497,1.01441844592613,1.09955643311136,1.13720575703509,1.26856292004374,1.25774989729188,0.761363147503239,70.6,1
0083,1.11738653957562,0.863383448881703,1.09900763626122,1.12559774574884,1.31892444435365,1.30752258555954,0.933675941180505,75.8,2
0085,1.13475873873386,1.64498323110215,1.27987508437459,1.40300220404692,1.6178920664913,1.3331126537738,1.31511774746641,72.8,1
0086,1.50160642584409,1.36626632369982,1.44629967994608,1.32724721966787,1.45714468139895,1.3348357481023,1.1177271105864,74.3,1
y ahora, el listado de variables es,
[osotolongo@brick03 fsga]$ head -n 1 f5cehbi_unc_covs.csv | sed 's/,/\n/g' | cat -n
1 Subject
2 braak_1
3 braak_2
4 braak_3
5 braak_4
6 braak_5
7 braak_6
8 extra
9 AGE
10 Gender
asi que para braak 1 hacemos,
[osotolongo@brick03 fsga]$ tail -n +2 f5cehbi_unc_covs.csv | awk -F',' '{print "Input f5cehbi_"$1" Main "$9" "$10" "$2}' | sed '1iVariables Age Gender SUVR' > braak1_unc_body.dat
[osotolongo@brick03 fsga]$ head braak1_unc_body.dat
Variables Age Gender SUVR
Input f5cehbi_0065 Main 74.7 2 1.89537957199309
Input f5cehbi_0070 Main 66.4 1 0.935908461890794
Input f5cehbi_0074 Main 72.7 2 1.09330958102051
Input f5cehbi_0075 Main 71.4 1 0.920501428843966
Input f5cehbi_0077 Main 74.1 1 1.01453341640456
Input f5cehbi_0081 Main 70.6 1 0.893958373845497
Input f5cehbi_0083 Main 75.8 2 1.11738653957562
Input f5cehbi_0085 Main 72.8 1 1.13475873873386
Input f5cehbi_0086 Main 74.3 1 1.50160642584409
[osotolongo@brick03 fsga]$ cat header_braak1.txt braak1_unc_body.dat > braak1_unc.fsgd
[osotolongo@brick03 fsga]$ head braak1_unc.fsgd
GroupDescriptorFile 1
Title FACEHBIV5_Braak_1
Class Main
Variables Age Gender SUVR
Input f5cehbi_0065 Main 74.7 2 1.89537957199309
Input f5cehbi_0070 Main 66.4 1 0.935908461890794
Input f5cehbi_0074 Main 72.7 2 1.09330958102051
Input f5cehbi_0075 Main 71.4 1 0.920501428843966
Input f5cehbi_0077 Main 74.1 1 1.01453341640456
Input f5cehbi_0081 Main 70.6 1 0.893958373845497
