source: tags/initial/GDE/PHYLIP/makeinf.c

# include <stdio.h>

char aliname[80], nucname[80], destname[80], taxname[11][255];
FILE *sourceali, *sourcenuc, *destali;
        /* files: their names and pointers
        */
int totnum, seqnum, blocknum, amslr, dpvou,
             nucoraa, taxnum, globcount, position;
        /* totnum is the total number of sequences in the alignment
           seqnum is the number of sequences that will be used
           counters that get incremented or decremented:
            amslr, bntms, count, dpvou, etc
           nucoraa: 1 for nuc alignments, 3 for amino acid alignments
           taxnum is the ID-number of the current taxon
           globcount is the counter for the main loop
                position is the number for the codon position to be analyzed
        */
int oldcheck, newcheck;   /* for checking the number of characters */
int nucmito, yorm, namebool;
        /* for nuclear genetic code: nucmito == 0
           for mitochondrial code : nucmito == 1
           if yorm == 0: no conversion of leu and arg first positions;
           if yorm == 1: conversion occurs
        */
int putnucbool;
        /* for exclusion of regions where homology is uncertain */
long begin, remember;
        /* current pointer position within sourceali and destali */
double  actualcount, count[9];


main()
{
long temp;                  /* holds a temporary position of pointer
                               within destali */

/* INITIALIZING GLOBAL VARIABLES */

position = 5;                   /* initialize to analyze all positions */
putnucbool = 0;                 /* initialize exclude-boolean */
actualcount = 0;                /* initialize base counter */
yorm = 2;                       /* initialize 'Y' or 'M' option */
nucmito = 2;                    /* initialize type of genetic code option */
for (amslr = 0; amslr < 9; amslr++) {
        count[amslr] = 0;
}                               /* initialize base counters */

/* FUNCTION CALLS START HERE */

parstuff();                     /* gets info from user */
findalignment();                /* finds the string ALIGNMENT */
findbeginning();                /* finds the start of the sequence alignment */
blocknum = countblocks();       /* counts the number of blocks in the
                                 * alignment */
goback();                       /* pointer back to last \n before alignment */

/* MAIN LOOP. ONE EXECUTION OF IT PROCESSES ONE SEQUENCE */

for (globcount = 0; globcount < seqnum; globcount++) {
        if (globcount > 0) {    /* if it's not the first sequence */
                jump(1);        /* puts pointer at beginning of sequence line */
                begin = ftell(sourceali);
        }
        findname();             /* finds the number of this sequence in
                                 * sourceali */
        goback();               /* goes back to beginning of line */
        findnucseq();           /* finds the corresponding sequence in
                                 * sourcenuc */

        for (dpvou = 0; dpvou < blocknum; dpvou++)
                /*
                 * loops through the alignment as many times as there are
                 * blocks
                 */
        {
                putnucs();      /* grabs nucleotides, puts them into destali */
                if (dpvou != (blocknum - 1))    /* if it's not the last
                                                 * block.. */
                        jump(seqnum - 1);       /* ..jump */
                else if (dpvou == (blocknum - 1)) {
                        tailcheck(globcount);   /* checks if there are bases
                                                 * remaining beyond the end
                                                 * of the sequence */
                        goback();       /* go back to beginning of previous
                                         * sequence in alignment */
                }
        }                       /* end of single sequence processing */

        if (globcount == 0) {
                oldcheck = newcheck;    /* length of first sequence is
                                         * standard */
                temp = ftell(destali);
                fseek(destali, remember, 0);
                fprintf(destali, "%10d", oldcheck);
                fseek(destali, temp, 0);
        } else if (oldcheck != newcheck) {
                printf("CRASH! Lengths of sequences do not agree. Check %s\n", aliname);
                exit(0);
        }
        if ((newcheck % 60) != 0)
                fputc('\n', destali);
        newcheck = 0;           /* resets counter */
}

/* END OF MAIN LOOP */

/* CALCULATE BASE COMPOSITION IF FIRST POSITIONS OF CODONS ARE INCLUDED */
if ((nucoraa == 3) && (position != 2) &&
    (position != 3) && (yorm == 1))
        summer();

fclose ( destali );
fclose ( sourcenuc );
fclose ( sourceali );

return 0;
}        /* end of main */


