User Tools

Site Tools


neuroimagen:tau_analysis_refact

Analisis de PET-tau

http://detritus.fundacioace.com/files/ADNI_UCBERKELEY_AV1451_Methods_2021-01-14.pdf

Voy a intentar entender el analisis de PET-tau a taves de un set de sujetos bajados de ADNI. Los sujetos deben tener una MRI y un PET-Tau hechos en fecha similar (+/- 6 meses?). En el documento de ADNI esta mas o menos explicado el analisis y con un poco de maña y otro de google podemos armar lo que falta.

Procesando

Basicamente se han de hacer algunas tareas que conducen a dos metricas distintas.

  • Registro de PET-Tau a espacio nativo MRI
  • Aplicar mascaras provenientes de la segmentacion FS a imagen PET-Tau
  • Calcular imagen SUVR segun valor en ROR
  • Calcular metricas sin corregir
  • Hacer PVC a imagen registrada
  • Calcular SUVR con correccion

Registro a T1w

El primer paso es colocar la imagen tau (4D de 6 slices) en espacio nativo T1. Debemos haber procesado el T1 con Freesurfer, asi que traemos primero el struc de FS. Ahora hacemos un split de la imagen tau y registramos cada imagen resultante al T1. La imagen tau final sera el valor medio de cada una de estas imagenes registradas.

Utilizamos las 6 imágenes de 5min para la correccion de movimiento

Y para ir preparando el wrapper

el call, usando el modulo SLURM.pm se hace como,

Informe de registro

Tras efectuar el registro se puede generar el informe para ver como ha ido. Aqui el procedimiento es muy parecido al de FBB asi que no vale la pena comentar nada

Hacer mascaras

Ahora debemos seleccionar las regiones de interes donde queremos medir el SUVR. Hay varias opciones, como las regiones de Braak donde normalmente se observa la progresion del AD, un combinacion de estas regiones o una meta-region que definamos.

Las definiciones de cada region se han tomado de los metodos de ADNI para el AV1451.

Ahora tomamos la segmentacion de FS (aparc+aseg), la llevamos a espacio nativo y seleccionamos, segun los LUT, los elementos de cada ROI. Se hace una mascara con cada ROI y simplemente se suman para tener la mascara de la ROI deseada.

Es bastante rapido de hacer

Esto en el wrapper ha de hacerse para cada roi que deseemos

Normalizacion

Mascara para normalizar

