/***************************************************************************** ****************************************************************************** ****************************************************************************** ***** ***** ***** PROGRAM CALCD v1.0 ***** ***** Copyright (c) 1990 by Joyce C. Miller. All Rights Reserved. ***** ***** A program to calculate genetic distance from equations ***** ***** from equations 5.53 - 5.55 of Nei (1987). ***** ***** ***** ***** This program takes in m(x), m(y), m(xy), and the number of base ***** ***** pairs in the restriction site, and calculates F (the proportion of ***** ***** fragments shared), and d (the distance). It uses equations ***** ***** 5.53 - 5.55 of Nei (1987). ***** ***** Data entry can be of two types: batch mode, where the data is ***** ***** in an ASCII file, or interactive, where it is entered by the user. ***** ***** In batch mode, the user must enter the name of an ASCII text file ***** ***** which holds the data, and the name of an output file that the F's ***** ***** and d's are to be written to. The ASCII text file should have the ***** ***** following format: ***** ***** Mx My Mxy r ex: 50 50 25 6 ***** ***** Mx My Mxy r 60 50 35 4 ***** ***** etc. etc. ***** ***** If no input and output file names are entered on the command line, ***** ***** the program goes into interactive mode. In this case, the user is ***** ***** asked for Mx, My, Mxy, and r, and F-hat and d-hat are printed to ***** ***** the screen. ***** ***** ***** ***** Reference: ***** ***** Nei, M. (1987) Molecular Evolutionary Genetics. Columbia ***** ***** University Press, New York. ***** ***** ***** ***** ***** ***** List of C functions used in this program: ***** ***** ***** ***** FUNCTION LIBRARY FUNCTION LIBRARY ***** ***** fabs math.h fclose stdio.h ***** ***** fgets stdio.h fopen stdio.h ***** ***** fprintf stdio.h log math.h ***** ***** printf stdio.h rewind stdio.h ***** ***** scanf stdio.h sqrt math.h ***** ***** sscanf stdio.h strlen string.h ***** ***** toupper ctype.h ***** ****************************************************************************** ****************************************************************************** *****************************************************************************/ #include #include #include #include #include #include /***************************************************************************** ** SYMBOLIC CONSTANTS ** *****************************************************************************/ #define VNUM 1.0 /* version number of this program */ #define YEAR 1990 /* calendar year this version was released */ #define TEST 0.000001 /* value used for testing end of iterative loop */ /***************************************************************************** ** TYPE DEFINITIONS ** *****************************************************************************/ typedef struct { /* structure to hold Mx, My, Mxy, and r */ float mx; /* number of fragments in OTU X */ float my; /* number of fragments in OTU Y */ float mz; /* number of fragments in common between OTU X and OTU Y */ float rv; /* number of base pairs in restriction site */ } ms; /***************************************************************************** ** FUNCTION PROTOTYPES ** *****************************************************************************/ void batch(int argc, char *argv[]); /* execute program in batch mode */ void getflnames(int argc, char *argv[]); /* get input & output file names */ FILE *computeFnd(FILE *fpw); /* compute F & d & print to file (batch mode) */ void screen(); /* execute program in interatrive mode */ void getm(); /* get Mx, My, Mxy & r from user */ /***************************************************************************** ** GLOBAL VARIABLES ** *****************************************************************************/ char infl[20]; /* name of input file */ char outfl[20]; /* name of output file */ ms m; /* structure to hold data */ /***************************************************************************** ****************************************************************************** ****************************************************************************** ***** ***** ***** MAIN PROGRAM ***** ***** ***** ****************************************************************************** ****************************************************************************** ****************************************************************************** ***** ***** ***** This section of the program looks at the command line, and decides ***** ***** if it is to execute a batch file, or to go into interactive mode, ***** ***** and begin asking for data. ***** ***** ***** ***** Functions called: ***** ***** batch -- computes F and d from data in ASCII text file. ***** ***** screen -- calculates F and d from data input by user. ***** *****************************************************************************/ main(int argc, char *argv[]) { printf("\r\n\nProgram CALCD v%3.1f\r\n",VNUM); printf("Copyright (c) %4.0d by Joyce C. Miller. ",YEAR); printf("All Rights Reserved.\r\n\n"); if (argc > 1) batch(argc,argv); /* execute as batch file */ else screen(); /* get data from user */ } /* END OF FUNCTION MAIN */ /***************************************************************************** ** ** ** FUNCTION BATCH ** ** ** ** This function is called from FUNCTION MAIN, and allows the program to ** ** execute as a batch file. First, the input file name and output file ** ** names are extracted from the command lines, the the input file is ** ** opened, and the output file is created. Then, each line in the input ** ** file is read. On each line, Mx, My, Mxy, and the r-value areread in, ** ** and then FUNCTION COMPUTEFND is called to compute F-hat and d-hat, and ** ** print it to the output file. ** ** ** ** Functions called: ** ** getflnames -- get input and output file names from command line. ** ** opntrdfl -- open text file for reading. ** ** opntwtfl -- open (create text file for writing. ** ** computeFnd -- compute F and d, and print them to file. ** *****************************************************************************/ void batch(int argc, char *argv[]) { FILE *fpr,*fpw; /* input and output file pointers */ char tempstr[255]; /* string read in from input file */ getflnames(argc,argv); /* get input file and output file */ fpr = opntrdfl(fpr,infl); /* open input file */ fpw = opntwtfl(fpw,outfl); /* open output file */ while ((fgets(tempstr,255,fpr)) != NULL) { /* get a string */ m.mx = 0; /* initialize structure to hold N's */ m.my = 0; m.mz = 0; m.rv = 0; if (strlen(tempstr) > 7) { /* get Mx, My, Mxy, and r */ sscanf(tempstr,"%f %f %f %f",&m.mx,&m.my,&m.mz,&m.rv); fpw = computeFnd(fpw); /* compute F & d & print to file */ } } rewind(fpr); /* close up input file */ fclose(fpr); rewind(fpw); /* close up output file */ fclose(fpw); } /* END OF FUNCTION BATCH */ /***************************************************************************** ** ** ** FUNCTION GETFLNAMES ** ** ** ** This function is called from FUNCTION BATCH, and reads the input and ** ** output file names from the command line. ** ** ** ** Functions called: none. ** *****************************************************************************/ void getflnames(int argc, char *argv[]) { int i; /* loop control variable */ if (argc < 3) { /* if either infl or outfl name is missing */ printf("\r\nError: Input and/or output file name is missing. "); printf("Both are required.\r\n"); exit(0); } else { /* get infl and outfl names from command line */ for (i=0; i (TEST*G2)); ++i) G1 = G2; /* compute G2 via eqn. 5.54 */ d = -(2*log(G2))/m.rv; /* compute d via eqn. 5.55 */ fprintf(fpw,"F = %8.6f d = %8.6f\n",F,d); /* print F-hat and d-hat */ } return(fpw); /* return file pointer */ } /* END OF FUNCTION COMPUTEFND */ /***************************************************************************** ** ** ** FUNCTION SCREEN ** ** ** ** This function is called from FUNCTION MAIN, and prompts the user for Mx, ** ** My, Mxy, and r. It then computes F, G1, G2, and d, and prints them to ** ** the screen. Unlike the identical computations in the batch mode, all of ** ** the G2 values are printed to the screen. ** ** ** ** Functions called: ** ** getn -- prompts user for Mx, My, Mxy, and r. ** *****************************************************************************/ void screen() { float F = 0; /* the proportion of fragments shared between OTUs X & Y */ float G1 = 0; /* value entering equation 5.54 */ float G2 = 0; /* value leaving equation 5.54 */ float d = 0; /* the number of nucleotide substitutions/base pair */ register int i; /* loop control variable */ getm(); /* get Mx, My, Mxy, and r from the user */ if ((m.mx <= 0) || (m.my <= 0)) { /* if invalid comparison */ printf("Mx and/or My are equal to zero.\r\nUnder this condition, "); printf("the comparison cannot be made.\r\n"); } else if (m.mz <= 0) { /* if no fragments in common */ printf("Mxy is equal to zero.\r\nUnder this condition, "); printf("F is zero, and d becomes infinite.\r\n"); } else { /* everything's okay, compute F and d */ F = ((2*m.mz)/(m.mx + m.my)); /* compute F-hat via eqn 5.54 */ G1 = sqrt(sqrt(F)); /* compute start value */ printf("\r\nF = %f\r\n",F); /* print F-hat */ printf("G1 = %f\r\n",G1); /* print start value */ for (i=0; ((fabs((G2=sqrt(sqrt(F*(3-2*G1)))) - G1)) > (TEST*G2)); ++i) { printf("G2 = %f\r\n",G2); /* compute G2 via eqn 5.54 */ G1 = G2; } d = -(2*log(G2))/m.rv; /* compute d via eqn. 5.55 */ printf("F = %f\r\n",F); /* print F-hat */ printf("d = %f\r\n",d); /* print d-hat */ } } /* END OF FUNCTION SCREEN */ /***************************************************************************** ** ** ** FUNCTION GETM ** ** ** ** This function is called from FUNCTION SCREEN, and prompts the user for ** ** Mx, My, Mxy, and r. It then puts them into the global structure M, and ** ** then returns to FUNCTION SCREEN. ** ** ** ** Functions called: none. ** *****************************************************************************/ void getm() { m.mx = 0; /* initialize structure M */ m.my = 0; m.mz = 0; m.rv = 0; /* get data from user: */ printf("\rEnter m(x): "); /* get Mx */ scanf("%f",&m.mx); printf("\rEnter m(y): "); /* get My */ scanf("%f",&m.my); printf("\rEnter m(xy): "); /* get Mxy */ scanf("%f",&m.mz); printf("\rEnter r-value: "); /* get # of base pairs in restriction site */ scanf("%f",&m.rv); printf("\rMx: %6.1f, My: %6.1f, Mxy: %6.1f, r: %3.1f\n",m.mx,m.my,m.mz,m.rv); } /* END OF FUNCTION GETM */ /****************************************************************************/