parstuff()
{                               /* gets info from user */
        int             bool, dumnum;
        char            dummy, nucchar;

        do {
                printf("Alignment file to be read:      ");
                scanf("%s", aliname);
                if (NULL == (sourceali = fopen(aliname, "r"))) {
                        bool = 0;
                        printf("Cannot open %s. Check path or spelling.\n", aliname);
                } else
                        bool = 1;
        }
        while (bool == 0);

        do {
                printf("Nucleotide file to be read:     ");
                scanf("%s", nucname);
                if (NULL == (sourcenuc = fopen(nucname, "r"))) {
                        bool = 0;
                        printf("Cannot open %s. Check path or spelling.\n", nucname);
                } else
                        bool = 1;
        }
        while (bool == 0);

        do {
                printf("Destination file to be written: ");
                scanf("%s", destname);
                if (NULL == (destali = fopen(destname, "w"))) {
                        bool = 0;
                        printf("Cannot open %s. Check path or spelling.\n", destname);
                } else
                        bool = 1;
        }
        while (bool == 0);

        do {
                printf("\nTotal number of sequences in alignment: ");
                scanf("%d", &totnum);
                if (totnum < 2) {
                        printf("%d; CRASH! That won't work.\n", &totnum);
                        exit(0);
                } else
                        bool = 1;
        }
        while (bool == 0);
        scanf("%c", &dummy);

        do {
                printf("Number of sequences to be used:         ");
                scanf("%d", &seqnum);
                if (seqnum < 2) {
                        printf("%d; CRASH! That won't work.\n", &seqnum);
                        exit(0);
                } else
                        bool = 1;
        }
        while ((bool == 0) || (seqnum > totnum));
        scanf("%c", &dummy);

        fprintf(destali, "  %d", seqnum);       /* Prints the number of
                                                 * sequences.. */
        remember = ftell(destali);      /* ..remembers where the pointer is.. */
        fprintf(destali, "    LENGTH\n");       /* ..and keeps ten spaces for
                                                 * the length to be filled in
                                                 * later  */
        do {
                printf("Nucleic acid or Protein coding sequence?          (n/p): ");
                scanf("%c", &nucchar);
                if ((nucchar == 'n') || (nucchar == 'N')) {
                        nucoraa = 1;
                        bool = 1;
                } else if ((nucchar == 'p') || (nucchar == 'P')) {
                        nucoraa = 3;
                        bool = 1;
                } else {
                        printf("Try again:\n");
                        scanf("%c", &dummy);
                        bool = 0;
                }
        }
        while (bool == 0);
        scanf("%c", &dummy);

        if (nucoraa == 3) {
                do {
                        printf("Nuclear or mitochondrial genetic code?            (n/m): ");
                        scanf("%c", &nucchar);
                        if ((nucchar == 'n') || (nucchar == 'N')) {
                                nucmito = 0;
                                bool = 1;
                        } else if ((nucchar == 'm') || (nucchar == 'M')) {
                                nucmito = 1;
                                bool = 1;
                        } else {
                                printf("Try again:\n");
                                scanf("%c", &dummy);
                                bool = 0;
                        }
                }
                while (bool == 0);
                scanf("%c", &dummy);

                do {
                        printf("\nEnter a number between 1 and 5, for the\n");
                        printf("   codon position you wish to analyze:\n");
                        printf("   1 for first, 2 for second, 3 for third\n");
                        printf("   4 for first plus second, 5 for all.            (1-5): ");
                        scanf("%d", &position);
                        if ((position < 1) || (5 < position)) {
                                printf("Bull....  Try again\n");
                        }
                }
                while ((position < 1) || (5 < position));
                scanf("%c", &dummy);

                if ((position != 2) && (position != 3)) {
                        do {
                                printf("Conversion of first positions to degenerate base? (y/n): ");
                                scanf("%c", &nucchar);
                                if ((nucchar == 'y') || (nucchar == 'Y')) {
                                        yorm = 1;
                                        bool = 1;
                                } else if ((nucchar == 'n') || (nucchar == 'N')) {
                                        yorm = 0;
                                        bool = 1;
                                } else {
                                        printf("Try again:\n");
                                        scanf("%c", &dummy);
                                        bool = 0;
                                }
                        }
                        while (bool == 0);
                        scanf("%c", &dummy);
                }
        }
        do {
                printf("Use nAmes or nUmbers as identifiers?              (a/u): ");
                scanf("%c", &nucchar);
                if ((nucchar == 'a') || (nucchar == 'A')) {
                        namebool = 0;
                        bool = 1;
                } else if ((nucchar == 'u') || (nucchar == 'U')) {
                        namebool = 1;
                        bool = 1;
                } else {
                        printf("Try again:\n");
                        scanf("%c", &dummy);
                        bool = 0;
                }
        }
        while (bool == 0);
        scanf("%c", &dummy);

        printf("\nNucleotide sequences in:           %12s\n", nucname);
        if (nucoraa == 1) {
                printf("Nucleotide alignment source:       %12s\n", aliname);
                printf("Nucleotide alignment destination:  %12s\n", destname);
        } else if (nucoraa == 3) {
                printf("Amino acid alignment source:       %12s\n", aliname);
                printf("Nucleotide alignment destination:  %12s\n", destname);
                if (position == 1) {
                        printf("First ");
                } else if (position == 2) {
                        printf("Second ");
                } else if (position == 3) {
                        printf("Third ");
                } else if (position == 4) {
                        printf("First plus second ");
                } else if (position == 5) {
                        printf("All ");
                }
                printf("codon positions will be used.\n");

                if ((yorm == 0) && (position != 2) && (position != 3)) {
                        if (nucmito == 0) {
                                printf("No conversion of L and R 1st positions.\n");
                        } else if (nucmito == 1) {
                                printf("No conversion of L 1st positions.\n");
                        }
                } else if ((yorm == 1) && (position != 2) && (position != 3)) {
                        if (nucmito == 0) {
                                printf("L and R 1st positions will be converted to Y and M.\n");
                        } else if (nucmito == 1) {
                                printf("L 1st positions will be converted to Y.\n");
                        }
                }               /* end of if yorm == 0 */
        }                       /* end of nucoraa == 3 */
        if (namebool == 0) {
                printf("Names ");
        } else if (namebool == 1) {
                printf("Numbers ");
        }
        printf("will be used to identify sequences.\n\n");

        return 0;
}