Vamos a normalizar por la sustancia gris del cerebelo. Para garantizar que no haya contaminacion en las mediciones de la materia gris del cerebelo tomamos solo la zona inferior (inferior cerebellum gray matter). Para ello se toma el atlas del cerebelo de SUIT en espacio MNI que proporciona AFNI. Se calcula el registro del T1w a espacio MNI y con la inversa se lleva este atlas a espacio nativo. De aqui se seleccionan las regiones inferiores del cerebelo y se hace una mascara del cerebelo calculado por FS. Ahora bien, dado que el cerebelo de SUIT no encaja exactamente con el segmentado de FS, tomamos la parte inferior y la parte superior (de SUIT) y a cada una se le hace un smooth con un kernel gaussiano de 8mm y se comparan. Si la inferior es mayor que la superior y pertenece a la segmentacion de FS, entonces se toma este punto como valido. La mascara resultante es la que se utilizara para normalizar. Esto parece muy complicado pero son dos operaciones matematicas, $B(I-S)$, donde \[ B(x) = \{ \begin{array}{cc} 1, & x>0 \\
0, & \textrm{otherwise} \end{array} \].

Un poco mas complicado pero similar a la extraccion de ROIs

Y en el wrapper se lanza como un proceso independiente en SLURM

Merge Masks

Ahora construimos un 4D con las mascaras que hemos escogido + la que utilizaremos para normalizar. Esto es un solo comando que usa el string que hemos venido construyendo y se lanza directamente tras terminar todas las mascaras

A cada mascara le queda una capa asignada dentro del mismo 4D,

braak 1 braak 2
braak 3 braak 4
braak 5 braak 6

Calculo de valores e imagen SUVR

Ahora, ya podemos aplicar las mascaras que tenemos y calcular cuanto vale el //uptake// en cada una. Ademas como voy a tener el valor de normalizacio, calculo la imagen SUVR

Esto se lanza justo cuando se terminan de hacer las mascaras

PVC correction

Usando la imagen registrada y las mascaras 4D que hemos construido, calcularemos la correccion de volumen parcial (PVC) correspondiente a cada mascara. Vamos a usar PETPVC, desarrollada en University College London para corregir usando el metodo GTM.

Esto es una sola orden simple

y los lanzamos desde el wrapper cuando tengamos el SUVR

PVC correction (metodos)

Por defecto hemos elegido el GMT como metodo de correccion pero es posible elegir otro metodo de manera relativamente sencilla.

Esta, por ejemplo, es la implementación del Multi-target correction (MTC),

que puede añadirse al wrapper con solo algunas lineas

Pero que tambien habria que añadir como opcion en el archivo de las metricas (Ver mas abajo)

Metricas

Cuando tenemos el output del PVC (y los SUVR sin corregir), hemos de escribir toda la informacion en una tabla que agrupe las metricas y los sujetos. Ha de ser una tabla para los SUVR corregidos y la misma tablas para los SUVR sin corregir.

Pero como tenemos los resultados en el mismo formato, solo hay que escribir un algoritmo

Esto ha de lanzarse tras hacer los PVC

Mas sobre metodos de PVC

Para extraer las metricas siguiendo un metodo extra de correccion habra que cambiar la linea,

my @subs = ("pvc", "unc");

y añadir la cabecera del metodo. Por ejemplo, para el MTC,

my @subs = ("pvc", "unc","mtc");

Wrapper

Al final nos ha quedado un algoritmo extremadamente paralelizado y mucho mas rapido (aunque mas complejo) que, por ejemplo el procesamiento PET-FBB. Se hacen muchisimos procesos cortos. Dependiendo de cuantas ROI se procesen puede llegar a dividirse un sujeto en mas de 10 tareas relativamente rapidas. Por ejemplo, el procesamiento por defecto de 57 sujetos se divide en 675 tareas independientes con su arbol de dependencias bien definido.

Todo el procesamiento Tau, del NIfTI a las metricas finales

Testing with ADNI

Procesamiento

Ahora tengo que probar que todo funciona correctamente y comprar los resultados contra los reportados en ADNI. el primer paso debe ser escoger los sujetos que voy a procesar.

Lo que hago es ir a R, instalar el paquete ADNIMERGE y encontrar la correspondencia entre los RID y los PTID ya que los resultados del AV1451 se reportan con el RID pero las imagenes se identifican por el PTID.

Cuando tenga las lista de los PTID debo ir a la pagina de ADNI y buscar, en es ta lista, aquellos que tienen MRI T1 y PET-Tau al mismo tiempo. Entonces simplemente escogo algunos sujetos y bajo las dos imagenes (MRI y PET) que matcheen bien (uy que palabreja). Los tags son muy distintos aqui asi que debo revisar bien el proceso de conversion a BIDS. Pero con unos 50 MRI que se convierten automaticamente tengo para empezar. Despues tocara revisar que los PETs esten OK.

Entonces sigo el procedimiento usual de cualquier proyecto, hago la DB, convierto a BIDS, paso el recon a todas las MRI. Primero vamos a construir el proyecto de la manera usual y a ejecutar la segmentacion de Freesurfer y sacamos al final los datos que necesitemos para comparar.

Data ADNI

A ver que hay dentro,

Ah, voy a sacar una medida de aqui a ver,

y a preprocesar un poco

Comparando con ADNI

Pegando en bash

Esto lo puedo mirar con R

Después de haber revisado todos los procedimientos, comparado de varias maneras con los datos de ADNI y usado varias maneras de normalización hemos de tomar decisiones.

  1. La normalización se hará por inferior cerebellum gray matter (ICGM) que parece ser la mas estable.

Aqui las comparaciones,

Uncorrected Braak 1 PVCed Braak 1
Uncorrected Braak 3+4 PVCed Braak 3+4
Uncorrected Braak 5+6 PVCed Braak 5+6
Uncorrected Meta Temporal PVCed Meta Temporal

$R^2$

Uncorrected Corrected
Braak 1 0.97 0.96
Braak 3+4 0.99 0.98
Braak 5+6 0.99 0.95
Meta temporal 0.99 0.99

$slope$

Uncorrected Corrected
Braak 1 1.11 0.90
Braak 3+4 1.05 0.98
Braak 5+6 1.02 0.92
Meta temporal 1.07 0.98

Comparando con otro metodo de correccion (MTC)

Comparando Geometric Matrix Transfer (GMT) com Multi-target Correction (MTC)

Braak 1 Meta-temporal

y el ajuste,

$R^2$ $slope$
Braak 1 0.99 1.02
Meta-temporal 1.0 0.98

How to cite

PET-Tau scans have been processing following ADNI method, freely available at ADNI website for registered researchers [1, 2]. First, the 4D 30min PET image was separated into 6 5min images for movement correction and registered into the T1w subject native space. The regions of interest for calculating SUVR, as well as for normalization, where built using the Freesurfer 7.2 segmentation of the T1w. Partial volume correction was calculated using GTM method [3] with PETPVC [4]. Corrected and uncorrected SUVR where calculated using inferior cerebellum gray matter as region of reference.

[1] http://loni.adni.usc.edu. Flortaucipir (AV-1451) processing methods. (https://adni.bitbucket.io/reference/docs/UCBERKELEYAV1451/UCBERKELEY_AV1451_Methods_2021-01-14.pdf)

[2] Maass, A., et al., Comparison of multiple tau-PET measures as biomarkers in aging and Alzheimer's disease. Neuroimage, 2017. 157: p. 448-463.

[3] Rousset, O.G., Y. Ma, and A.C. Evans, Correction for partial volume effects in PET: principle and validation. J Nucl Med, 1998. 39(5): p. 904-11.

[4] BA Thomas, V Cuplov, A Bousse, A Mendes, K Thielemans, BF Hutton, K Erlandsson. PETPVC: a toolbox for performing partial volume correction techniques in positron emission tomography. Physics in Medicine and Biology 61 (22), 7975.

Plotting with ggseg

Para mostrar los valores de $\tau$-SUVR voy a usar ggseg. Esto es un script python que voy a generar a partir del output del analisis.

Para empezar, necesito saber: proyecto, radiotracer, sujeto y modalidad que voy a graficar. Asi que estos valores los voy a meter como inputs de un script Perl y a partir d ahi decido que hacer,

#!/usr/bin/perl
use strict;
use warnings;
use NEURO4 qw(check_pet check_subj load_project);
use Data::Dump qw(dump);
my $tracer;
my $subject;
my $mod = 'unc';
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
	$_ = shift;
	if (/^-tracer/) {$tracer = shift; chomp($tracer);}
	if (/^-id/) {$subject = shift; chomp($subject);}
	if (/^-mod/) {$mod = shift; chomp($mod);}
}
my $study = shift;
die "Should supply project name\n" unless $study;
die "Should supply subject ID\n" unless $subject;

