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
Next revision Both sides next revision
neuroimagen:centiloid [2019/04/11 08:46]
osotolongo [Procesando GAAIN]
neuroimagen:centiloid [2019/04/16 09:43]
osotolongo [Procesando GAAIN]
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.
  
 + 
 ==== 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.]]+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 ** 
 + 
 +Primero hago un script sencillo que registre con  ANTS al espacio MNI, 
 + 
 +<code bash fbb2std.sh> 
 +#!/bin/sh 
 + 
 +id=$1 
 +shift 
 + 
 +wdir=$1 
 +shift 
 + 
 +echo "I need the FBB image at MNI space"                                                                                                                                                
 +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_2mm.nii.gz ${wdir}/${id}_fbb_t1_mniWarp.nii.gz ${wdir}/${id}_fbb_t1_mniAffine.txt 
 + 
 +</code> 
 + 
 +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 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 $ref_roi =  "voi_WhlCbl_2mm.nii"; 
 +my $ctx_roi = "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 (/^-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 
 + my $roi_mask = $roi_paths.$ctx_roi; 
 + my $order = "fslstats ".$w_dir."/".$dg.$subject."_fbb_mni -k ".$roi_mask." -M"; 
 + my $ctx = qx/$order/; 
 + my $roi_mask = $roi_paths.$ctx_roi; 
 + $order = "fslstats ".$w_dir."/".$dg.$subject."_fbb_mni -k ".$roi_mask." -M"; 
 + print "$order\n"; 
 + my $norm = qx/$order/; 
 + $subjects{$subject}{$npf.'mean'} = $ctx/$norm;  
 +
 + 
 +open OF, ">$ofile"; 
 +print OF "Subject; SUVR; Centilod"; 
 +print OF "\n"; 
 +foreach my $subject (sort keys %subjects){ 
 + print OF "$subject"; 
 + my $mean = $subjects{$subject}{$npf.'mean'}; 
 + print OF ";$mean"; 
 + print OF ";", centiloid_fbb($mean); 
 + print OF "\n"; 
 +
 +close OF; 
 +</code> 
 + 
 +Ahora hago un wrapper que lance todo, 
 + 
 +<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 $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{ 
 + @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>
neuroimagen/centiloid.txt · Last modified: 2020/08/04 10:58 (external edit)