findalignment()
{                               /* finds the string "ALIGNMENT" in the
                                 * alignment */
        char            testchar[9], dummy, seqname[200];
        int             bntms, count, eqwpv;

        for (count = 0; count < totnum; count++) {
                fscanf(sourceali, "%d", &taxnum);
                dummy = fgetc(sourceali);
                fgets(seqname, 200, sourceali);
                bntms = 0;
                while ((seqname[bntms] != '\n') && (bntms < 10)) {
                        if (seqname[bntms] != '\n') {
                                taxname[bntms][taxnum] = seqname[bntms];
                        }
                        bntms++;
                }
                if (bntms < 10) {
                        for (eqwpv = bntms; eqwpv < 10; eqwpv++)
                                taxname[eqwpv][taxnum] = ' ';
                }
        }

        testchar[0] = '\n';
        while (strncmp(testchar, "ALIGNMENT", 9) != 0) {
                for (count = 0; count < 9; count++) {
                        testchar[count] = fgetc(sourceali);
                        if (testchar[count] == '\n') {
                                break;
                        }
                }
                dummy = testchar[count];
                if (dummy == '\n') {
                        continue;
                }
                while ((dummy = fgetc(sourceali)) != '\n');
        }
        return 0;
}

findbeginning()
{                               /* finds the first line of sequence alignment */
        char            testchar;

        while ((testchar = fgetc(sourceali)) != '\n');
        while ((testchar != 45 && testchar < 65) ||
               (testchar > 91 && testchar != '{'))
                testchar = fgetc(sourceali);
        begin = ftell(sourceali) - 1;
        /* pointer will be put to last \n before first amino acid */
        return 0;
}

countblocks()
{                               /* counts the number of seqeunce blocks in
                                 * the alignment */
        char            dummy;
        int             dumnuc;

        dumnuc = -1;
        do {
                dummy = fgetc(sourceali);
                while ((dummy != '*') && (dummy != EOF))
                        dummy = fgetc(sourceali);
                dumnuc++;
                while ((dummy != '\n') && (dummy != EOF))
                        dummy = fgetc(sourceali);
        }
        while (dummy != EOF);
        return dumnuc;
}