Ahora cargo las propiedades del proyecto y leo los datos en un hash,

my %std = load_project($study);
my $ifile = $std{'DATA'}.'/'.$study.'_tau_suvr_'.$tracer.'_'.$mod.'.csv';
open IDF, "<$ifile";
my $seeker;
my %taud;
my @levels;
while(<IDF>){
	unless ($seeker){
		@levels = split ", ", $_;
		chomp @levels;
	}else{
		my @values = split ", ", $_;
		chomp @values;
		if ($values[0] eq $subject){
			for (my $i = 1; $i < scalar @values; $i++){
				$taud{$levels[$i]} = $values[$i] unless $levels[$i] eq 'extra';
			}
		}
	}
	$seeker++;
}
close IDF;

Ahora, segun estos datos toca escribir el script python. En principio no sabemos que ROIs han sido definidas para calcular el SUVR pero si que estan definidas dentro del pipeline. Asi que lo que hago es aislar cada una y traducirlas a las ROI de FS.

my $odata = "data = {\n";
foreach my $tag (sort keys %taud){
	my $libfile = $ENV{'PIPEDIR'}.'/lib/tau/'.$tag.'.roi';
	open ILF, "<$libfile";
	while (<ILF>){
		my ($hand, $roi) = /\d+,(\w)_(\w+)/;
                unless ($hand eq 'R'){
		      $odata.="'".$roi."_left': ".$taud{$tag}.",\n";
		      $odata.="'".$roi."_right': ".$taud{$tag}.",\n";
                }
	}
	close ILF;
}
$odata.="}\n";

Eso lo vuelco dentro del script python,

my $pfile = $std{'DATA'}.'/'.$study.'_'.$tracer.'_'.$subject.'_'.$mod.'.py';
open OPF, ">$pfile";
print OPF '#!/usr/bin/python3'."\n";
print OPF "import ggseg\n";
print OPF "$odata\n";
print OPF "ggseg.plot_dk(data, cmap='Spectral', figsize=(15,15),";
print OPF "background='k', edgecolor='w', bordercolor='gray',";
print OPF "ylabel='PET-".$tracer." SUVR', title='".$subject."')\n";
close OPF;

El call es algo como,

./generate_tau_pyplot.pl -tracer ro948 -id 0075 -mod unc f5cehbi

que nos deja un script python como,

#!/usr/bin/python3
import ggseg
data = {
'entorhinal_left': 0.880505023804403,
'entorhinal_right': 0.880505023804403,
....
'paracentral_left': 1.31710059291554,
'paracentral_right': 1.31710059291554,
}
 
ggseg.plot_dk(data, cmap='Spectral', figsize=(15,15),background='k', edgecolor='w', bordercolor='gray',ylabel='PET-ro948 SUVR', title='0075')

y cuando se ejecuta podemos salvar una imagen con los distintos valores de $\tau$-SUVR por las ROIs de FS,

PET-Tau FACEHBI project

La particularidad de este proyecto es que hay dos marcadores PI2620 y RO948. En principio esta diferencia se indica en los tags. Los json deberían indicar el marcador,

         "SeriesDescription": "[DY_NAC] Dynamic Brain",
        "ProtocolName": "Head/FACEHBI_PI2620_",

Así que añadimos un par de reglas al archivo de conversion,

                {
                        dataType: pet,
                        modalityLabel: tau,
                        customLabels: pi2620,
                        criteria: {
                                SeriesDescription: *DY_NAC*,
                                ProtocolName : *PI2620*
                        }
                },
                {
                        dataType: pet,
                        modalityLabel: tau,
                        customLabels: ro948,
                        criteria: {
                                SeriesDescription: *DY_NAC*,
                                ProtocolName : *RO948*
                        }
                }

Ahora, la configuracion del proyecto es la siguiente,

[osotolongo@brick03 f5cehbi]$ cat ~/.config/neuro/f5cehbi.cfg
DATA = /old_nas/f5cehbi
SRC = /nas/data/f5cehbi/raw
PET = /nas/data/f5cehbi/raw_tau/PI2620
WORKING = /old_nas/f5cehbi/working
BIDS = /old_nas/f5cehbi/bids

y dependiendo del marcador que vayamos a procesar cambiamos el path de PET al correcto,

  1. /nas/data/f5cehbi/raw_tau/PI2620
  2. /nas/data/f5cehbi/raw_tau/RO948

y ejecutamos,

$ pet2bids.pl f5cehbi

Pero, hay alguna imagen mal identificada, y para encontrarla hacemos,

