neuroimagen:centiloid
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revisionLast revisionBoth sides next revision | ||
neuroimagen:centiloid [2019/04/12 08:36] – [Implementando pipeline Centiloid] osotolongo | neuroimagen:centiloid [2019/04/23 10:01] – [Implementando pipeline Centiloid] osotolongo | ||
---|---|---|---|
Line 144: | Line 144: | ||
</ | </ | ||
- | Ahora, las imagenes | + | Ahora, las imagenes |
<code bash fbbtemp_reg.sh> | <code bash fbbtemp_reg.sh> | ||
Line 213: | Line 213: | ||
</ | </ | ||
- | Lo voy a hacer con // | + | A ver que tal ese modelo, |
< | < | ||
+ | > d <- read.csv(" | ||
+ | > m <- lm(d$Global ~ d$WC_suvr) | ||
+ | > summary(m) | ||
+ | |||
+ | Call: | ||
+ | lm(formula = d$Global ~ d$WC_suvr) | ||
+ | |||
+ | Residuals: | ||
+ | Min 1Q Median | ||
+ | -0.121024 -0.042402 | ||
+ | |||
+ | Coefficients: | ||
+ | Estimate Std. Error t value Pr(> | ||
+ | (Intercept) | ||
+ | d$WC_suvr | ||
+ | --- | ||
+ | Signif. codes: | ||
+ | |||
+ | Residual standard error: 0.07003 on 22 degrees of freedom | ||
+ | Multiple R-squared: | ||
+ | F-statistic: | ||
+ | |||
+ | </ | ||
+ | Voy a hacer con //gnuplot// que es mas facil, | ||
+ | |||
+ | < | ||
+ | gnuplot> f(x) = m*x +n | ||
gnuplot> fit f(x) " | gnuplot> fit f(x) " | ||
iter chisq | iter chisq | ||
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
iter chisq | iter chisq | ||
- | After 5 iterations the fit converged. | + | After 4 iterations the fit converged. |
- | final sum of squares of residuals : 0.134007 | + | final sum of squares of residuals : 0.0793291 |
- | rel. change during last iteration : -1.72206e-08 | + | rel. change during last iteration : -3.5012e-11 |
degrees of freedom | degrees of freedom | ||
- | rms of residuals | + | rms of residuals |
- | variance of residuals (reduced chisquare) = WSSR/ | + | variance of residuals (reduced chisquare) = WSSR/ |
Final set of parameters | Final set of parameters | ||
======================= | ======================= | ||
- | m | + | m |
- | n | + | n = 0.0186989 |
correlation matrix of the fit parameters: | correlation matrix of the fit parameters: | ||
Line 246: | Line 272: | ||
{{: | {{: | ||
- | Y con pendiente de $1.08 \pm 0.04$ creo que estamos bien. | ||
- | Resumen: Hemos validado correctamente | + | Tengo un ajuste con $r^2 = 0.97$ y una pendiente de $0.85 \pm 0.02$. La diferencia en los procedimientos esta en que en nuestra pipeline, siguiendo |
+ | Vamos a comprobarlo. Cambio la normalizacion para que normalice por el cerebelo completo y de paso incluyo todos los sujetos, | ||
+ | |||
+ | <code bash> | ||
+ | [osotolongo@detritus centiloid]$ awk -F ";" | ||
+ | [osotolongo@detritus centiloid]$ awk -F ";" | ||
+ | [osotolongo@detritus centiloid]$ join calcs_wcb.dat reference.dat | ||
+ | < | ||
+ | |||
+ | Reviso el modelo en R, | ||
+ | |||
+ | < | ||
+ | > d <- read.csv(" | ||
+ | > m <- lm(d$Global ~ d$WC_suvr) | ||
+ | > summary(m) | ||
+ | |||
+ | Call: | ||
+ | lm(formula = d$Global ~ d$WC_suvr) | ||
+ | |||
+ | Residuals: | ||
+ | Min 1Q Median | ||
+ | -0.073212 -0.028196 | ||
+ | |||
+ | Coefficients: | ||
+ | Estimate Std. Error t value Pr(> | ||
+ | (Intercept) | ||
+ | d$WC_suvr | ||
+ | --- | ||
+ | Signif. codes: | ||
+ | |||
+ | Residual standard error: 0.04017 on 32 degrees of freedom | ||
+ | Multiple R-squared: | ||
+ | F-statistic: | ||
+ | </ | ||
+ | y miro el ajuste entre los datos en gnuplot, | ||
+ | < | ||
+ | gnuplot> f(x) = m*x +n | ||
+ | gnuplot> fit f(x) " | ||
+ | |||
+ | After 4 iterations the fit converged. | ||
+ | final sum of squares of residuals : 0.0531383 | ||
+ | rel. change during last iteration : -3.75201e-12 | ||
+ | |||
+ | degrees of freedom | ||
+ | rms of residuals | ||
+ | variance of residuals (reduced chisquare) = WSSR/ | ||
+ | |||
+ | Final set of parameters | ||
+ | ======================= | ||
+ | |||
+ | m = 1.00834 | ||
+ | n = -0.0311848 | ||
+ | |||
+ | |||
+ | correlation matrix of the fit parameters: | ||
+ | |||
+ | | ||
+ | m | ||
+ | n -0.965 | ||
+ | </ | ||
+ | |||
+ | Tengo una pendiente de $1.01 \pm 0.02$ y un ajuste con $R^2 = 0.9878$, lo cual esta muy bien. Este resultado valida el metodo de obtencion de SUVR de nuestra pipeline. | ||
+ | |||
+ | {{: | ||
+ | |||
+ | **Nota:** El SUVR se debe calcular tomando como referencia la materia gris del cerebelo si los valores se van a comparar con datos de ADNI ya que esta es la metodologia que usan pero si se va a calcular centiloides debe tomarse como referencia el cerebelo completo. | ||
==== Implementando pipeline Centiloid ===== | ==== Implementando pipeline Centiloid ===== | ||
- | Con el objetivo de compara con otros estudios hemos de implementar el metodo usando las plantillas originales de [[https:// | + | Con el objetivo de compara con otros estudios hemos de implementar el metodo usando las plantillas originales de [[https:// |
+ | |||
+ | De cualquier manera voy a reescribir el proceso de registro para integrarlo en el cluster, | ||
+ | |||
+ | <code perl fbb_correct.pl> | ||
+ | # | ||
+ | |||
+ | use strict; use warnings; | ||
+ | use Data::Dump qw(dump); | ||
+ | use File:: | ||
+ | use File:: | ||
+ | |||
+ | use NEURO qw(populate check_subj load_study print_help get_petdata get_pair achtung shit_done check_or_make); | ||
+ | |||
+ | #my $study; | ||
+ | my $cfile; | ||
+ | my $sok = 0; | ||
+ | |||
+ | @ARGV = (" | ||
+ | while (@ARGV and $ARGV[0] =~ /^-/) { | ||
+ | $_ = shift; | ||
+ | last if /^--$/; | ||
+ | #if (/^-e/) { $study = shift; chomp($study); | ||
+ | if (/^-cut/) { $cfile = shift; chomp($cfile); | ||
+ | if (/^-h/) { print_help $ENV{' | ||
+ | } | ||
+ | my $study = shift; | ||
+ | unless ($study) { print_help $ENV{' | ||
+ | |||
+ | my %std = load_study($study); | ||
+ | |||
+ | my $subj_dir = $ENV{' | ||
+ | |||
+ | #Run this script ONLY on " | ||
+ | #or change the paths to the correct ones | ||
+ | |||
+ | my $w_dir=$std{' | ||
+ | #my $pet_dir = $spect? | ||
+ | my $pet_dir = $std{' | ||
+ | #my $petnc_dir = $std{' | ||
+ | my $db = $std{' | ||
+ | my $data_dir=$std{' | ||
+ | my $outdir = " | ||
+ | check_or_make($outdir); | ||
+ | |||
+ | my %plist = populate($db); | ||
+ | my @ok_pets; | ||
+ | |||
+ | print " | ||
+ | my %pets; | ||
+ | |||
+ | if ($cfile){ | ||
+ | my %cuts = get_pair($data_dir."/" | ||
+ | foreach my $cut (keys %cuts){ | ||
+ | if(grep {/$cut/} %plist){ | ||
+ | $pets{$cut} = $plist{$cut}; | ||
+ | } | ||
+ | } | ||
+ | }else{ | ||
+ | %pets = %plist; | ||
+ | } | ||
+ | |||
+ | my $err_file = $data_dir."/ | ||
+ | open ERRF, "> | ||
+ | |||
+ | foreach my $subject (sort keys %pets){ | ||
+ | my @names = find(file => ' | ||
+ | dump @names; | ||
+ | if(@names){ | ||
+ | my @mri = find(file => ' | ||
+ | if(@mri){ | ||
+ | push @ok_pets, $subject; | ||
+ | }else{ | ||
+ | print ERRF " | ||
+ | } | ||
+ | } | ||
+ | } | ||
+ | |||
+ | foreach my $subject (@ok_pets){ | ||
+ | #my @mri = find(file => ' | ||
+ | my $orderfile = $outdir.'/' | ||
+ | my $order = $ENV{' | ||
+ | open ORD, "> | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD "srun $order\n"; | ||
+ | close ORD; | ||
+ | my $jobid = `sbatch $orderfile`; | ||
+ | $jobid = ( split ' ', $jobid )[ -1 ]; | ||
+ | $order = $ENV{' | ||
+ | $orderfile = $outdir.'/' | ||
+ | open ORD, "> | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD "srun $order\n"; | ||
+ | close ORD; | ||
+ | my $xorder = ' | ||
+ | system(" | ||
+ | } | ||
+ | |||
+ | my $order = $ENV{' | ||
+ | my $orderfile = $outdir.'/ | ||
+ | open ORD, "> | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD "srun $order\n"; | ||
+ | close ORD; | ||
+ | my $xorder = ' | ||
+ | exec($xorder); | ||
+ | </ | ||
+ | |||
+ | **Y ahora la parte nueva ** | ||
- | Voy a hacer un script nuevo para esto, reutilizando el antiguo. | + | Primero hago un script sencillo que registre con ANTS al espacio MNI, |
<code bash fbb2std.sh> | <code bash fbb2std.sh> | ||
Line 266: | Line 479: | ||
echo "I need the FBB image at MNI space" | echo "I need the FBB image at MNI space" | ||
- | ANTS 3 -m CC[${FSLDIR}/ | + | ANTS 3 -m CC[${FSLDIR}/ |
- | WarpImageMultiTransform 3 ${wdir}/ | + | WarpImageMultiTransform 3 ${wdir}/ |
</ | </ | ||
- | y despues hago un //wrapper// que lo lance en paralelo y haga las metricas, | + | la relacion reportada |
- | <code perl parallel_fbb_cl_metrics.pl> | + | |
- | # | + | |
+ | <code perl fbb_cl_masks.pl> | ||
+ | # | ||
+ | # Copyright 2019 O. Sotolongo < | ||
use strict; use warnings; | use strict; use warnings; | ||
Line 282: | Line 496: | ||
use File:: | use File:: | ||
use File:: | use File:: | ||
- | use Parallel:: | ||
- | my %ROI_Comps | + | my $ref_roi |
- | " | + | my $ctx_roi |
- | " | + | |
- | " | + | |
- | " | + | |
- | " | + | |
- | ); | + | |
- | my $roi_paths = $ENV{' | + | my $roi_paths = $ENV{' |
- | my $attach = 1; | + | |
- | my $reduce = 0; | + | |
my $withstd = 0; | my $withstd = 0; | ||
my $cfile; | my $cfile; | ||
Line 303: | Line 509: | ||
$_ = shift; | $_ = shift; | ||
last if /^--$/; | last if /^--$/; | ||
- | if (/^-l$/) { $attach = 0;} | ||
- | if (/^-r$/) { $reduce = 1;} | ||
- | if (/^-std$/) {$withstd = 1;} | ||
if (/^-cut/) { $cfile = shift; chomp($cfile); | if (/^-cut/) { $cfile = shift; chomp($cfile); | ||
- | if (/^-h$/) { print_help $ENV{' | ||
} | } | ||
Line 315: | Line 517: | ||
my $w_dir=$std{' | my $w_dir=$std{' | ||
my $data_dir=$std{' | my $data_dir=$std{' | ||
- | my $max_processes = 20; | ||
- | # Redirect ouput to logfile (do it only when everything is fine) | ||
- | #my $debug = " | ||
- | #open STDOUT, "> | ||
- | #open STDERR, ">& | ||
- | |||
- | my $pm = new Parallel:: | ||
our %subjects; | our %subjects; | ||
- | $pm-> | ||
- | sub { my ($pid, $exit_code, $ident, $exit_signal, | ||
- | foreach my $tag (sort keys %{$data}){ | ||
- | $subjects{$ident}{$tag}=${$data}{$tag} | ||
- | } | ||
- | } | ||
- | ); | ||
my @plist = find(file => ' | my @plist = find(file => ' | ||
my $ofile = $data_dir."/" | my $ofile = $data_dir."/" | ||
Line 348: | Line 536: | ||
@pets = @plist; | @pets = @plist; | ||
} | } | ||
- | |||
foreach my $pet (sort @pets){ | foreach my $pet (sort @pets){ | ||
Line 357: | Line 544: | ||
} | } | ||
} | } | ||
+ | my %sdata; | ||
foreach my $subject (sort keys %subjects){ | foreach my $subject (sort keys %subjects){ | ||
my $care; | my $care; | ||
my $norm; | my $norm; | ||
my $dg = $subjects{$subject}{' | my $dg = $subjects{$subject}{' | ||
- | my @care; | + | # Apply masks to FBB |
- | my %sdata; | + | my $roi_mask = $roi_paths.$ctx_roi; |
- | $pm-> | + | my $order = "fslstats ".$w_dir."/ |
- | + | my $ctx = qx/$order/; | |
- | # Get FBB image into MNI space | + | my $roi_mask = $roi_paths.$ctx_roi; |
- | my $order = "fbb2std.sh " | + | $order = " |
- | print | + | print " |
- | system($order); | + | my $norm = qx/ |
- | # Apply masks to FBB | + | $subjects{$subject}{$npf.' |
- | foreach my $npf (sort keys %ROI_Comps){ | + | |
- | my $roi_mask = $roi_paths.$ROI_Comps{$npf}; | + | |
- | # get mean and std for mask | + | |
- | $order = " | + | |
- | print " | + | |
- | (my $mean, my $std) = map{/ | + | |
- | $sdata{$npf.' | + | |
- | $sdata{$npf.' | + | |
- | } | + | |
- | $pm-> | + | |
- | # remove temp dir | + | |
- | #remove( \1, ($mdir)); | + | |
} | } | ||
- | $pm-> | ||
open OF, "> | open OF, "> | ||
- | + | print OF " | |
- | print OF " | + | |
- | + | ||
- | foreach my $npf (sort keys %ROI_Comps){ | + | |
- | if($withstd){ | + | |
- | print OF ";$npf"," | + | |
- | }else{ | + | |
- | print OF "; | + | |
- | } | + | |
- | } | + | |
print OF " | print OF " | ||
foreach my $subject (sort keys %subjects){ | foreach my $subject (sort keys %subjects){ | ||
print OF " | print OF " | ||
- | foreach my $npf (sort keys %ROI_Comps){ | + | my $mean = $subjects{$subject}{$npf.' |
- | my $mean = $subjects{$subject}{$npf.' | + | print OF "; |
- | my $std = $subjects{$subject}{$npf.' | + | print OF ";", |
- | if($withstd){ | + | |
- | print OF "; | + | |
- | print OF ";", | + | |
- | }else{ | + | |
- | print OF "; | + | |
- | print OF ";", | + | |
- | } | + | |
- | } | + | |
print OF " | print OF " | ||
} | } | ||
close OF; | close OF; | ||
+ | </ | ||
- | my $zfile = $ofile.' | + | Ahora hago un wrapper que lance todo, |
- | system(" | + | |
- | if ($attach){ | + | <code perl fbb_cl_metrics.pl> |
- | shit_done basename($ENV{_}), $study, $zfile; | + | # |
+ | # Copyright 2019 O. Sotolongo < | ||
+ | use strict; use warnings; | ||
+ | |||
+ | use File:: | ||
+ | use NEURO qw(print_help get_pair load_study achtung shit_done get_lut check_or_make centiloid_fbb); | ||
+ | use Data::Dump qw(dump); | ||
+ | use File:: | ||
+ | use File:: | ||
+ | |||
+ | my $wcut = 0; | ||
+ | my $cfile; | ||
+ | |||
+ | @ARGV = (" | ||
+ | |||
+ | while (@ARGV and $ARGV[0] =~ /^-/) { | ||
+ | $_ = shift; | ||
+ | last if /^--$/; | ||
+ | | ||
+ | if (/^-h$/) { print_help $ENV{' | ||
+ | } | ||
+ | |||
+ | my $study = shift; | ||
+ | unless | ||
+ | my %std = load_study($study); | ||
+ | my $w_dir=$std{' | ||
+ | my $data_dir=$std{' | ||
+ | my $outdir = " | ||
+ | check_or_make($outdir); | ||
+ | our %subjects; | ||
+ | |||
+ | my @plist = find(file => ' | ||
+ | my $ofile = $data_dir."/" | ||
+ | my $patt = '([A-Z,a-z]{1, | ||
+ | |||
+ | my @pets; | ||
+ | |||
+ | if ($cfile){ | ||
+ | my %cuts = get_pair($data_dir."/" | ||
+ | foreach my $cut (keys %cuts){ | ||
+ | if(grep {/$cut/} @plist){ | ||
+ | $pets[$cut] = $plist[$cut]; | ||
+ | } | ||
+ | } | ||
}else{ | }else{ | ||
- | achtung basename($ENV{_}), | + | @pets = @plist; |
} | } | ||
+ | |||
+ | foreach my $pet (sort @pets){ | ||
+ | (my $dg,my $subject) = $pet =~ /$patt/; | ||
+ | if($subject){ | ||
+ | $subjects{$subject}{' | ||
+ | $subjects{$subject}{' | ||
+ | } | ||
+ | } | ||
+ | |||
+ | foreach my $subject (sort keys %subjects){ | ||
+ | my $norm; | ||
+ | my $dg = $subjects{$subject}{' | ||
+ | my $subj = $dg.$subject; | ||
+ | #Making sbatch scripts | ||
+ | # Get FBB image into MNI space | ||
+ | my $order = $ENV{' | ||
+ | my $orderfile = $outdir.'/' | ||
+ | open ORD, "> | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD "srun $order\n"; | ||
+ | close ORD; | ||
+ | system(" | ||
+ | } | ||
+ | my $order = $ENV{' | ||
+ | my $orderfile = $outdir.'/ | ||
+ | open ORD, "> | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD '# | ||
+ | print ORD "srun $order\n"; | ||
+ | close ORD; | ||
+ | my $xorder = ' | ||
+ | exec($xorder); | ||
+ | </ | ||
+ | |||
+ | |||
+ | |||
+ | <code bash> | ||
+ | [osotolongo@detritus centiloid]$ fbb_cl_metrics.pl centiloid | ||
+ | ..... | ||
+ | [osotolongo@detritus centiloid]$ squeue | ||
+ | JOBID PARTITION | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | </ | ||
+ | |||
+ | A ver esos resultados, | ||
+ | |||
+ | <code bash> | ||
+ | [osotolongo@detritus centiloid]$ sed 's/;/ /g; s/ / /g' centiloid_fbb_cl.csv > centiloid_fbb_cl.dat | ||
+ | [osotolongo@detritus centiloid]$ join centiloid_fbb_cl.dat reference.dat > toreview.dat | ||
+ | </ | ||
+ | |||
+ | Compruebo el modelo, | ||
+ | |||
+ | < | ||
+ | > d <- read.csv(" | ||
+ | > m <- lm(d$SUVR ~ d$WC_suvr) | ||
+ | > summary(m) | ||
+ | |||
+ | Call: | ||
+ | lm(formula = d$SUVR ~ d$WC_suvr) | ||
+ | |||
+ | Residuals: | ||
+ | Min 1Q Median | ||
+ | -0.077423 -0.020244 -0.007479 | ||
+ | |||
+ | Coefficients: | ||
+ | Estimate Std. Error t value Pr(> | ||
+ | (Intercept) | ||
+ | d$WC_suvr | ||
+ | --- | ||
+ | Signif. codes: | ||
+ | |||
+ | Residual standard error: 0.03411 on 32 degrees of freedom | ||
+ | Multiple R-squared: | ||
+ | F-statistic: | ||
+ | </ | ||
+ | y el ajuste, | ||
+ | < | ||
+ | gnuplot> f(x) = m*x +n | ||
+ | gnuplot> fit f(x) " | ||
+ | iter chisq | ||
+ | 0 3.4203670136e+01 | ||
+ | 1 3.2574055507e-01 | ||
+ | 2 4.1098565887e-02 | ||
+ | 3 4.1072773521e-02 | ||
+ | 4 4.1072773521e-02 | ||
+ | iter chisq | ||
+ | |||
+ | After 4 iterations the fit converged. | ||
+ | final sum of squares of residuals : 0.0410728 | ||
+ | rel. change during last iteration : -5.9486e-12 | ||
+ | |||
+ | degrees of freedom | ||
+ | rms of residuals | ||
+ | variance of residuals (reduced chisquare) = WSSR/ | ||
+ | |||
+ | Final set of parameters | ||
+ | ======================= | ||
+ | m = 1.04544 | ||
+ | n = -0.0618459 | ||
+ | |||
+ | correlation matrix of the fit parameters: | ||
+ | m n | ||
+ | m | ||
+ | n -0.967 | ||
</ | </ | ||
+ | Con lo que tengo una pendiente de $1.05 \pm 0.02$ y una $R^2 = 0.9906$ que esta bastante bien. | ||
+ | {{: | ||
neuroimagen/centiloid.txt · Last modified: 2020/08/04 10:58 by 127.0.0.1