jump(howmany)
{                               /* jumps howmany lines in sourceali */
        char            linedump, aminocheck;
        int             eqwpv;
        long            lastpos;
        eqwpv = 0;

        do {
                linedump = fgetc(sourceali);
                /*
                 * checks if the first character on the line is an amino acid
                 * or a [
                 */
                if ((linedump < 65 || 91 < linedump) && (linedump != 45))
                        /* if it's not an amino acid.. */
                {
                        --eqwpv;/* ..do not increase counter */
                        if (linedump == '{')    /* if it's an ignore
                                                 * character.. */
                                while ((linedump = fgetc(sourceali)) != '}');
                        /* ..go to resume char */
                }
                if (linedump != '\n')   /* if not at end of line.. */
                        while ((linedump = fgetc(sourceali)) != '\n');
                /* ..goes to end of line */
                eqwpv++;
        }
        while (eqwpv < howmany);

        do {                    /* this loop ties up the end by ensuring that
                 * the next line starts with an amino acid or a [ */
                lastpos = ftell(sourceali);     /* remembers the position in
                                                 * the sequence */
                aminocheck = fgetc(sourceali);  /* amino acid, hyphen, or '['
                                                 * ? */
                if ((aminocheck < 65 || 91 < aminocheck) && aminocheck != 45) {
                        if (aminocheck == '{')  /* if it's an ignore
                                                 * character.. */
                                while ((aminocheck = fgetc(sourceali)) != '}');
                        /* ..goes to resume char */
                        if (aminocheck != '\n') /* if not at end of line ... */
                                while ((aminocheck = fgetc(sourceali)) != '\n');
                        /* ..goes to end of line */
                }
        }
        while ((aminocheck < 65 || 91 < aminocheck) && (aminocheck != 45));

        fseek(sourceali, lastpos, 0);   /* loop is ended by amino acid; now
                                         * go back to the beginning of the
                                         * current line */
        return 0;
}                               /* end of jump */


findname()
{                               /* finds the number of the sequence line that
                                 * it's looking at and prints it or the
                                 * corresponding name at the beginning of the
                                 * destfile's corresponding sequence */
        char            dummy;
        int             count, totnumloop, testcount, namefield;
        long            currpos;

        amslr = 0;
        dummy = fgetc(sourceali);
        if (dummy == '{') {     /* if it's an ignore char.. */
                do
                        dummy = fgetc(sourceali);
                while (dummy != '}');
                do
                        dummy = fgetc(sourceali);
                while (dummy != '\n');
                begin = ftell(sourceali);       /* for putnucs (\n) */
                dummy = fgetc(sourceali);       /* pointer on first aa */
        }                       /* ..read until on the first amino acid after
                                 * the resume char */
        while (dummy == 32 || dummy == 45 ||
               (65 <= dummy && dummy <= 91) || dummy == 93) {
                dummy = fgetc(sourceali);
                amslr++;
        }

/* USER: IF YOUR ALIGNMENT ONLY HAS THE TAXON NUMBER AT THE END OF EACH
 * LINE, AND NOT ALSO THE POSITION NUMBER, DELETE ONE OF THE FOLLOWING
 * TWO LINES BEFORE COMPILING! */

        fscanf(sourceali, "%d", &taxnum);       /* dumps the length number */
        fscanf(sourceali, "%d", &taxnum);       /* correct number of taxon */

        if (namebool == 1) {    /* if use numbers */
                fprintf(destali, "%d", taxnum);
                if (taxnum < 10) {
                        testcount = 1;
                }
                 /* writes the number ... */
                else if ((taxnum < 100) && (taxnum > 9)) {
                        testcount = 2;
                } else if ((taxnum < 512) && (taxnum > 99)) {
                        testcount = 3;
                } else {
                        printf("Number of sequence out of range.\n");
                        exit(0);
                }

                for (namefield = testcount; namefield < 10; namefield++) {
                        fputc(' ', destali);
                }               /* ..and fills up the namefield to 10 chars */
        } else if (namebool == 0) {     /* if use names, not numbers.. */
                for (count = 0; count < 10; count++) {
                        fprintf(destali, "%c", taxname[count][taxnum]);
                }               /* ..writes the name */
        }
        while (dummy != '\n') {
                dummy = fgetc(sourceali);
        }

        return 0;
}

findnucseq()
{                               /* first puts file-pointer at > of sequence
                                 * and then puts it at first nucleotide */
        char            dummy;
        int             countseq;

        fseek(sourcenuc, 0, 0);

        countseq = -1;          /* sequence numbers start with 0 */
        while (countseq < taxnum) {
                countseq++;
                while ((dummy = fgetc(sourcenuc)) != '>');      /* marker */
                while ((dummy = fgetc(sourcenuc)) != '\n');     /* marker line end */
                while ((dummy = fgetc(sourcenuc)) != '\n');     /* name line end */
        }
        return 0;
}