[osotolongo@brick03 f5cehbi]$ ls /nas/data/f5cehbi/raw_tau/PI2620 > upload_pi2620.list
[osotolongo@brick03 f5cehbi]$ find bids/ -name "*pi2620_*tau.json" | awk -F"-" '{print $2}' | awk -F"/" '{print $1}' > onbids.list
[osotolongo@brick03 f5cehbi]$ grep -f upload_pi2620.list f5cehbi_pet.csv | awk -F";" '{print $1}' > subjects_uploaded.list
[osotolongo@brick03 f5cehbi]$ grep -v "`cat onbids.list`" subjects_uploaded.list
0071
0078

y ahora estos los he de copiar manualmente,

[osotolongo@brick03 f5cehbi]$ grep "DY_NAC" bids/tmp_dcm2bids/sub-0071/*.json
bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.json:	"SeriesDescription": "[DY_NAC] Dynamic Brain",
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.json bids/sub-0071/pet/sub-0071_pi2620_tau.json
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0071/237420_FPT002_Head_FACEHBI_TAU_30M_20220117132341.nii.gz bids/sub-0071/pet/sub-0071_pi2620_tau.nii.gz
[osotolongo@brick03 f5cehbi]$ grep "DY_NAC" bids/tmp_dcm2bids/sub-0078/*.json
bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.json:	"SeriesDescription": "[DY_NAC] Dynamic Brain",
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.json bids/sub-0078/pet/sub-0078_pi2620_tau.json
[osotolongo@brick03 f5cehbi]$ cp bids/tmp_dcm2bids/sub-0078/260730_FPT001_Head_FACEHBI_TAU_30M_20220117140057.nii.gz bids/sub-0078/pet/sub-0078_pi2620_tau.nii.gz

Y ya estamos. Ahora a procesar, pero antes hay que hacer limpieza,

[osotolongo@brick03 f5cehbi]$ rm -rf working/taus working/.tmp_* working/*_tau.nii.gz working/*_suvr* working/*_masks* working/*_pvc* working/*_mtc* working/*_unc*

Ahora si,

[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer pi2620 f5cehbi

y al terminar tendremos los resultados en los archivos,

[osotolongo@brick03 f5cehbi]$ ls f5cehbi_tau_suvr_pi2620_*
f5cehbi_tau_suvr_pi2620_mtc.csv  f5cehbi_tau_suvr_pi2620_pvc.csv  f5cehbi_tau_suvr_pi2620_unc.csv
[osotolongo@brick03 f5cehbi]$ head f5cehbi_tau_suvr_pi2620_unc.csv
Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6, extra
0071, 0.642265450913569, 0.445293067511959, 0.579153982059088, 0.747920179741428, 0.72149021446756, 0.610043587665485, 0.456761589685176
0075, 0.712122676066806, 0.38863874053704, 0.572840374589894, 0.800212087107015, 0.729010277863309, 0.662967042884597, 0.394067798859667
0076, 0.83589891768702, 0.763072917986403, 0.98202394592715, 1.1291513524068, 1.07060320066878, 0.942356703562845, 0.809746693748801
0078, 0.862191018953226, 0.627670886199243, 0.862667901143999, 1.10188426829118, 1.04589664481027, 0.884161149229649, 0.63109502246773
0087, 0.830211995236038, 0.456955305688661, 0.528028370642198, 1.02278147223258, 0.706801020407484, 0.549391337405959, 0.484003816604253
0092, 0.745907131688215, 0.688052328589026, 0.991095394357976, 1.35968730905804, 1.30310065337532, 1.1900620384451, 0.807880424703717
0098, 0.702882998845167, 0.711627059244375, 0.856521933131609, 0.990771318090343, 0.936773228442926, 0.842062696773501, 0.764783413964935
0099, 0.823298330211384, 0.705876423260448, 0.875380807683978, 1.1361288650671, 0.961009113607904, 0.86318920811104, 0.92071927514608
0100, 0.857882757785999, 0.656176690905404, 1.12023931254886, 1.17050947142544, 1.30016270374768, 1.19007631440357, 0.857547667880239

Y ahora voy a hacer el otra marcador. Edito el config

[osotolongo@brick03 f5cehbi]$ cat ~/.config/neuro/f5cehbi.cfg
DATA = /old_nas/f5cehbi
SRC = /nas/data/f5cehbi/raw
PET = /nas/data/f5cehbi/raw_tau/RO948
WORKING = /old_nas/f5cehbi/working
BIDS = /old_nas/f5cehbi/bids

Hago los bids,

[osotolongo@brick03 f5cehbi]$ pet2bids.pl f5cehbi

Limpio el directorio de trabajo,

[osotolongo@brick03 f5cehbi]$ cat clean_tau.sh
#!/bin/bash
rm -rf working/taus working/.tmp_* working/*_tau.nii.gz working/*_suvr* working/*_masks* working/*_pvc* working/*_mtc* working/*_unc*
[osotolongo@brick03 f5cehbi]$ ./clean_tau.sh

Lanzo los procesos,

[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer ro948 f5cehbi

y espero los resultados,

[osotolongo@brick03 f5cehbi]$ head f5cehbi_tau_suvr_ro948_pvc.csv
Subject, braak_1, braak_2, braak_3, braak_4, braak_5, braak_6, extra
0075, 1.04088050314465, 0.64527944539737, 1.15851915380217, 1.52799456832476, 2.21913236134934, 2.39967124070898, 1.30373785020011
0115, 0.908310815175864, 0.492099523944246, 1.10060882563261, 1.43477011430239, 1.83175427770736, 1.77909186294406, 0.796637090468415
0120, 1.46854534539977, 1.07698127739582, 1.8255003454564, 1.66630044512906, 2.06223198812989, 2.256736399778, 1.52880200251447

Integrando el proyecto en XNAT

Subir un PET a XNAT es trivial. El problema aqui es que hay dos PETs y los DICOM no contienen suficiente información m( para que XNAT distinga ambos. Asi que tiende a agruparlos en el mismo PET.

Lo que voy a hacer es fijar el experiment_id y asignar a cada uno el PET que le corresponde. Ejemplo,

$ xnatapic upload_dicom --project_id tau0 --subject_id ${sbj} --experiment_id ${sbj}_PI2620 /nas/data/f5cehbi/raw_tau/PI2620/${pet}

O en contexto,

$ while read -r line; do pet=$(echo ${line} | awk -F"," '{print $3}'); sbj=$(echo ${line} | awk -F"," '{print $1}'); if [ -d /nas/data/f5cehbi/raw_tau/PI2620/${pet} ]; then xnatapic upload_dicom --project_id tau0 --subject_id ${sbj} --experiment_id ${sbj}_PI2620 /nas/data/f5cehbi/raw_tau/PI2620/${pet}; else echo "No ${pet}"; fi; done < fpt_matches.csv

Lo que quedaría como,

Pero como llego a tener el archivo necesario para esto (fpt_matches.csv)?

Lo primero es obtener la lista de los PET del archivo compartido y arreglarle los endline,

$ tr -d '\r' < pet_tau_list.csv > pet_tau_list_proper.csv

Ahora tomamos los codigos del proyecto (f5cehbi),

$ sed 's/;/,/' f5cehbi_mri.csv > f5cehbi_codes.csv
$ sort -t, -k 2 f5cehbi_codes.csv > f5cehbi_codes_sorted.csv

Y todo esto lo juntamos,

$ tail -n +2 pet_tau_list_proper.csv | awk -F"," '{print $1","$3}' | sort -t, > fpt_list.csv
$ join -t, -1 2 -2 1 f5cehbi_codes_sorted.csv fpt_list.csv > fpt_matches.csv

Por supuesto que este archivo es necesario solo una vez, para ambos tracers. Pero debe actualizarse cuando tenemos nuevos sujetos.

Extracción y presentación de datos

Ya hemos calculado los SVUR pero ahora necesitamos ponerlo de manera comprensible en un archivo para entregar. Primero sacamos las fechas, aprovechando el formato de almacenamiento en XNAT,

$ xnatapic list_experiments --project_id tau0 --modality PET --label --date > xnat_tau0_experiments.csv
$ grep RO948 xnat_tau0_experiments.csv | awk -F"," '{print $3","$2}' | sed 's/_.*,/,/' | sed '1iPSubject,Date' > ro948_dates.csv
$ grep PI2620 xnat_tau0_experiments.csv | awk -F"," '{print $3","$2}' | sed 's/_.*,/,/' | sed '1iPSubject,Date' > pi2620_dates.csv

Ahora hay que unirlo con los resultados pero necesitamos unificar los codigos,

$ sed 's/;/,/' f5cehbi_mri.csv | sort -t, -k 1 | sed '1iSubject,PSubject' > codes.csv
$ join -t, codes.csv f5cehbi_tau_suvr_pi2620_unc.csv | sed 's/ //g' > pi2620_unc_wcodes.csv
$ (head -n 1 pi2620_unc_wcodes.csv && tail -n +2 pi2620_unc_wcodes.csv | sort -t, -k 2) > pi2620_unc_wcodes_sorted.csv
$ join -t, -1 1 -2 2 pi2620_dates.csv pi2620_unc_wcodes_sorted.csv > tau_tmp/PI2620_UNC.csv

La ultima parte hay que hacerla para cada uno de los sheets que se quieran exportar,

$ ls -l tau_tmp/
total 28
-rw-r--r-- 1 osotolongo osotolongo  625 Apr 19 10:48 INFO.csv
-rw-rw-r-- 1 osotolongo osotolongo 4758 Apr 19 10:49 PI2620_PVC.csv
-rw-rw-r-- 1 osotolongo osotolongo 4870 Apr 19 10:47 PI2620_UNC.csv
-rw-rw-r-- 1 osotolongo osotolongo 1055 Apr 19 10:51 RO948_PVC.csv
-rw-rw-r-- 1 osotolongo osotolongo 1071 Apr 19 10:52 RO948_UNC.csv

y ahora voy a hacer la conversion a xls,

dir2xls.pl
#!/usr/bin/perl
#
# To convert a bunch of CSV files into a single XLS
#
# Copyleft 2022 O. Sotolongo <asqwerty@gmail.com>
#
use strict;
use warnings;
use Text::CSV qw( csv );
use Spreadsheet::Write;
 
my $idir;
my $ofile;
@ARGV = ("-h") unless @ARGV;
while (@ARGV and $ARGV[0] =~ /^-/) {
        $_ = shift;
        last if /^--$/;
        if (/^-i/) { $idir = shift; chomp($idir);}
        if (/^-o/) { $ofile = shift; chomp($ofile);}
}
die "Should supply results directory\n" unless $idir;
die "Should output filename\n" unless $ofile;
 
$ofile =~ s/\.(\w*)?$/.xls/;
opendir (DIR, $idir);
my @ifiles = grep(/\.csv/, readdir(DIR));
close DIR;
 
my $workbook = Spreadsheet::Write->new(file => $ofile, sheet => 'INFO');
my $infof = $idir.'/INFO.csv';
if( -e $infof ){
        my $inf_data =  csv (in => $infof);
        for my $i (0 .. $#{$inf_data}) {
                my $row = $inf_data->[$i];
                $workbook->addrow($row);
        }
}
foreach my $ifile (@ifiles){
        unless( $ifile eq 'INFO.csv'){
                my $ipath = $idir.'/'.$ifile;
                my $idata = csv (in => $ipath); # as array of array
                (my $shname = $ifile) =~ s/\.csv$//;
                $workbook->addsheet($shname);
                for my $i (0 .. $#{$idata}) {
                        my $row = $idata->[$i];
                        $workbook->addrow($row);
                }
        }
}
$workbook->close();

Voy a dejar una copia de esto en github. :TO-DO:

$ ./dir2xls.pl -i tau_tmp/ -o tau_results_20220419.xls

Thinking around

$ sed 's/ //g' f5cehbi_tau_suvr_ro948_unc.csv | awk '{print $0",ro948"}' | sed 's/extra,ro948/extra,tracer/' > f5cehbi_tau_suvr_ro948_unc_tracer.csv
$ sed 's/ //g' f5cehbi_tau_suvr_pi2620_unc.csv | awk '{print $0",pi2620"}' | sed 's/extra,pi2620/extra,tracer/' > f5cehbi_tau_suvr_pi2620_unc_tracer.csv
$ tail -n +2 f5cehbi_tau_suvr_pi2620_unc_tracer.csv > tmp.csv
$ cat f5cehbi_tau_suvr_ro948_unc_tracer.csv tmp.csv > tau_suvr_unc.csv

Y esto voy a compararlo, a ver como se distinguen, en media, los uptakes de los tracers,

comp_tracers.r
library("ggplot2")
input_file="tau_suvr_unc.csv"
unc <- read.csv(input_file)
rois = c("braak_1","braak_2","braak_3","braak_4","braak_5","braak_6")
for (roi in rois){
	output_fig <- paste0(trimws(roi),'.ps')
	postscript(output_fig, width=1024, height=600, bg="white")
	myp1 <- ggplot(unc, aes_string( x="tracer", y=trimws(roi))) + geom_boxplot()
	print(myp1)
	dev.off()
}
$ Rscript comp_tracers.r
$ for x in braak_*.ps; do convert "${x}" -rotate 90 ${x%.ps}_mean_suvr.png; done
braak 1 braak 2
braak 3 braak 4
braak 5 braak 6

O todo junto,

$ echo "Subject,SUVR,Tracer,ROI" > header.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$2","$9",braak_1"}' > braak1.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$3","$9",braak_2"}' > braak2.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$4","$9",braak_3"}' > braak3.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$5","$9",braak_4"}' > braak4.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$6","$9",braak_5"}' > braak5.data
$ tail -n +2 tau_suvr_unc.csv | awk -F',' '{print $1","$7","$9",braak_6"}' > braak6.data
$ cat header.data braak1.data braak2.data braak3.data braak4.data braak5.data braak6.data > tau_suvr_unc_grouped.csv

y

comp_tau_all.r
library("ggplot2")
theme_set(theme_minimal())
input_file="tau_suvr_unc_grouped.csv"
unc <- read.csv(input_file)
output_fig <- 'output_tau_all_rois.ps'
postscript(output_fig, width=1024, height=600, bg="white")
pp <- ggplot(unc, aes_string( x="ROI", y="SUVR", fill="Tracer")) + geom_boxplot() + geom_hline(yintercept=1, color = "red", size=1)
print(pp)
dev.off()

O puedo comparar los valores de SUVR de tau con los valores de SUVR de FBB. lo mas facil es con los globales,

$ xnat_pullcl.pl -x f5cehbi -o f5cehbi_fbb_cl.csv
$ sed 's/;/,/g;s/Subject/PSubject/' f5cehbi_fbb_cl.csv > f5cehbi_fbb_cl_lone.csv
$ join -t, -1 2 -2 1 codes_sorted.csv f5cehbi_fbb_cl_lone.csv > f5cehbi_fbb_ok.csv
$ (head -n 1 f5cehbi_fbb_ok.csv && tail -n +2 f5cehbi_fbb_ok.csv | sort -t, -k 2) > f5cehbi_fbb_resorted.csv
$ join -t, -1 1 -2 2 f5cehbi_tau_suvr_pi2620_unc_tracer.csv f5cehbi_fbb_resorted.csv > comp_fbb_pi2620.csv
$ join -t, -1 1 -2 2 f5cehbi_tau_suvr_ro948_unc_tracer.csv f5cehbi_fbb_resorted.csv > comp_fbb_ro948.csv
$ tail -n +2 comp_fbb_ro948.csv > tmp.csv
$ cat comp_fbb_pi2620.csv tmp.csv > comp_fbb_tau.csv
comp_tracers.r
library("ggplot2")
theme_set(theme_minimal())
input_file="tau_suvr_unc.csv"
unc <- read.csv(input_file)
rois = c("braak_1","braak_2","braak_3","braak_4","braak_5","braak_6")
for (roi in rois){
	output_fig <- paste0(trimws(roi),'.ps')
	postscript(output_fig, width=1024, height=600, bg="white")
	myp1 <- ggplot(unc, aes_string( x="tracer", y=trimws(roi))) + geom_boxplot(aes(fill=factor(tracer)), show.legend=F)+labs(x="Radiotracer", y="SUVR", title=paste0("SUVR for ",roi))
	print(myp1)
	dev.off()
}
input_file="comp_fbb_tau.csv"
uac <- read.csv(input_file)
for (roi in rois){
	output_fig <- paste0('fbb_tau_',roi,'.ps')
	postscript(output_fig, width=1024, height=600, bg="white")
	myp1 <- ggplot(uac, aes_string(x=roi, y="SUVR", color="factor(tracer)")) + geom_point() + labs(x="tau SUVR", y="FBB SUVR", title=paste0("Comparing FBB Global SUVR and tau ",roi," SUVR"))
	myp1 <- myp1 + guides(color = guide_legend(title = "tracer"))
	print(myp1)
	dev.off
}
$ Rscript comp_tracers.r
$ for x in fbb_tau_braak_*.ps; do convert ${x} -rotate 90 ${x\.ps}.png; done
braak 1 braak 2
braak 3 braak 4
braak 5 braak 6

Con un poco de esfuerzo podemos usar el mismo metodo para medir el SUVR del FBB en las regiones de Braak.

Solo habria que reutilizar el procedimiento de PET-tau y engañarlo para que use los PET-FBB en su lugar. Como voy a suponer que ya he hecho el analisis de FBB, las imagenes puedo bajarmelas de XNAT, asi que me ahorro el registro a espacio nativo del PET.

Solo tengo que sacar la URI del fbb registrado y guardado y bajarme la imagen.

foreach my $subject (@pets){
        my $psubject = $pet_data{$subject};
        my $fake_tau = $w_dir.'/'.$subject.'_tau.nii.gz';
        #print "$fake_tau\n";
        my $xrd = xget_pet($xconf_data{'HOST'}, $jsession, $xprj, $psubject);
        if ($xrd){
                my $xld = 'curl -f -X GET -b "JSESSIONID='.$jsession.'" "'.$xconf_data{'HOST'}.'/data/experiments/'.$xrd.'/files?format=json" 2>/dev/null';
                #print "$xld\n";
                my $jres = qx/$xld/;
                my $xfres = decode_json $jres;
                foreach my $xres (@{$xfres->{'ResultSet'}{'Result'}}){
                        if ($xres->{'file_content'} eq 'PET_reg'){
                                my $xuri = $xres->{'URI'};
                                my $grd = 'curl -f -b "JSESSIONID='.$jsession.'" -X GET "'.$xconf_data{'HOST'}.$xuri.'" -o '.$fake_tau.' 2>/dev/null';
                                system($grd);
                        }
                }
        }
....
....
}

El resto del procedimiento es el mimo excepto que la region de referencia debe ser la materia gris del cerebelo. Por lo que debe ser mas rapido todo.

De momento voy a mirar solo los uncorrected.

pi2620 <- read.csv("f5cehbi_tau_suvr_pi2620_unc.csv")
fbb <- read.csv("f5cehbi_fbb_suvr_fake.csv")
ro948 <- read.csv("f5cehbi_tau_suvr_ro948_unc.csv")
merge(fbb, pi2620, by=c("Subject")) -> c1
for (roi in rois){
        output_fig <- paste0('fbb_pi2620_local_',roi,'.ps')
        postscript(output_fig, width=1024, height=600, bg="white")
        myp1 <- ggplot(c1, aes_string( x=paste0(roi,".y"), y=paste0(roi,".x"))) + geom_point() + labs(x="PI2620 SUVR", y="FBB local SUVR", title=paste0("Comparing FBB ",roi, " SUVR and PI2620 ",roi," SUVR"))
        print(myp1)
        dev.off()
}
merge(fbb, ro948, by=c("Subject")) -> c1

for (roi in rois){
        output_fig <- paste0('fbb_ro948_local_',roi,'.ps')
        postscript(output_fig, width=1024, height=600, bg="white")
        myp1 <- ggplot(c1, aes_string( x=paste0(roi,".y"), y=paste0(roi,".x"))) + geom_point() + labs(x="RO948 SUVR", y="FBB local SUVR", title=paste0("Comparing FBB ",roi, " SUVR and RO948 ",roi," SUVR"))
        print(myp1)
        dev.off()
}
$ Rscript comp_tracers.r
$ for x in *_local_*.ps; do convert ${x} -rotate 90 ${x%.ps}.png; done

RO948

braak 1 braak 2
braak 3 braak 4
braak 5 braak 6

PI2620

braak 1 braak 2
braak 3 braak 4
braak 5 braak 6

Si quiero comparar los tracers entre ellos,

$ join -t, -j 1 f5cehbi_tau_suvr_ro948_unc.csv f5cehbi_tau_suvr_pi2620_unc.csv > tau_ropi.csv
$ sed 's/,/ /g' tau_ropi.csv > tau_ropi.dat

y voy a pintar con gnuplot porque con R es demasiado complicado,

tau_plot.plt
set terminal postscript eps color enhanced size 10,7 "Times-Roman, 32"
set output "braak_1_comp.eps"
m(x) = a*x +b
fit m(x) "tau_ropi.dat" u 9:2 via a,b
unset key
f(x) = x
set xlabel "PI2620 SUVR" font "Times-Roman, 24"
set ylabel "RO948 SUVR" font "Times-Roman, 24"
set title "BRAAK 1 ROI"
plot [0.5:1.5][0.5:1.5] "tau_ropi.dat" u 9:2 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
 
set title "BRAAK 2 ROI"
set output "braak_2_comp.eps"
fit m(x) "tau_ropi.dat" u 10:3 via a,b
plot [0.2:1.1][0.2:1.1] "tau_ropi.dat" u 10:3 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
 
set title "BRAAK 3 ROI"
set output "braak_3_comp.eps"
fit m(x) "tau_ropi.dat" u 11:4 via a,b
plot [0.3:1.8][0.3:1.8] "tau_ropi.dat" u 11:4 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
 
set title "BRAAK 4 ROI"
set output "braak_4_comp.eps"
fit m(x) "tau_ropi.dat" u 12:5 via a,b
plot [0.6:2.0][0.6:2.0] "tau_ropi.dat" u 12:5 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
 
set title "BRAAK 5 ROI"
set output "braak_5_comp.eps"
fit m(x) "tau_ropi.dat" u 13:6 via a,b
plot [0.5:1.8][0.5:1.8] "tau_ropi.dat" u 13:6 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output
 
set title "BRAAK 6 ROI"
set output "braak_6_comp.eps"
fit m(x) "tau_ropi.dat" u 14:7 via a,b
plot [0.4:1.6][0.4:1.6] "tau_ropi.dat" u 14:7 w p ps 2 pt 7 lc 3, f(x) w l lw 5 lc 1, m(x) w l lc 4
set output

esto puede extenderse a la combinacion de ROIs

Para intentar caracterizar el off-target binding he descomentado la parte en que se calcula la retencion en eroded WM,

                #Hacer mascara de Eroded WM
                $ptask{'output'} = $outdir.'/tau_ewm_'.$subject;
                $ptask{'filename'} = $outdir.'/'.$subject.'_ewm.sh';
                $ptask{'command'} = $ENV{'PIPEDIR'}.'/bin/get_tref_ewm.sh '.$study.'_'.$subject.' '.$w_dir.'/.tmp_'.$subject;
                $mask_chain.= $w_dir.'/.tmp_'.$subject.'/rois/ewm.nii.gz ';
                send2slurm(\%ptask);

y entonces,

[osotolongo@brick03 f5cehbi]$ tau_proc.pl -tracer ro948 -r meta -p f5cehbi
...
...
...
[osotolongo@brick03 f5cehbi]$ for x in `find working/ -name "*_unc.csv"`; do id=$(echo ${x} | sed 's/.*\/\(.*\)_.*/\1/'); ewm=$(grep "^3" ${x} | awk '{print $2}'); icgm=$(grep "^4" ${x} | awk '{print $2}'); echo "${id} ${ewm} ${icgm}"; done | sort -n > ewm_icgm_ro948.csv
[osotolongo@brick03 f5cehbi]$ cat ewm_icgm_ro948.csv
0065 29.982110 20.723816
0074 20.095760 19.340429
0083 16.295862 13.545850
0089 35.377898 42.095356
0090 36.817361 45.491038
0092 39.262153 46.518434
0094 15.278639 22.369623
0097 21.845054 24.733724
0115 25.884322 34.812070
0118 20.809293 16.964745
0121 40.504989 36.016138
0141 30.426140 28.183759
0142 14.925220 15.207533
0149 42.510422 25.144353

