User Tools

Site Tools


neuroimagen:centiloid

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision
Previous revision
Last revisionBoth sides next revision
neuroimagen:centiloid [2019/04/12 08:36] – [Implementando pipeline Centiloid] osotolongoneuroimagen:centiloid [2019/04/23 10:01] – [Implementando pipeline Centiloid] osotolongo
Line 144: Line 144:
 </code> </code>
  
-Ahora, las imagenes PEt vieen en formato distinto al usado habitualmente (1 solo slice de 20 min) por lo que hay que retocar el script //fbbtemp_reg.sh// un poco para añadirle la posibilidad de proesar esto correctamente.+Ahora, las imagenes PET vienen en formato distinto al usado habitualmente (1 solo slice de 20 min) por lo que hay que retocar el script //fbbtemp_reg.sh// un poco para añadirle la posibilidad de proesar esto correctamente.
  
 <code bash fbbtemp_reg.sh> <code bash fbbtemp_reg.sh>
Line 213: Line 213:
 </code> </code>
  
-Lo voy a hacer con //gnuplot// que es mas facil,+A ver que tal ese modelo,
  
 <code> <code>
 +> d <- read.csv("toreview.dat", sep = " ", header = TRUE)
 +> m <- lm(d$Global ~ d$WC_suvr)
 +> summary(m)
 +
 +Call:
 +lm(formula = d$Global ~ d$WC_suvr)
 +
 +Residuals:
 +      Min        1Q    Median        3Q       Max 
 +-0.121024 -0.042402  0.000411  0.047290  0.126656 
 +
 +Coefficients:
 +            Estimate Std. Error t value Pr(>|t|)    
 +(Intercept)  0.01999    0.05799   0.345    0.734    
 +d$WC_suvr    1.15167    0.03919  29.390   <2e-16 ***
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +Residual standard error: 0.07003 on 22 degrees of freedom
 +Multiple R-squared:  0.9752,    Adjusted R-squared:  0.974 
 +F-statistic: 863.8 on 1 and 22 DF,  p-value: < 2.2e-16
 +
 +</code>
 +Voy a hacer con //gnuplot// que es mas facil,
 +
 +<code>
 +gnuplot> f(x) = m*x +n
 gnuplot> fit f(x) "toreview.dat" u 2:3 via m,n gnuplot> fit f(x) "toreview.dat" u 2:3 via m,n
 iter      chisq       delta/lim  lambda                           iter      chisq       delta/lim  lambda                          
-   2.2640655922e+00   0.00e+00  1.03e+00    8.467408e-01   1.869890e-02 +   3.6936711672e+01   0.00e+00  1.41e+00    1.000000e+00   1.000000e+00 
-   1.5434204471e-01  -1.37e+06  1.03e-01    1.015137e+00   1.868981e-02 +   3.6992605848e-01  -9.88e+06  1.41e-01    6.105751e-01   4.583953e-01 
-   1.5172756904e-01  -1.72e+03  1.03e-02    1.021604e+00   1.342157e-02 +   7.9414981509e-02  -3.66e+05  1.41e-02    8.423745e-01   2.635659e-02 
-   1.3453722070e-01  -1.28e+04  1.03e-03    1.072796e+00  -7.770453e-02 +   7.9329087727e-02  -1.08e+02  1.41e-03    8.467400e-01   1.870027e-02 
-   1.3400723332e-01  -3.95e+02  1.03e-04    1.083478e+00  -9.671955e-02 +   7.9329087724e-02  -3.50e-06  1.41e-04    8.467408e-01   1.869890e-02
-   1.3400723101e-01  -1.72e-03  1.03e-05    1.083501e+00  -9.675931e-02+
 iter      chisq       delta/lim  lambda                           iter      chisq       delta/lim  lambda                          
  
-After iterations the fit converged. +After 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    (FIT_NDF)                        : 22 degrees of freedom    (FIT_NDF)                        : 22
-rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.0780464 +rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.0600489 
-variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00609124+variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00360587
  
 Final set of parameters            Asymptotic Standard Error Final set of parameters            Asymptotic Standard Error
 =======================            ========================== =======================            ==========================
