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 revision Previous revision
Next revision
Previous revision
neuroimagen:centiloid [2019/04/16 09:16]
osotolongo [Implementando pipeline Centiloid]
neuroimagen:centiloid [2020/08/04 10:58] (current)
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 275: Line 275:
 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. 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 =====
  
Line 603: Line 667:
 </code> </code>
  
-**Problemas:** Ekl procedimiento de registro demora bastante asi que estaria bien paralelizarlo en el cluster y hacer las mascaras luego. Para ello hay que separar los procesos en distintos scripts, pero no deberia ser dificil. 
  
-===== Usando el cluster ===== 
- 
-<code perl fbb_cl_metrics.pl> 
-#!/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 $withstd = 0; 
-my $cfile; 
- 
-@ARGV = ("-h") unless @ARGV; 
- 
-while (@ARGV and $ARGV[0] =~ /^-/) { 
-    $_ = shift; 
-    last if /^--$/; 
-    if (/^-std$/) {$withstd = 1;} 
-    if (/^-cut/) { $cfile = shift; chomp($cfile);} 
-    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{ 
- @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 
- ##### me quede aqui ---> quitaresto y hacer como en el proc de FS 
- 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." ".($withstd?"-std":""); 
-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 perl fbb_cl_masks.pl> 
-#!/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 %ROI_Comps = ( 
- "Cgray" => "voi_CerebGry_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/2mm/'; 
-my $withstd = 0; 
-my $cfile; 
- 
-@ARGV = ("-h") unless @ARGV; 
- 
-while (@ARGV and $ARGV[0] =~ /^-/) { 
-    $_ = shift; 
-    last if /^--$/; 
-    if (/^-std$/) {$withstd = 1;} 
-    if (/^-cut/) { $cfile = shift; chomp($cfile);} 
-} 
-  
-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'}; 
- 
-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{ 
- @pets = @plist; 
-} 
- 
-foreach my $pet (sort @pets){    
-    (my $dg,my $subject) = $pet =~ /$patt/; 
-    if($subject){ 
- $subjects{$subject}{'dg'} = $dg; 
- $subjects{$subject}{'pet'} = $pet; 
- } 
-} 
-my %sdata; 
-foreach my $subject (sort keys %subjects){ 
- my $care; 
- my $norm; 
- my $dg = $subjects{$subject}{'dg'}; 
- # Apply masks to FBB 
- foreach my $npf (sort keys %ROI_Comps){ 
- my $roi_mask = $roi_paths.$ROI_Comps{$npf}; 
- # get mean and std for mask 
- my $order = "fslstats ".$w_dir."/".$dg.$subject."_fbb_mni -k ".$roi_mask." -M -S"; 
- print "$order\n"; 
- (my $mean, my $std) = map{/(\d+\.\d*)\s*(\d+\.\d*)/} qx/$order/; 
-            $subjects{$subject}{$npf.'mean'} = $mean; 
-            $subjects{$subject}{$npf.'std'} = $std; 
- }  
-} 
- 
-open OF, ">$ofile"; 
-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"; 
-foreach my $subject (sort keys %subjects){ 
- print OF "$subject"; 
- foreach my  $npf (sort keys %ROI_Comps){ 
- my $mean = $subjects{$subject}{$npf.'mean'}; 
- my $std = $subjects{$subject}{$npf.'std'}; 
- 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"; 
-} 
-close OF; 
-</code> 
  
 <code bash> <code bash>
Line 841: Line 700:
              15943     devel fbb2std_ osotolon  R       0:01      1 brick01              15943     devel fbb2std_ osotolon  R       0:01      1 brick01
 </code> </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>
 +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.1555406209.txt.gz · Last modified: 2020/08/04 10:46 (external edit)