/*
@author: joao lopes
@workplace: Reading University
@date: 8th May 2009
*/
#include "interface.h"
/*
It creates an ABC .len file given a Nexus file (.nex)
@arg length format filename
@arg output filename
*/
int createFreqTab4(char *input,char *output){
int cloc, cpop, cdna, csamp, csite, i, j, //iterators
nr1, //retrieve info from input file
nr2, //retrieve info from input file
reachEOF, //check if EOF has been reached
endSets, //check if end of Sets has been reached
endPop, //check if end of pop has been reached
firstInd, //check if using '-' to define populations
equalDna, //auxiliar to build the .len file
outsize, //size of output file name
npop, //number of populations
nsites, //number of bases
ndna, //number of different dna data
nsamp, //total number of samples
*lsamp, //size of sample per pop
*id, //identify the belonging population of each used individual
**freq; //freq of dna data by no_pop by no. of diff dna data
char c1, //gets the value of a char temporarly
info1[MAXCHAR], //retrieve info from the input file
info2[MAXCHAR], //retrieve info from the input file
aux[MAXCHAR], //auxiliar
*outname, //name of the output file
*mainseq, //main sequence of a loci
**lpop=NULL, //list of the populations names
**indiv, //list of samples
**valS, //list of all the diff dna sequencies
***genotS; //list of all the individual's haplotype per pop per loci
FILE *inp, //pntr to the input file
*outp; //pntr to the output file
time_t startClock; //time when the program starts
const struct tm *startTime; //struct time when the program starts
inp = fopen(input,"r");
if(inp == NULL)
return 1; //cannot open .nex file
outsize = strlen(output) + 5;
outname = (char *)malloc(outsize*sizeof(char));
strcpy(outname,output);
outp = fopen(strcat(outname,".len"),"w");
if(outp == NULL)
return 2; //cannot create output file
time( &startClock ); // Get time in seconds
startTime = localtime( &startClock ); // Convert time to struct tm form
fgets(aux,MAXCHAR,inp); //#NEXUS
reachEOF=0;
while(!reachEOF){
//reads input file
while(isspace(c1=getc(inp))||isendline(c1)||c1=='\t');
if(c1 == '[')
while((c1 = getc(inp))!=']');
else if(c1==EOF)
reachEOF=1;
else
ungetc(c1,inp);
if(!reachEOF){
while(isspace(c1=getc(inp))||isendline(c1)||c1=='\t');
ungetc(c1,inp);
fscanf(inp,"%s",aux);
if(strcmp(aux,"begin")==0||strcmp(aux,"BEGIN")==0||strcmp(aux,"Begin")==0){
fscanf(inp,"%s",aux);
if(strcmp(aux,"data;")==0||strcmp(aux,"DATA;")==0||strcmp(aux,"Data;")==0){
fscanf(inp,"%s",aux); //DIMENSIONS
fscanf(inp,"%s",aux);
i=0; j=0;
while(!isdigit(aux[i]))
i++;
while(isdigit(aux[i+j])){
info1[j]=aux[i+j];
j++;
}
nsamp=atoi(info1);
fscanf(inp,"%s",aux);
i=0; j=0;
while(!isdigit(aux[i]))
i++;
while(isdigit(aux[i+j])){
info1[j]=aux[i+j];
j++;
}
nsites=atoi(info1);
fscanf(inp,"%s",aux); //FORMAT
fscanf(inp,"%s",aux);
i=0; j=0;
while(aux[i]!='=')
i++;
while(isalpha(aux[1+i+j])){
info1[j]=aux[1+i+j];
j++;
}
info1[j]='\0';
if(!(strcmp(info1,"DNA")==0||strcmp(info1,"dna")==0||strcmp(info1,"Dna")==0))
return 3; //program can't deal with specified data
//MATRIX: go through all the individuals
fscanf(inp,"%s",aux); //MATRIX
indiv=(char**)malloc(nsamp*sizeof(char*));
for(csamp=0; csamp<nsamp; csamp++)
indiv[csamp]=(char*)malloc((nsites+1)*sizeof(char));
csamp=0;
while(csamp<nsamp){
while(isspace(c1=getc(inp))||isendline(c1)||c1=='\t');
if(c1 == '[')
while((c1 = getc(inp))!=']');
else
ungetc(c1,inp);
while(isspace(c1=getc(inp))||isendline(c1)||c1=='\t');
if(c1 == '[')
ungetc(c1,inp);
else{
ungetc(c1,inp);
fscanf(inp,"%s",aux);
while(isspace(c1=getc(inp))||c1=='\t');
ungetc(c1,inp);
csite=0;
while(csite<nsites){
if((c1=getc(inp))=='T'||c1=='A'||c1=='C'||c1=='G'||
c1=='t'||c1=='a'||c1=='c'||c1=='g'||c1=='N'||c1=='-'){
indiv[csamp][csite]=c1;
csite++;
}
else if(isspace(c1)||isendline(c1)||c1=='\t'){
while(isspace(c1=getc(inp))||c1=='\t');
ungetc(c1,inp);
}
else
return 4; //problem reading MATRIX
}
csamp++;
}
}
while(isspace(c1=getc(inp))||isendline(c1)||c1=='\t');
ungetc(c1,inp);
fscanf(inp,"%s",aux); //;
fscanf(inp,"%s",aux); //END;
}
else if(strcmp(aux,"sets;")==0||strcmp(aux,"SETS;")==0||strcmp(aux,"Sets;")==0){
id=(int*)malloc(nsamp*sizeof(int));
for(csamp=0; csamp<nsamp; csamp++)
id[csamp]=-1;
//SETS: go through all the populations
cpop=0;
endSets=0;
while(!endSets){
fscanf(inp,"%s",aux);
if(strcmp(aux,"end;")==0||strcmp(aux,"END;")==0||strcmp(aux,"End;")==0){
endSets=1;
}
else if(strcmp(aux,"TAXSET")==0||strcmp(aux,"taxset")==0||strcmp(aux,"Taxset")==0){
lpop=(char**)myAlloc(lpop,(cpop+1)*sizeof(char*));
lpop[cpop]=(char*)malloc(MAXCHAR*sizeof(char));
fscanf(inp,"%s",aux);
strcpy(lpop[cpop],aux);
fscanf(inp,"%s",aux); //=
//get individuals belonging to pop
endPop=0;
while(!endPop){
fscanf(inp,"%s",aux);
i=0;
firstInd=1;
while(i<strlen(aux)){
if(aux[i]==';'){
endPop=1;
i++;
}
else if(aux[i]=='-'){
firstInd=0;
i++;
}
else if(isdigit(aux[i])&&firstInd){
j=0;
while(isdigit(aux[i])){
info1[j]=aux[i];
i++;
j++;
}
info1[j]='\0';
nr1=atoi(info1);
}
else if(isdigit(aux[i])&&(!firstInd)){
j=0;
while(isdigit(aux[i])){
info2[j]=aux[i];
i++;
j++;
}
info2[j]='\0';
nr2=atoi(info2);
}
else
return 5; //problem reading TAXSET
}
if(firstInd)
id[nr1-1]=cpop;
else{
if(nr1>=nr2)
return 6; //problem defining TAXSET: n1>=n2
for(i=nr1;i<=nr2;i++)
id[i-1]=cpop;
}
}
cpop++;
}
else
return 7; //unrecognized specifier after BEGIN SETS (use TAXSET)
}
npop=cpop;
}
else
return 8; //unrecognized specifier after BEGIN (use DATA or SETS)
}
else if(aux[0]=='['){
for(i=strlen(aux)-1; i>=0; i--)
ungetc(aux[i],inp);
}
else
return 9; //unrecognized command (only BEGIN accepted)
}
}
//fill lsamp and genotS
genotS = (char***)malloc(npop*sizeof(char**));
lsamp = (int *)malloc(npop*sizeof(int));
for(cpop=0; cpop<npop; cpop++){
lsamp[cpop]=0;
genotS[cpop] = (char**)malloc(nsamp*sizeof(char*));
for(csamp=0; csamp<nsamp; csamp++){
genotS[cpop][csamp] = (char *)malloc(nsites*sizeof(char));
}
}
for(csamp=0; csamp<nsamp; csamp++){
if(id[csamp]!=-1){
for(csite=0; csite<nsites; csite++){
genotS[id[csamp]][lsamp[id[csamp]]][csite] = indiv[csamp][csite];
}
lsamp[id[csamp]]++;
}
}
//allocate memory to freq, IDsort, valS and ldna
freq = (int **)malloc(npop*sizeof(int *));
valS = (char **)malloc(nsamp*sizeof(char *));
for(cpop=0;cpop<npop;++cpop){
freq[cpop] = (int *)malloc(nsamp*sizeof(int));
}
//fill freq[cpop][cdna] and valS[cdna][csite]
ndna=0;
for(cpop=0;cpop<npop;cpop++){
//run through every individual
for(csamp=0;csamp<lsamp[cpop];++csamp){
//check if the current haplotype is already listed
for(cdna=0;cdna<ndna;++cdna){
equalDna = 1;
for(csite=0 ; csite<nsites ; csite++){
if(genotS[cpop][csamp][csite] != valS[cdna][csite]){
if(valS[cdna][csite]=='N'||valS[cdna][csite]=='-'){
valS[cdna][csite] = genotS[cpop][csamp][csite];
}
else if(genotS[cpop][csamp][csite]=='A'||genotS[cpop][csamp][csite]=='a'||
genotS[cpop][csamp][csite]=='T'||genotS[cpop][csamp][csite]=='t'||
genotS[cpop][csamp][csite]=='G'||genotS[cpop][csamp][csite]=='g'||
genotS[cpop][csamp][csite]=='C'||genotS[cpop][csamp][csite]=='c'){
equalDna = 0;
break;
}
}
}
if(equalDna){
++freq[cpop][cdna]; //increase freq of current haplotype
break;
}
}
//current haplotype is not listed yet
if(cdna>=ndna){
for(j=0;j<npop;++j){
freq[j][ndna] = 0;//initiate as 0 the freq of current haplotype in other populations
}
freq[cpop][ndna] = 1; //initiate as 1 the freq of current haplotype in its population
valS[cdna] = (char *)malloc((nsites+1)*sizeof(char));
for(csite=0; csite<nsites ; csite++){
valS[cdna][csite] = genotS[cpop][csamp][csite]; //adds current haplotype to valS[][][]
}
valS[cdna][csite] = '\0';
++ndna; //increase by one the number of different haplotypes of a particular locus
}
}
}
//transforms dna sequences into an array of segregating sites
mainseq = (char *)malloc(nsites*sizeof(char));
for(csite=0; csite<nsites; csite++)
mainseq[csite] = valS[0][csite];
for(cdna=0;cdna<ndna;++cdna){
for(csite=0;csite<nsites;++csite){
if(valS[cdna][csite]!=mainseq[csite]&&(valS[cdna][csite]=='A'||
valS[cdna][csite]=='a'||
valS[cdna][csite]=='G'||
valS[cdna][csite]=='g'||
valS[cdna][csite]=='T'||
valS[cdna][csite]=='t'||
valS[cdna][csite]=='C'||
valS[cdna][csite]=='c')){
valS[cdna][csite] = '1';
}
else{
valS[cdna][csite] = '0';
}
}
}
//writes the output file
fprintf(outp,"# Nexus file converted to .len file with nexus2table1.0\n");
fprintf(outp,"# input file: %s\n",input);
fprintf(outp,"# output file: %s.len\n",output);
fprintf(outp,"# date: %s#\n",asctime(startTime));
for(cpop=0;cpop<npop;++cpop){
fprintf(outp,"#Population %d - %s\n",cpop+1,lpop[cpop]);
}
fprintf(outp,"#Locus 1 - locus1\n");
fprintf(outp,"\n%d\n",npop); //outp: npop
fprintf(outp,"1\ns\n\n"); //outp: nloc, ltype
fprintf(outp,"%d\n",ndna); //outp: ndna
for(cpop=0;cpop<npop;++cpop){
for(cdna=0;cdna<ndna;++cdna){
fprintf(outp,"%4d ",freq[cpop][cdna]); //(S)outp: freq[][]
}
fprintf(outp,"\n");
}
fprintf(outp,"\n");
for(cdna=0;cdna<ndna;++cdna)
fprintf(outp,"%4d ",cdna); //(S)outp: cdna
fprintf(outp,"\n\n");
fprintf(outp,"%d \n",nsites); //(S)outp: nsites
for(cdna=0;cdna<ndna;++cdna)
fprintf(outp,"%d %s\n",cdna,valS[cdna]); //(S)outp: cdna, valS[]
fprintf(outp,"\n");
//free stuff
for(csamp=0 ; csamp<nsamp ; csamp++){
free(indiv[csamp]);
}
free(indiv);
for(cpop=0 ; cpop<npop ; cpop++){
for(csamp=0 ; csamp<nsamp ; csamp++){
free(genotS[cpop][csamp]);
}
free(genotS[cpop]);
}
free(genotS);
for(cpop=0; cpop<npop ; cpop++){
free(freq[cpop]);
free(lpop[cpop]);
}
free(freq);
free(lpop);
for(cdna=0 ; cdna<ndna; cdna++)
free(valS[cdna]);
free(valS);
free(mainseq);
free(lsamp);
free(outname);
free(id);
//close files
fclose(inp);
fclose(outp);
return 0;
} //end of main
