Detritus Wiki

User Tools

Site Tools

Trace:
neuroimagen:centiloid

Differences

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

Link to this comparison view

Next revision		Previous revision
neuroimagen:centiloid [2019/04/09 09:48] – created osotolongoneuroimagen:centiloid [2020/08/04 10:58] (current) – external edit 127.0.0.1
Line 11: Line 11:
 Esto es sencillo de implementar pero antes hay que calibrar el metodo de obtencion de //SUVR// de la pipeline con las imagenes procedentes de [[http://www.gaain.org/centiloid-project|GAAIN]]. Esto es sencillo de implementar pero antes hay que calibrar el metodo de obtencion de //SUVR// de la pipeline con las imagenes procedentes de [[http://www.gaain.org/centiloid-project|GAAIN]].
  
-==== Comparando con GAAIN =====+==== Procesando GAAIN =====
 Basicamente descargamos las imagenes y los valores de centiloid calculados en GAAIN, Basicamente descargamos las imagenes y los valores de centiloid calculados en GAAIN,
  
Line 144: Line 144:
 </code> </code>
  
 +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>
 +#!/bin/sh
 +
 +study=$1
 +shift
 +
 +id=$1
 +shift
 +
 +tdir=$1
 +shift
 +
 +wdir=$1
 +shift
 +
 +items=(`ls ${tdir}/${id}* | grep -v "_" | grep -v ".json"`)
 +#shift
 +
 +debug=0
 +
 +#Now get the uncorrected PETs and register to user space MRI
 +for i in ${!items[*]}; do
 +        tf=`printf "${id}s%04d" $i`
 +        #${FSLDIR}/bin/fslreorient2std ${tdir}/${tf} ${tdir}/${id}_tmp
 +        ${FSLDIR}/bin/imcp ${tdir}/${tf} ${tdir}/${id}_tmp
 +        ${FSLDIR}/bin/flirt -ref ${wdir}/${id}_struc -in ${tdir}/${id}_tmp -omat ${tdir}/${tf}_pet2struc.mat -out ${tdir}/${tf}_reg
 +        #${FSLDIR}/bin/flirt -ref ${wdir}/${id}_brain -in ${tdir}/${id}_tmp -init ${tdir}/${tf}_pet2struc.mat -out ${tdir}/${tf}_reg
 +done
 +if [ ${#items[@]} -gt 1 ]; then
 +echo ${#items[@]}
 +a=`for i in ${!items[*]}; do printf " ${tdir}/${id}s%04d_reg " $i; done`
 +${FSLDIR}/bin/fslmerge -t ${wdir}/${id}_tmp_mvc $a
 +#${FSLDIR}/bin/fslmaths ${dir}/${id}_tmp_pet_in_struc -thr 0 -mas ${dir}/${id}_brain ${dir}/${id}_pet_in_struc
 +#${FSLDIR}/bin/fslmaths ${dir}/${id}_tmp_pet_in_struc -mas ${dir}/${id}_brain ${dir}/${id}_pet_in_struc
 +
 +${FSLDIR}/bin/mcflirt -in ${wdir}/${id}_tmp_mvc -out ${wdir}/${id}_tmp_corr
 +${PIPEDIR}/bin/4dmean.pl ${wdir}/${id}_tmp_corr
 +${FSLDIR}/bin/flirt -ref ${wdir}/${id}_struc -in ${wdir}/${id}_mean -omat ${wdir}/${id}_fbb2struc.mat -out ${wdir}/${id}_fbb
 +else
 +tf=`printf "${id}s%04d" ${item[0]}`
 +${FSLDIR}/bin/mcflirt -in ${tdir}/${tf}_reg -out ${wdir}/${id}_tmp_corr
 +${FSLDIR}/bin/imcp ${wdir}/${id}_tmp_corr ${wdir}/${id}_fbb
 +fi
 +
 +if [ $debug = 0 ] ; then
 +    rm ${tdir}/${id}_tmp*
 +    rm ${wdir}/${id}_tmp*
 +fi
 +</code>
 +De aqui se puede hacer lo usual,
 +<code bash>
 +[osotolongo@detritus centiloid]$ fbb_correct.pl centiloid
 +</code>
 +Y podemos revisar el report,
 +{{:neuroimagen:fbb_report_crop.png| registration report}}
 +Todo parece ir bien asi que,
 +<code bash>
 +[osotolongo@detritus ~]$ parallel_fbb_rois_metrics.pl centiloid
 +</code>
 +
 +Ahora solo hay que compara los valores globales de SUVR con los de la tabla de GAAIN. 
 +
 +<code bash>
 +[osotolongo@detritus centiloid]$ awk -F ";" '{print $1,$2,$4}' FBB_suvr_centiloid.csv | grep -v Y | sort -n > reference.dat
 +[osotolongo@detritus centiloid]$ join calcs.dat reference.dat > toreview.dat
 +</code>
 +
 +A ver que tal ese modelo,
 +
 +<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
 +iter      chisq       delta/lim  lambda   m             n            
 +   0 3.6936711672e+01   0.00e+00  1.41e+00    1.000000e+00   1.000000e+00
 +   1 3.6992605848e-01  -9.88e+06  1.41e-01    6.105751e-01   4.583953e-01
 +   2 7.9414981509e-02  -3.66e+05  1.41e-02    8.423745e-01   2.635659e-02
 +   3 7.9329087727e-02  -1.08e+02  1.41e-03    8.467400e-01   1.870027e-02
 +   4 7.9329087724e-02  -3.50e-06  1.41e-04    8.467408e-01   1.869890e-02
 +iter      chisq       delta/lim  lambda   m             n            
 +
 +After 4 iterations the fit converged.
 +final sum of squares of residuals : 0.0793291
 +rel. change during last iteration : -3.5012e-11
 +
 +degrees of freedom    (FIT_NDF)                        : 22
 +rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.0600489
 +variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00360587
 +
 +Final set of parameters            Asymptotic Standard Error
 +=======================            ==========================
 +m               = 0.846741         +/- 0.02881      (3.403%)
 +n               = 0.0186989        +/- 0.0497       (265.8%)
 +
 +correlation matrix of the fit parameters:
 +                m      n      
 +m               1.000 
 +n              -0.969  1.000 
 +</code>
 +
 +{{:neuroimagen:fit.png|}}
 +
 +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:
 +
 +               m      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 =====
 +
 +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>
 +
 +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   m             n            
 +   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   m             n            
 +
 +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.1554803328.txt.gz · Last modified: 2020/08/04 10:45 (external edit)