goback()
{
        fseek(sourceali, begin, 0);
        return 0;
}


countnuc(nuc)
{
        actualcount++;
        switch (nuc) {
        case 'G':
                count[0]++;
                break;
        case 'A':
                count[1]++;
                break;
        case 'T':
        case 'U':
                count[2]++;
                break;
        case 'C':
                count[3]++;
                break;
        case 'Y':
                count[4]++;
                break;
        case 'R':
                count[5]++;
                break;
        case 'M':
                count[6]++;
                break;
        case 'N':
                count[7]++;
                break;
        default:
                count[8]++;
                break;
        }
        return 0;
}


putnucs()
{
char            base, dummy;
int             zahl;

while ((dummy = fgetc(sourceali)) != '\n') {
 if (dummy == 91) {     /* if it's a start-exclude character ([) */
  putnucbool = 1;
 } else if (dummy == 93) { /* if it's a stop-exclude character (]) */
  putnucbool = 0;
 } else if (65 <= dummy && dummy <= 90) { /* if it's an authentic
                                           * sequence character */
  for (zahl = 0; zahl < nucoraa; zahl++) {
   base = fgetc(sourcenuc);
   if ((65 <= base && base <= 90) || (97 <= base && base <= 122)) {
    if (putnucbool == 0) {
     if (nucoraa == 3) { /* if it's an amino acid sequence */
      if (yorm == 1) { /* if conversion is desired */
       if ((dummy == 'L') && (zahl == 0)) {
         base = 'Y'; /* if amino acid is leucine */
       }
       if (nucmito == 0) { /* if code is nuclear */
        if ((dummy == 'R') && (zahl == 0)) {
         base = 'M'; /* if amino acid is arginine */
        }
       }
      }
      if (97 <= base && base <= 122) {
       base -= 32;
      }
     }
     if ((position == 5) || ((position == 1) && (zahl == 0)) ||
        ((position == 2) && (zahl == 1)) ||
        ((position == 3) && (zahl == 2)) ||
        ((position == 4) && ((zahl == 0) || (zahl == 1)))) {
      countnuc(base);
      newcheck++;
      fputc(base, destali);
      if (newcheck % 60 == 0)
       fputc('\n', destali);
     }
    } /* end of if conditional on putnucbool */
   } /* end of if conditional on authentic nucleotide */
   else if (base == '*') {
    printf("CRASH! not enough nucleotides.\n");
    exit(0);
   } else
     --zahl;
  } /* end of loop counted by nucoraa and zahl */
 } /* end of if conditional on authentic sequence character in alignment */
 else if (dummy == 45) { /* if it is a hyphen */
  if (putnucbool == 0) {
   for (zahl = 0; zahl < nucoraa; zahl++) {
    if ((position == 5) || ((position == 1) && (zahl == 0)) ||
       ((position == 2) && (zahl == 1)) ||
       ((position == 3) && (zahl == 2)) ||
       ((position == 4) && ((zahl == 0) || (zahl == 1)))) {
     newcheck++;
     fputc('-', destali);
     if (newcheck % 60 == 0)
      fputc('\n', destali);
    }
   }
  }
 }
}  /* end of while */
return 0;
}  /* end of putnucs */


tailcheck(seq)
{
        char            checkbool, pacman;
        checkbool = 0;
        do {
                pacman = fgetc(sourcenuc);
                if ((65 <= pacman && pacman <= 90) ||
                    (97 <= pacman && pacman <= 122)) {
                        checkbool = 1;
                        if (checkbool == 1)
                                printf("%c", pacman);
                }
        }
        while (pacman != '*');
        if (checkbool == 1) {
                printf("\nWARNING. Sequence # %d in %s contains the above\n", seq, nucname);
                printf(" bases beyond the number of bases it should contain\n");
                printf(" according to the alignment in %s.\n", aliname );
                printf(" Please compare the two files.\n\n");
        }
}

summer()
{
        int             cou;

        count[1] += (2 * count[6] / 3);
        count[2] += (count[4] / 3);
        count[3] += ((2 * count[4] + count[6]) / 3);
        printf("Frequencies of A, C, G, T:\n");
        printf("%10.5lf", count[1] / actualcount);
        printf("%10.5lf", count[3] / actualcount);
        printf("%10.5lf", count[0] / actualcount);
        printf("%10.5lf\n", count[2] / actualcount);
        return 0;
}
/* end of program */