Las columnas son: ID, EWM, ICGM

[osotolongo@brick03 f5cehbi]$ sed '1iid,ewm,icgm' ewm_icgm_ro948.csv | sed 's/ /,/g' > ewm_icgm_ro948_fixed.csv

Ahora hacemos lo mismo para pi2620 y,

[osotolongo@brick03 f5cehbi]$ for x in `find working/ -name "*_unc.csv"`; do id=$(echo ${x} | sed 's/.*\/\(.*\)_.*/\1/'); ewm=$(grep "^3" ${x} | awk '{print $2}'); icgm=$(grep "^4" ${x} | awk '{print $2}'); echo "${id} ${ewm} ${icgm}"; done | sort -n > ewm_icgm_pi2620.csv
[osotolongo@brick03 f5cehbi]$ awk '{print $0" r0948"}' ewm_icgm_ro948.csv | sed 's/ /,/g;1iid,ewm,icgm,tracer' > ewm_icgm_1.csv
[osotolongo@brick03 f5cehbi]$ awk '{print $0" pi2620"}' ewm_icgm_pi2620.csv | sed 's/ /,/g' > ewm_icgm_2.csv
[osotolongo@brick03 f5cehbi]$ cat ewm_icgm_1.csv ewm_icgm_2.csv > ewm_icgm.csv

Ahora, aqui hay una cosa rara. Estas son las regiones de referencia y supuestamente la relacion entre ambas deberia ser cercana a 1. Al ser menor que 1 implica que se esta reteniendo el trazador en el cerebelo.

neuroimagen/tau_analysis_refact.txt · Last modified: 2022/06/27 09:52 by osotolongo