-m               1.0835           +/- 0.03745      (3.456%) +m               0.846741         +/- 0.02881      (3.403%) 
-n               -0.0967593       +/- 0.0646       (66.76%)+n               = 0.0186989        +/- 0.0497       (265.8%)
  
 correlation matrix of the fit parameters: correlation matrix of the fit parameters:
Line 246: Line 272:
  
 {{:neuroimagen:fit.png|}} {{:neuroimagen:fit.png|}}
-Y con pendiente de $1.08 \pm 0.04$ creo que estamos bien. 
  
-Resumen: Hemos validado correctamente el metodo que se usapara sacar los SUVR.+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 el metodo de ADNI, se toma como region de referencia la materia gris del cerebelo y en el caso de GAAIN se toma todo el cerebelo. No obstante la correspondencia entre las mediciones es buena.
  
 +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 ";" '{print $1,$3}' centiloid_fbb_fs_suvr_roisi_wcb.csv > calcs_wcb.dat
 +[osotolongo@detritus centiloid]$ awk -F ";" '{print $1,$4}' FBB_suvr_centiloid.csv | sed 's/Y1/7/' > reference.dat
 +[osotolongo@detritus centiloid]$ join calcs_wcb.dat reference.dat 
 +<code>
 +
 +Reviso el modelo en R,
 +
 +<code>
 +> d <- read.csv("toreview.dat", sep = " ", header = TRUE)
 +> m <- lm(d$Global ~ d$WC_suvr)
 +> summary(m)
 +
 +Call:
 +lm(formula = d$Global ~ d$WC_suvr)
 +
 +Residuals:
 +      Min        1Q    Median        3Q       Max 
 +-0.073212 -0.028196  0.000229  0.026218  0.079672 
 +
 +Coefficients:
 +            Estimate Std. Error t value Pr(>|t|)    
 +(Intercept)  0.04626    0.02573   1.798   0.0817 .  
 +d$WC_suvr    0.98001    0.01894  51.731   <2e-16 ***
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +Residual standard error: 0.04017 on 32 degrees of freedom
 +Multiple R-squared:  0.9882,    Adjusted R-squared:  0.9878 
 +F-statistic:  2676 on 1 and 32 DF,  p-value: < 2.2e-16
 +</code>
 + y miro el ajuste entre los datos en gnuplot,
 +<code>
 +gnuplot> f(x) = m*x +n
 +gnuplot> fit f(x) "toreview.dat" u 2:3 via m,n
 +
 +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    (FIT_NDF)                        : 32
 +rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.0407501
 +variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00166057
 +
 +Final set of parameters            Asymptotic Standard Error
 +=======================            ==========================
 +
 +m               = 1.00834          +/- 0.01949      (1.933%)
 +n               = -0.0311848       +/- 0.02683      (86.03%)
 +
 +
 +correlation matrix of the fit parameters:
 +
 +                    n      
 +m               1.000 
 +n              -0.965  1.000 
 +</code>
 +
 +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. 
 +
 +{{:neuroimagen:wcbl_fit.png|}}
 +
 +**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://www.ncbi.nlm.nih.gov/pmc/articles/PMC4300247/|Klunk et. al.]]. La arte buena es que el preproceso y el registro al espacio de sujeto lo tenemos hecho. Solo habria que registrar al espacio MNI y sacar las metricas usando las plantillas.+Con el objetivo de compara con otros estudios hemos de implementar el metodo usando las plantillas originales de [[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4300247/|Klunk et. al.]]. La parte buena es que el preproceso y el registro al espacio de sujeto lo tenemos hecho. Solo habria que registrar al espacio MNI y sacar las metricas usando las plantillas. 
 + 
 +De cualquier manera voy a reescribir el proceso de registro para integrarlo en el cluster, 
 + 
 +<code perl fbb_correct.pl> 
 +#!/usr/bin/perl 
 + 
 +use strict; use warnings; 
 +use Data::Dump qw(dump); 
 +use File::Find::Rule; 
 +use File::Copy::Recursive qw(dirmove); 
 + 
 +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 = ("-h") unless @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{'PIPEDIR'}.'/doc/fbb_reg.hlp'; exit;} 
 +
 +my $study = shift; 
 +unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/fbb_reg.hlp'; exit;} 
 + 
 +my %std = load_study($study); 
 + 
 +my $subj_dir = $ENV{'SUBJECTS_DIR'}; 
 + 
 +#Run this script ONLY on "Detritus" 
 +#or change the paths to the correct ones 
 + 
 +my $w_dir=$std{'WORKING'}; 
 +#my $pet_dir = $spect?$std{'SPECT'}:$std{'PET-FDG'}; 
 +my $pet_dir = $std{'PET-FBB'}; 
 +#my $petnc_dir = $std{'FBB-NC'}; 
 +my $db = $std{'DATABASE'}; 
 +my $data_dir=$std{'DATA'}; 
 +my $outdir = "$std{'DATA'}/slurm"; 
 +check_or_make($outdir); 
 + 
 +my %plist = populate($db); 
 +my @ok_pets; 
 + 
 +print "Collecting needed files\n"; 
 +my %pets; 
 + 
 +if ($cfile){ 
 + my %cuts = get_pair($data_dir."/".$cfile); 
 + foreach my $cut (keys %cuts){ 
 + if(grep {/$cut/} %plist){ 
 + $pets{$cut} = $plist{$cut}; 
 +
 +
 +}else{ 
 + %pets = %plist; 
 +
 + 
 +my $err_file = $data_dir."/.fbb_errors.debug"; 
 +open ERRF, ">$err_file" or die "Could not open error file!"; 
 + 
 +foreach my $subject (sort keys %pets){ 
 + my @names = find(file => 'name' => "$pets{$subject}$subject*", '!name' => '*_*', 'name' => '*.nii.gz', in => $pet_dir);  
 + dump @names; 
 + if(@names){ 
 + my @mri = find(file => 'name' => "brain.mgz", in => $subj_dir."/".$study."_".$pets{$subject}.$subject."/mri"); 
 + if(@mri){ 
 + push @ok_pets, $subject; 
 + }else{ 
 + print ERRF "Subject $subject has no MRI brain extracted!!!\n"; 
 +
 + }  
 +
 + 
 +foreach my $subject (@ok_pets){ 
 + #my @mri = find(file => 'name' => $pets{$subject}.$subject."_brain.nii.gz", in => "$w_dir/"); 
 + my $orderfile = $outdir.'/'.$subject.'_get_subj.sh'; 
 + my $order = $ENV{'PIPEDIR'}."/bin/get_fs_subj.sh ".$study." ".$pets{$subject}.$subject." ".$w_dir; 
 + open ORD, ">$orderfile"; 
 + print ORD '#!/bin/bash'."\n"; 
 + print ORD '#SBATCH -J fbb_correct_'.$study."\n"; 
 + print ORD '#SBATCH --time=1:0:0'."\n"; #si no ha terminado en X horas matalo 
 + print ORD '#SBATCH --mail-type=FAIL,STAGE_OUT'."\n"; #no quieres que te mande email de todo 
 + print ORD '#SBATCH -o '.$outdir.'/fbb_get_subj-%j'."\n"; 
 + print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n"; 
 + print ORD "srun $order\n"; 
 + close ORD; 
 + my $jobid = `sbatch $orderfile`; 
 + $jobid = ( split ' ', $jobid )[ -1 ]; 
 + $order = $ENV{'PIPEDIR'}."/bin/fbbtemp_reg.sh ".$study." ".$pets{$subject}.$subject." ".$pet_dir." ".$w_dir; 
 + $orderfile = $outdir.'/'.$subject.'_fbb_reg.sh'; 
 + open ORD, ">$orderfile"; 
 + print ORD '#!/bin/bash'."\n"; 
 + print ORD '#SBATCH -J fbb_correct_'.$study."\n"; 
 + print ORD '#SBATCH --time=1:0:0'."\n"; #si no ha terminado en X horas matalo 
 + print ORD '#SBATCH --mail-type=FAIL,STAGE_OUT'."\n"; #no quieres que te mande email de todo 
 + print ORD '#SBATCH -o '.$outdir.'/fbb_reg-%j'."\n"; 
 + print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n"; 
 + print ORD "srun $order\n"; 
 + close ORD; 
 + my $xorder = 'sbatch --depend=afterok:'.$jobid.' '.$orderfile; 
 + system("$xorder"); 
 +
 + 
 +my $order = $ENV{'PIPEDIR'}."/bin/make_fbb_report.pl ".$study; 
 +my $orderfile = $outdir.'/fbb_report.sh'; 
 +open ORD, ">$orderfile"; 
 +print ORD '#!/bin/bash'."\n"; 
 +print ORD '#SBATCH -J fbb_correct_'.$study."\n"; 
 +print ORD '#SBATCH --time=2:0:0'."\n"; #si no ha terminado en X horas matalo 
 +print ORD '#SBATCH --mail-type=FAIL,END'."\n"; #email cuando termine o falle 
 +print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n"; 
 +print ORD '#SBATCH -o '.$outdir.'/fbb_rep-%j'."\n"; 
 +print ORD "srun $order\n"; 
 +close ORD; 
 +my $xorder = 'sbatch --dependency=singleton'.' '.$orderfile; 
 +exec($xorder); 
 +</code> 
 + 
 +**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,+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}/data/standard/MNI152_T1_1mm.nii.gz, ${wdir}/${id}_struc.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o ${wdir}/${id}_fbb_t1_mni.nii.gz +ANTS 3 -m CC[${FSLDIR}/data/standard/MNI152_T1_2mm.nii.gz, ${wdir}/${id}_struc.nii.gz, 1, 4] -r Gauss[0,3] -t Elast[1.5] -i 30x20x10 -o ${wdir}/${id}_fbb_t1_mni.nii.gz 
-WarpImageMultiTransform 3 ${wdir}/${id}_fbb.nii.gz ${wdir}/${id}_fbb_mni.nii.gz -R ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz ${wdir}/${id}_fbb_t1_mniWarp.nii.gz ${wdir}/${id}_fbb_t1_mniAffine.txt+WarpImageMultiTransform 3 ${wdir}/${id}_fbb.nii.gz ${wdir}/${id}_fbb_mni.nii.gz -R ${FSLDIR}/data/standard/MNI152_T1_2mm.nii.gz ${wdir}/${id}_fbb_t1_mniWarp.nii.gz ${wdir}/${id}_fbb_t1_mniAffine.txt
  
 </code> </code>
  
-y despues hago un //wrapper// que lo lance en paralelo y haga las metricas, +la relacion reportada en el paper de referencia es tomando como ROI de referencia el cerebelo completo. Asi que solo vale la pena calcular esa
-<code perl parallel_fbb_cl_metrics.pl> +
-#!/usr/bin/perl+
  
 +<code perl fbb_cl_masks.pl>
 +#!/usr/bin/perl
 +# Copyright 2019 O. Sotolongo <asqwerty@gmail.com>
 use strict; use warnings; use strict; use warnings;
    
Line 282: Line 496:
 use File::Remove 'remove'; use File::Remove 'remove';
 use File::Basename qw(basename); use File::Basename qw(basename);
-use Parallel::ForkManager; 
  
-my %ROI_Comps +my $ref_roi  "voi_WhlCbl_2mm.nii"; 
- "Cgray" => "voi_CerebGry_2mm.nii", +my $ctx_roi = "voi_ctx_2mm.nii";
- "WCbl" => "voi_WhlCbl_2mm.nii", +
- "WcblBS" => "voi_WhlCblBrnStm_2mm.nii", +
- "Pons" => "voi_Pons_2mm.nii", +
- "ctx" => "voi_ctx_2mm.nii"+
- );+
  
-my $roi_paths = $ENV{'PIPEDIR'}.'/lib/Centiloid_Std_VOI/nifti/1mm/'+my $roi_paths = $ENV{'PIPEDIR'}.'/lib/Centiloid_Std_VOI/nifti/2mm/';
-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{'PIPEDIR'}.'/doc/pet_metrics.hlp'; exit;} 
 } }
    
Line 315: Line 517:
 my $w_dir=$std{'WORKING'}; my $w_dir=$std{'WORKING'};
 my $data_dir=$std{'DATA'}; my $data_dir=$std{'DATA'};
-my $max_processes = 20; 
  
-# Redirect ouput to logfile (do it only when everything is fine) 
-#my $debug = "$data_dir/.debug_pet_fs_metrics.log"; 
-#open STDOUT, ">$debug" or die "Can't redirect stdout"; 
-#open STDERR, ">&STDOUT" or die "Can't dup stdout"; 
- 
-my $pm = new Parallel::ForkManager($max_processes); 
 our %subjects; our %subjects;
  
-$pm->run_on_finish( 
- sub { my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data) = @_; 
- foreach my $tag (sort keys %{$data}){ 
- $subjects{$ident}{$tag}=${$data}{$tag} 
- } 
- } 
-); 
 my @plist = find(file => 'name' => "*_fbb.nii.gz", '!name' => "*tmp*", in => $w_dir); my @plist = find(file => 'name' => "*_fbb.nii.gz", '!name' => "*tmp*", in => $w_dir);
 my $ofile = $data_dir."/".$study."_fbb_cl.csv"; my $ofile = $data_dir."/".$study."_fbb_cl.csv";
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}{'dg'};  my $dg = $subjects{$subject}{'dg'};
- my @care; + # Apply masks to FBB 
- my %sdata; + my $roi_mask = $roi_paths.$ctx_roi
- $pm->start($subject) and next+ my $order = "fslstats ".$w_dir."/".$dg.$subject."_fbb_mni -k ".$roi_mask.-M"; 
-  + my $ctx = qx/$order/
- # Get FBB image into MNI space + my $roi_mask = $roi_paths.$ctx_roi
- my $order = "fbb2std.sh ".$dg.$subject." ".$w_dir; + $order = "fslstats ".$w_dir."/".$dg.$subject."_fbb_mni -k ".$roi_mask." -M"; 
- print "$order\n"; + print "$order\n"; 
- system($order)+ my $norm = qx/$order/; 
- # Apply masks to FBB + $subjects{$subject}{$npf.'mean'} = $ctx/$norm;
- foreach my $npf (sort keys %ROI_Comps){ +
- my $roi_mask = $roi_paths.$ROI_Comps{$npf}+
- # get mean and std for mask +
- $order = "fslstats ".$dg.$subject."_fbb_mni -k ".$roi_mask." -M -S"; +
- print "$order\n"; +
- (my $mean, my $std) map{/(\d+\.\d*)\s*(\d+\.\d*)/qx/$order/; +
-            $sdata{$npf.'mean'= $mean; +
-            $sdata{$npf.'std'} = $std; +
- }  +
- $pm->finish($subject, \%sdata);  +
- # remove temp dir +
- #remove( \1, ($mdir));+
 } }
