metabcor.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate TACs of authentic (unchanged) tracer and radioactive",
  "metabolite(s) from measured plasma TAC and fractions of authentic tracer.",
  " ",
  "Usage: @P [Options] plasmafile fractionfile",
  " ",
  "Fraction file can have either of two formats:",
  " 1) it can contain the sample times (min), and the fraction(s) of authentic",
  "    tracer, and optionally the fractions of different metabolites; or",
  " 2) parameters of a mathematical function fitted to the fraction curves.",
  " ",
  "By default, result TACs are written in files named as *_pure.* and",
  "*_met.*, and if several metabolites exist, the following as *_met2.*",
  "and *_met3.*",
  " ",
  "Options:",
  " -fnpure=<filename>",
  "     Replace default filename parent tracer TAC.",
  " -fnmet=<filename>",
  "     Replace default filename for the plasma metabolite TAC.",
  " -fnmet2=<filename>",
  "     Replace default filename the second plasma metabolite TAC.",
  " -fnmet3=<filename>",
  "     Replace default filename  the 3rd plasma metabolite TAC.",
  " -pure=<id>",
  "     The result file is named by adding specified id text before file",
  "     extension, instead of the default '_pure'.",
  " -met=<id>",
  "     The result file is named by adding specified id text before file",
  "     name extension, instead of the default '_met'.",
  " -met2=<id>",
  "     The result file is named by adding specified id text before file",
  "     name extension, instead of the default '_met2'.",
  " -met3=<id>",
  "     The result file is named by adding specified id text before file",
  "     name extension, instead of the default '_met3'.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fit_ppf, fit_fexp, fit2dat, tac2svg, taccalc",
  " ",
  "Keywords: input, plasma, TAC, modelling, metabolite correction",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  int    ri, fi, fj, n, filetype, ret;
  double v, v_max;
  DFT    plasma, fract;
  FIT    fit;
  char  *cptr, tmp[FILENAME_MAX], plfile[FILENAME_MAX], ratfile[FILENAME_MAX];
  char   output_id0[128], output_id1[128], output_id2[128], output_id3[128];
  char  *outfile[4];
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  plfile[0]=ratfile[0]=(char)0;
  for(n=0; n<4; n++) outfile[n]=(char*)NULL;
  strcpy(output_id0, "_pure");
  strcpy(output_id1, "_met");
  strcpy(output_id2, "_met2");
  strcpy(output_id3, "_met3");
  dftInit(&plasma); dftInit(&fract); fitInit(&fit);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "PURE=", 5)==0) {
      cptr+=5; if(strlen(cptr)>0 && strlen(cptr)<128) {
        strlcpy(output_id0, cptr, 128); continue;}
    } else if(strncasecmp(cptr, "MET=", 4)==0) {
      cptr+=4; if(strlen(cptr)>0 && strlen(cptr)<128) {
        strlcpy(output_id1, cptr, 128); continue;}
    } else if(strncasecmp(cptr, "MET2=", 5)==0) {
      cptr+=5; if(strlen(cptr)>0 && strlen(cptr)<128) {
        strlcpy(output_id2, cptr, 128); continue;}
    } else if(strncasecmp(cptr, "MET3=", 5)==0) {
      cptr+=5; if(strlen(cptr)>0 && strlen(cptr)<128) {
        strlcpy(output_id3, cptr, 128); continue;}
    } else if(strncasecmp(cptr, "FNPURE=", 7)==0) {
      cptr+=7; if(strlen(cptr)>0 && strlen(cptr)<FILENAME_MAX) {
        outfile[0]=cptr; continue;}
    } else if(strncasecmp(cptr, "FNMET=", 6)==0) {
      cptr+=6; if(strlen(cptr)>0 && strlen(cptr)<FILENAME_MAX) {
        outfile[1]=cptr; continue;}
    } else if(strncasecmp(cptr, "FNMET2=", 7)==0) {
      cptr+=7; if(strlen(cptr)>0 && strlen(cptr)<FILENAME_MAX) {
        outfile[2]=cptr; continue;}
    } else if(strncasecmp(cptr, "FNMET3=", 7)==0) {
      cptr+=7; if(strlen(cptr)>0 && strlen(cptr)<FILENAME_MAX) {
        outfile[3]=cptr; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    if(!plfile[0]) {strcpy(plfile, argv[ai]); continue;}
    else if(!ratfile[0]) {strcpy(ratfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!ratfile[0]) {
    fprintf(stderr, "Error: missing file name for fraction data.\n");
    return(1);
  }
  if(verbose>5) MATHFUNC_TEST=verbose-5; else MATHFUNC_TEST=0;
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("plfile := %s\n", plfile);
    printf("ratfile := %s\n", ratfile);
    printf("output_id0 := %s\n", output_id0);
    printf("output_id1 := %s\n", output_id1);
    printf("output_id2 := %s\n", output_id2);
    printf("output_id3 := %s\n", output_id3);
    for(n=0; n<4; n++) if(outfile[n]) printf("outfile%d := %s\n", n, outfile[n]);
  }
 
 
  /*
   *  Read plasma TAC
   */
  if(verbose>1) printf("reading %s\n", plfile);
  if(dftRead(plfile, &plasma)) {
    fprintf(stderr, "Error in reading '%s': %s\n", plfile, dfterrmsg);
    return(2);
  }
  /* Delete other than the first column */
  if(plasma.voiNr>1) {
    fprintf(stderr, 
       "Warning: plasma file contains several TACs; only the first is used.\n");
    plasma.voiNr=1;
  }
  if(0) dftPrint(&plasma);
  /* Make sure that times are in minutes */
  /* Save the original unit, so that units can be converted back */
  if(plasma.timeunit==TUNIT_UNKNOWN) {
    plasma.timeunit=TUNIT_MIN;
    fprintf(stderr, "Warning: assuming time_unit := %s\n",
      petTunit(plasma.timeunit));
  }
  int orig_plasma_time_unit=plasma.timeunit;
  if(plasma.timeunit!=TUNIT_MIN) {
    if(verbose>2)
      fprintf(stdout, 
        "Note: plasma sample times are converted to minutes.\n");
    dftTimeunitConversion(&plasma, TUNIT_MIN);    
  }
  /* Make sure that plasma data has memory for unchanged tracer curve and
     all metabolite curves */
  ret=dftAddmem(&plasma, 4); if(ret) {
    fprintf(stderr, "Error in memory allocation.\n");
    dftEmpty(&plasma); return(2);
  }
  /* Set names for the TACs */
  if(strlen(plasma.voi[0].voiname)<1 || strcmp(plasma.voi[0].voiname, ".")==0)
    strcpy(plasma.voi[0].voiname, "Plasma");
  for(ri=1; ri<5; ri++) {
    strcpy(plasma.voi[ri].voiname, plasma.voi[0].voiname);
    strcpy(plasma.voi[ri].place, plasma.voi[0].place);
    plasma.voi[ri].size=plasma.voi[0].size;
  }
  cptr=output_id0; if(cptr[0]=='_') cptr++; strcpy(plasma.voi[1].hemisphere, cptr);
  cptr=output_id1; if(cptr[0]=='_') cptr++; strcpy(plasma.voi[2].hemisphere, cptr);
  cptr=output_id2; if(cptr[0]=='_') cptr++; strcpy(plasma.voi[3].hemisphere, cptr);
  cptr=output_id3; if(cptr[0]=='_') cptr++; strcpy(plasma.voi[4].hemisphere, cptr);
  for(ri=1; ri<5; ri++)
    sprintf(plasma.voi[ri].name, "%s %s", plasma.voi[ri].voiname, plasma.voi[ri].hemisphere);
  /* Remove NA's (only after dftAddmem()) */
  fi=0; n=plasma.frameNr;
  while(fi<plasma.frameNr) {
    if(!isnan(plasma.voi[0].y[fi])) {fi++; continue;}
    for(fj=fi+1; fj<plasma.frameNr; fj++) {
      plasma.x1[fj-1]=plasma.x1[fj];
      plasma.x2[fj-1]=plasma.x2[fj];
      plasma.x[fj-1]=plasma.x[fj];
      plasma.voi[0].y[fj-1]=plasma.voi[0].y[fj];
      plasma.w[fj-1]=plasma.w[fj];
    }
    plasma.frameNr--;
  }
  if(verbose>2) printf("%d non-available sample(s) removed from plasma TAC.\n",
                       n-plasma.frameNr);
  if(plasma.frameNr<1) {
    fprintf(stderr, "Error in %s: invalid contents.\n", plfile);
    dftEmpty(&plasma); return(2);
  }
 
 
  /*
   *  Read fraction curve(s)
   */
  if(verbose>1) printf("Reading fractions.\n");
  /* Check that file can be opened for read */
  {
    FILE *fp;
    fp=fopen(ratfile, "r");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open '%s'.\n", ratfile);
      dftEmpty(&plasma); return(3);
    }
    fclose(fp);
  }
  /* Determine file format */
  filetype=dftFormat(ratfile);
  if(filetype==DFT_FORMAT_UNKNOWN) {
    fprintf(stderr, "Error: unknown fraction file format.\n");
    dftEmpty(&plasma); return(3);
  }
  if(verbose>2) printf("Fractions filetype=%d\n", filetype);
 
  if(filetype==DFT_FORMAT_FIT) { /* function fitted to fractions */
 
    /* read FIT file */
    if(verbose>1) printf("reading %s\n", ratfile);
    if(fitRead(ratfile, &fit, verbose-5)) {
      fprintf(stderr, "Error in reading '%s': %s\n", ratfile, fiterrmsg);
      dftEmpty(&plasma); return(6);
    }
    if(verbose>10) fitPrint(&fit);
    /* Unknown time units are assumed to be in minutes */
    if(fit.timeunit==TUNIT_UNKNOWN) fit.timeunit=TUNIT_MIN;
    /* Check that file does not contain more than 4 fraction curves */
    if(fit.voiNr>4) {
      fprintf(stderr, "Error in '%s': too many metabolites.\n", ratfile);
      dftEmpty(&plasma); fitEmpty(&fit); return(6);
    }
    /* Check that fraction data extends to most of the plasma sample range */
    if(fit.voi[0].end>0.0 && plasma.x[plasma.frameNr-1]>1.33*fit.voi[0].end) {
      fprintf(stderr, "Warning: excessive extrapolation needed for fractions.\n");
    }
    /* Evaluate function values at plasma sample times */
    for(ri=0; ri<fit.voiNr; ri++) {
      ret=fitEvaltac(&fit.voi[ri], plasma.x, plasma.voi[ri+1].y2, plasma.frameNr);
      if(ret) {
        fprintf(stderr, "Error %d in function evaluation.\n", ret);
        dftEmpty(&plasma); fitEmpty(&fit); return(6);
      }
    }
    plasma.voiNr+=fit.voiNr;
    fitEmpty(&fit); /* no need for original fitted fractions any more */
    /* Check for percentages */
    for(fi=0, v_max=v=0.0; fi<plasma.frameNr; fi++) {
      for(ri=1, v=0.0; ri<plasma.voiNr; ri++) v+=plasma.voi[ri].y2[fi];
      if(v>v_max) v_max=v;
    }
    if(verbose>2) printf("v_max=%g\n", v_max);
    if(v_max>101.0) {
      fprintf(stderr, "Error: invalid fraction values.\n");
      dftEmpty(&plasma); dftEmpty(&fract); return(6);
    }
    if(v_max>1.01) {
      fprintf(stderr, "Warning: converting percentages to fractions.\n");
      for(fi=0; fi<plasma.frameNr; fi++) for(ri=1; ri<plasma.voiNr; ri++)
        plasma.voi[ri].y2[fi]/=100.0;
    }
 
  } else {   /* try to read as discrete fractions */
 
    /* read file */
    if(verbose>1) printf("reading %s\n", ratfile);
    if(dftRead(ratfile, &fract)) {
      fprintf(stderr, "Error in reading '%s': %s\n", ratfile, dfterrmsg);
      dftEmpty(&plasma); return(4);
    }
    /* Unknown time units are assumed to be in minutes */
    if(fract.timeunit==TUNIT_UNKNOWN) fract.timeunit=TUNIT_MIN;
    /* Check that file does not contain more than 4 fraction columns */
    if(fract.voiNr>4) {
      fprintf(stderr, "Error in '%s': too many metabolites.\n", ratfile);
      dftEmpty(&plasma); dftEmpty(&fract); return(4);
    }
    /* Check for percentages (before setting (0,1)) */
    for(fi=0, v_max=v=0.0; fi<fract.frameNr; fi++) {
      for(ri=0, v=0.0; ri<fract.voiNr; ri++)
        if(!isnan(fract.voi[ri].y[fi])) v+=fract.voi[ri].y[fi];
      if(v>v_max) v_max=v;
    }
    if(verbose>2) printf("v_max=%g\n", v_max);
    if(v_max>101.0) {
      fprintf(stderr, "Error: invalid fraction values.\n");
      dftEmpty(&plasma); dftEmpty(&fract); return(4);
    }
    if(v_max>1.011) {
      fprintf(stderr, "Warning: converting percentages to fractions.\n");
      for(fi=0; fi<fract.frameNr; fi++) for(ri=0; ri<fract.voiNr; ri++)
        if(!isnan(fract.voi[ri].y[fi])) fract.voi[ri].y[fi]/=100.0;
    }
    /* Assume that fractions at zero time are 1, 0, 0,... */
    if(fract.x[0]>0) {
      if(verbose>2) printf("adding zero fraction sample\n");
      ret=dftAddnullframe(&fract);
      if(ret) {
        fprintf(stderr, "Error in processing fractions.\n");
        dftEmpty(&plasma); dftEmpty(&fract); return(3);
      }
      fract.voi[0].y[0]=1.0;
    }
    /* Remove NA's in parent fractions, but do not allow NA's in 
       metabolite fractions */
    for(ri=n=0; ri<fract.voiNr; ri++) for(fi=0; fi<fract.frameNr; fi++) {
      if(isnan(fract.voi[ri].y[fi])) {
        if(ri>0) { // 0=parent
          fprintf(stderr, "Error: missing metabolite fraction(s).\n");
          dftEmpty(&plasma); dftEmpty(&fract); return(3);
        }
        n++;
      }
    }
    if(n>0) {
      if(verbose>1) printf("replacing NAs in fractions\n");
      ret=dftNAfill(&fract);
      if(ret) {
        fprintf(stderr, "Error in replacing missing fractions.\n");
        dftEmpty(&plasma); dftEmpty(&fract); return(3);
      }
    }
    if(verbose>10) dftPrint(&fract);
    /* Check that fraction data extends to most of the plasma sample range */
    if(plasma.x[plasma.frameNr-1]>1.5*fract.x[fract.frameNr-1]) {
      fprintf(stderr, "Error: excessive extrapolation needed for fractions.\n");
      dftEmpty(&plasma); dftEmpty(&fract); return(3);
    }
    if(plasma.x[plasma.frameNr-1]>1.2*fract.x[fract.frameNr-1]) {
      fprintf(stderr, "Warning: excessive extrapolation needed for fractions.\n");
    }
    /* Interpolate fraction curves to plasma sample times */
    if(verbose>1) printf("interpolation of fractions\n");
    for(ri=0; ri<fract.voiNr; ri++) {
      if(plasma.timetype==DFT_TIME_STARTEND)
        ret=interpolate4pet(fract.x, fract.voi[ri].y, fract.frameNr,
                 plasma.x1, plasma.x2, plasma.voi[ri+1].y2, NULL, NULL, 
                 plasma.frameNr);
      else
        ret=interpolate(fract.x, fract.voi[ri].y, fract.frameNr,
                 plasma.x, plasma.voi[ri+1].y2, NULL, NULL, plasma.frameNr);
      if(ret) {
        fprintf(stderr, "Error in interpolation.\n");
        dftEmpty(&plasma); dftEmpty(&fract); return(4);
      }
    }
    plasma.voiNr+=fract.voiNr;
    dftEmpty(&fract); /* no need for original fractions any more */
  }
  if(verbose>8) dftPrint(&plasma);
 
  /*
   *  Calculate metabolite fractions, if...
   */
  if(plasma.voiNr==2) { /* only unchanged fraction was given */
    if(verbose>1) 
      printf("calculating metabolite fractions for extra metabolite\n");
    for(fi=0; fi<plasma.frameNr; fi++)
      plasma.voi[2].y2[fi]=1.0-plasma.voi[1].y2[fi];
    plasma.voiNr=3;
  } else {
    /* Check if the sum of fractions is clearly <1 */
    if(verbose>1) printf("checking the sum of fractions\n");
    for(fi=0, v_max=0.0, n=0; fi<plasma.frameNr; fi++) {
      for(ri=1, v=0.0; ri<plasma.voiNr; ri++) v+=plasma.voi[ri].y2[fi];
      if(plasma.voiNr<5) plasma.voi[plasma.voiNr].y2[fi]=1.0-v;
      if(v>0.0) n++; else continue;
      if((1.0-v)>v_max) v_max=1.0-v;
    }
    if(verbose>2 && (v_max>0 || n>0))
      printf("sum of fractions was less than 1.0; vmax=%g n=%d\n", v_max, n);
    if(v_max>0.02 && n>1) { /* yes there are, this is not a rounding up error */
      if(plasma.voiNr>4) { /* there's no room for more */
        fprintf(stderr, "Error: sum of fractions was less than 1.0.\n");
        dftEmpty(&plasma); return(7);
      }
      printf("sum of fractions was less than 1.0; assuming another metabolite.\n");
      plasma.voiNr++;
    }
  }
  if(verbose>7) dftPrint(&plasma);
 
 
  /*
   *  Calculate the curves of authentic tracer and metabolites
   */
  for(ri=1; ri<plasma.voiNr; ri++) {
    /* Multiply fractions with total plasma concentration */
    for(fi=0; fi<plasma.frameNr; fi++)
      plasma.voi[ri].y[fi]=plasma.voi[0].y[fi]*plasma.voi[ri].y2[fi];
  }
 
  /* Change times back to original units if necessary */
  if(plasma.timeunit!=orig_plasma_time_unit) {
    if(verbose>3) printf("converting plasma %s to %s\n",
      petTunit(plasma.timeunit), petTunit(orig_plasma_time_unit) );
    if(verbose>4) printf("last_t := %g\n", plasma.x[plasma.frameNr-1]);
    ret=dftTimeunitConversion(&plasma, orig_plasma_time_unit);
    if(verbose>4) printf("last_t := %g\n", plasma.x[plasma.frameNr-1]);
  }
 
 
  /*
   *  Save the curves of authentic tracer and metabolites
   */
  n=plasma.voiNr;
  for(ri=1; ri<n; ri++) {
    /* Copy one curve at a time to the first place */
    ret=dftCopyvoi(&plasma, ri, 0);
    if(ret) {
      fprintf(stderr, "Error in processing TACs.\n");
      dftEmpty(&plasma); return(9);
    }
    /* Construct filename */
    if(outfile[ri-1]) {
      strcpy(tmp, outfile[ri-1]);
    } else {
      strcpy(tmp, plfile); cptr=strrchr(tmp, '.'); if(cptr!=NULL) *cptr=(char)0;
      switch(ri) {
        case 1: strcat(tmp, output_id0); break;
        case 2: strcat(tmp, output_id1); break;
        case 3: strcat(tmp, output_id2); break;
        case 4: strcat(tmp, output_id3); break;
        default:
          fprintf(stderr, "Error in programmer.\n");
          dftEmpty(&plasma); return(99);
      }
      cptr=strrchr(plfile, '.'); if(cptr!=NULL) strcat(tmp, cptr);
    }
    /* Save TAC */
    if(verbose>0) fprintf(stdout, " writing %s\n", tmp);
    plasma.voiNr=1; ret=dftWrite(&plasma, tmp); plasma.voiNr=n;
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, tmp, dfterrmsg);
      dftEmpty(&plasma); return(11);
    }
  }
 
  /* Free memory */
  dftEmpty(&plasma);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/