-$pm->wait_all_children; 
  
 open OF, ">$ofile"; open OF, ">$ofile";
- +print OF "Subject; SUVRCentilod";
-print OF "Subject"; +
- +
-foreach my  $npf (sort keys %ROI_Comps){ +
-        if($withstd){ +
-                print OF ";$npf","_Mean;","$npf","_STD",";$npf","_c_Mean;","$npf","_c_STD"; +
-        }else{   +
-                print OF ";$npf",";$npf","_c"; +
-        } +
-}+
 print OF "\n"; print OF "\n";
 foreach my $subject (sort keys %subjects){ foreach my $subject (sort keys %subjects){
  print OF "$subject";  print OF "$subject";
- foreach my  $npf (sort keys %ROI_Comps){ + my $mean = $subjects{$subject}{$npf.'mean'}; 
- my $mean = $subjects{$subject}{$npf.'mean'}; + print OF ";$mean"; 
- my $std = $subjects{$subject}{$npf.'std'}; + print OF ";", centiloid_fbb($mean);
- if($withstd){ +
- print OF ";$mean",";$std"; +
- print OF ";", centiloid_fbb($mean), ";", centiloid_fbb($std);  +
- }else{ +
- print OF ";$mean"; +
- print OF ";", centiloid_fbb($mean); +
-+
- }+
  print OF "\n";  print OF "\n";
 } }
 close OF; close OF;
 +</code>
  
-my $zfile = $ofile.'.gz'; +Ahora hago un wrapper que lance todo,
-system("gzip -c $ofile > $zfile");+
  
-if ($attach){ +<code perl fbb_cl_metrics.pl> 
- shit_done basename($ENV{_}), $study, $zfile;+#!/usr/bin/perl 
 +# Copyright 2019 O. Sotolongo <asqwerty@gmail.com> 
 +use strict; use warnings; 
 +  
 +use File::Find::Rule; 
 +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::Remove 'remove'; 
 +use File::Basename qw(basename); 
 + 
 +my $wcut = 0; 
 +my $cfile; 
 + 
 +@ARGV = ("-h") unless @ARGV; 
 + 
 +while (@ARGV and $ARGV[0] =~ /^-/) { 
 +    $_ = shift; 
 +    last if /^--$/; 
 +    if (/^-cut/) { $cfile = shift; chomp($cfile); $wcut = 1;} 
 +    if (/^-h$/) { print_help $ENV{'PIPEDIR'}.'/doc/pet_metrics.hlp'; exit;} 
 +
 +  
 +my $study = shift; 
 +unless ($study) { print_help $ENV{'PIPEDIR'}.'/doc/pet_metrics.hlp'; exit;} 
 +my %std = load_study($study)
 +my $w_dir=$std{'WORKING'}; 
 +my $data_dir=$std{'DATA'}; 
 +my $outdir = "$std{'DATA'}/slurm"; 
 +check_or_make($outdir); 
 +our %subjects; 
 + 
 +my @plist = find(file => 'name' => "*_fbb.nii.gz"'!name' => "*tmp*", in => $w_dir); 
 +my $ofile = $data_dir."/".$study."_fbb_cl.csv"; 
 +my $patt = '([A-Z,a-z]{1,4})(\d{1,6})_fbb'; 
 + 
 +my @pets; 
 + 
 +if ($cfile){ 
 + my %cuts = get_pair($data_dir."/".$cfile); 
 + foreach my $cut (keys %cuts){ 
 + if(grep {/$cut/} @plist){ 
 + $pets[$cut] = $plist[$cut]; 
 +
 + }
 }else{ }else{
- achtung basename($ENV{_}), $ofile, $study;+ @pets = @plist;
 } }
  
 +
 +foreach my $pet (sort @pets){   
 +    (my $dg,my $subject) = $pet =~ /$patt/;
 +    if($subject){
 + $subjects{$subject}{'dg'} = $dg;
 + $subjects{$subject}{'pet'} = $pet;
 + }
 +}
 +
 +foreach my $subject (sort keys %subjects){
 + my $norm;
 + my $dg = $subjects{$subject}{'dg'};
 + my $subj = $dg.$subject;
 + #Making sbatch scripts
 + # Get FBB image into MNI space
 + my $order = $ENV{'PIPEDIR'}."/bin/fbb2std.sh ".$subj." ".$w_dir;
 + my $orderfile = $outdir.'/'.$subj.'_fbb_reg.sh';
 + open ORD, ">$orderfile";
 + print ORD '#!/bin/bash'."\n";
 + print ORD '#SBATCH -J fbb2std_'.$study."\n";
 + print ORD '#SBATCH --time=4:0:0'."\n"; #si no ha terminado en X horas matalo
 + print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
 + print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
 + print ORD '#SBATCH -o '.$outdir.'/fbb2std-'.$subj.'-%j'."\n";
 + print ORD "srun $order\n";
 + close ORD;
 + system("sbatch $orderfile");
 +}
 +my $order = $ENV{'PIPEDIR'}."/bin/fbb_cl_masks.pl ".$study." ".($wcut?"-cut $cfile":"");
 +my $orderfile = $outdir.'/fbb_masks.sh';
 +open ORD, ">$orderfile";
 +print ORD '#!/bin/bash'."\n";
 +print ORD '#SBATCH -J fbb2std_'.$study."\n";
 +print ORD '#SBATCH --time=24:0:0'."\n"; #si no ha terminado en X horas matalo
 +print ORD '#SBATCH --mail-type=FAIL,END'."\n"; #email cuando termine o falle
 +print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
 +print ORD '#SBATCH -o '.$outdir.'/fbbmasks-%j'."\n";
 +print ORD "srun $order\n";
 +close ORD;
 +my $xorder = 'sbatch --dependency=singleton'.' '.$orderfile;
 +exec($xorder);
 +</code>
 +
 +
 +
 +<code bash>
 +[osotolongo@detritus centiloid]$ fbb_cl_metrics.pl centiloid
 +.....
 +[osotolongo@detritus centiloid]$ squeue
 +             JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
 +             15944     devel fbb2std_ osotolon PD       0:00      1 (Dependency)
 +             15920     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15921     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15922     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15923     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15924     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15925     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15926     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15927     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15928     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15929     devel fbb2std_ osotolon  R       0:04      1 brick01
 +             15930     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15931     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15932     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15933     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15934     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15935     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15936     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15937     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15938     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15939     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15940     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15941     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15942     devel fbb2std_ osotolon  R       0:01      1 brick01
 +             15943     devel fbb2std_ osotolon  R       0:01      1 brick01
 +</code>
 +
 +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
 +</code>
 +
 +Compruebo el modelo,
 +
 +<code>
 +> d <- read.csv("toreview.dat", sep = " ", header = TRUE)
 +> m <- lm(d$SUVR ~ d$WC_suvr)
 +> summary(m)
 +
 +Call:
 +lm(formula = d$SUVR ~ d$WC_suvr)
 +
 +Residuals:
 +      Min        1Q    Median        3Q       Max 
 +-0.077423 -0.020244 -0.007479  0.020566  0.074633 
 +
 +Coefficients:
 +            Estimate Std. Error t value Pr(>|t|)    
 +(Intercept)  0.07059    0.02185   3.231  0.00286 ** 
 +d$WC_suvr    0.94779    0.01609  58.921  < 2e-16 ***
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +Residual standard error: 0.03411 on 32 degrees of freedom
 +Multiple R-squared:  0.9909,    Adjusted R-squared:  0.9906 
 +F-statistic:  3472 on 1 and 32 DF,  p-value: < 2.2e-16
 +</code>
 +y el ajuste,
 +<code>
 +gnuplot> f(x) = m*x +n
 +gnuplot> fit f(x) "toreview.dat" u 2:4 via m,n
 +iter      chisq       delta/lim  lambda                          
 +   0 3.4203670136e+01   0.00e+00  1.19e+00    1.000000e+00   1.000000e+00
 +   1 3.2574055507e-01  -1.04e+07  1.19e-01    7.952201e-01   2.956239e-01
 +   2 4.1098565887e-02  -6.93e+05  1.19e-02    1.042966e+00  -5.844827e-02
 +   3 4.1072773521e-02  -6.28e+01  1.19e-03    1.045444e+00  -6.184558e-02
 +   4 4.1072773521e-02  -5.95e-07  1.19e-04    1.045444e+00  -6.184591e-02
 +iter      chisq       delta/lim  lambda                          
 +
 +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    (FIT_NDF)                        : 32
 +rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.0358263
 +variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00128352
 +
 +Final set of parameters            Asymptotic Standard Error
 +=======================            ==========================
 +m               = 1.04544          +/- 0.01774      (1.697%)
 +n               = -0.0618459       +/- 0.02406      (38.9%)
 +
 +correlation matrix of the fit parameters:
 +                m      n      
 +m               1.000 
 +n              -0.967  1.000 
 </code> </code>
 +Con lo que tengo una pendiente de $1.05 \pm 0.02$ y una $R^2 = 0.9906$ que esta bastante bien.
 +{{:neuroimagen:c_fit.png|}}
  
neuroimagen/centiloid.txt · Last modified: 2020/08/04 10:58 by 127.0.0.1