#include #include #include #include #include #include #include #include #include #include //#include "/work/bjornw/source/c/pdb/molecule.h" //#include "/work/bjornw/source/c/pdb/quality_measure/lgscore.h" #include "/afs/pdc.kth.se/home/b/bjornw/source/c/pdb/molecule.h" #include "/afs/pdc.kth.se/home/b/bjornw/source/c/pdb/quality_measure/lgscore.h" //#include "/afs/pdc.kth.se/home/b/bjornw/bjorn/source/c/pdb/quality_molecule.h" #define MAXFILE 10000 #define STRING_BUFFER 1000 void usage(); int SameMethod(char *file1, char *file2); int ispdb(char *filename); char *basename(char *fname); main(int argc,char *argv[]) /* Main routine */ { DIR *dip; FILE *fp; struct dirent *dit; char dir[PATH_MAX]=""; char file2[PATH_MAX]; char tempfilename[PATH_MAX]; char infilename[PATH_MAX]=""; char matrix_outfile[PATH_MAX]=""; char filenames[MAXFILE][100]; char filenames_with_path[MAXFILE][PATH_MAX]; //char **filenames='\0'; //[MAXFILE][100]; //char **filenames_with_path='\0';//[MAXFILE][STRING_BUFFER]; int i,j,k; int files=0; int L=4; int step=2; double minsim=121.0; double d0=2.23606798; //sqrt(5); double factor=0.5; char temp[PATH_MAX]; lgscore LG[1]; // double LG_average[MAXFILE]={0}; //double S_average[MAXFILE]={0}; //double Sstr[MAXFILE][MAXRES]={{0}}; // int number_of_comparisons[MAXFILE]={0}; double *LG_average,*S_average; double **Sstr; int *number_of_comparisons; int total_number_of_comparisons=0; int maxlen=0; int userdef_len=-1; double rmsd=0; double sqrt5=2.23606798; //sqrt(5); char target[STRING_BUFFER]="T0XXX"; int verbose=0; int fastmode=0; int memorymode=1; int bytes_freed=0; int target_defined_by_user=0; int caspoutput=0; int lgscoreoutput=0; int force_compare_all=0; int only_pairwise=0; int output_similarity_matrix=0; int maxfile=MAXFILE; //molecule m[MAXFILE]; dyn_molecule *dm; static clock_t st_time; static clock_t en_time; static struct tms st_cpu; static struct tms en_cpu; /* Parse command line for PDB files and options */ //i=1; //printf("%d %d %d %d\n",sizeof(dm),sizeof(dyn_molecule),sizeof(atm),sizeof(double)); //exit(0); i=1; if(argc > 2) { while (id_name); // strcpy(filenames[files],dir); if(dit->d_name[0] != '.') { strcpy(tempfilename,dir); strcat(tempfilename,dit->d_name); if(ispdb(tempfilename)) { // printf("%d %s\n",files,tempfilename); //filenames=(char**)realloc(filenames,sizeof(char)*20); //(files+1)); ////filenames[files]=(char*)realloc(filenames[files],sizeof(char)*(strlen(dit->d_name)+1)); //filenames[files]=(char*)realloc(filenames[files],sizeof(char)*20000); //printf("%s\n","first"); // //filenames_with_path=(char**)realloc(filenames_with_path,sizeof(char)*20); //(files+1)); //// filenames_with_path[files]=(char*)realloc(filenames_with_path[files],sizeof(char)*(strlen(tempfilename)+1)); //filenames_with_path[files]=(char*)realloc(filenames_with_path[files],sizeof(char)*2000); //printf("%s\n","second"); strcpy(filenames[files],dit->d_name); strcpy(filenames_with_path[files],tempfilename); //strcat(filenames_with_path[files],dir); //strcat(filenames_with_path[files],filenames[files]); if(filenames[files][0] == 'T' && filenames[files][1] == '0' && !target_defined_by_user) { strncpy(target,filenames[files],5); target[5]='\0'; } files++; if(files>MAXFILE) { fprintf(stderr,"Maximum number of files %d exceded. Increase MAXFILE.\n",MAXFILE); exit(1); } } } } if (closedir(dip) == -1) { perror("closedir"); return 0; } } if(strcmp(argv[i],"-i")==0) { i++; strcpy(infilename,argv[i]); // printf("%s\n",infilename); fp=fopen(infilename,"r"); if(fp!=NULL) { files=0; //while(fscanf(fp,"%s",tempfilename)!=EOF) while(fscanf(fp,"%s",tempfilename)!=EOF) { //printf("%s %d\n",tempfilename,strlen(tempfilename)); if(ispdb(tempfilename)) { strcpy(filenames_with_path[files],tempfilename); // printf("%s %d\n",filenames_with_path[files],strlen(filenames_with_path[files])); /// printf("%s %d\n",filenames_with_path[files],strlen(filenames_with_path[files])); strcpy(filenames[files],basename(tempfilename)); // printf("%s\n", basename(tempfilename)); // filenames[files]=tempfilename[basename(tempfilename)]; // strcpy(filenames[files],tempfilename); // strcpy(filenames[i],buff); //tmp=ftell(fp); //filenames[files]; // printf("%s %s %d %s\n",filenames_with_path[files],filenames[files],files,basename(tempfilename)); files++; } } //fclose(fp); } else { printf("Cannot open %s",tempfilename); exit(0); } fclose(fp); } if(strcmp(argv[i],"-m")==0) { i++; strcpy(matrix_outfile,argv[i]); output_similarity_matrix=1; force_compare_all=1; } if(strcmp(argv[i],"-L")==0) { i++; userdef_len=atoi(argv[i]); } if(strcmp(argv[i],"-t")==0) { i++; target_defined_by_user=1; strcpy(target,argv[i]); } if(strcmp(argv[i],"-f")==0) { fastmode=1; L=4; minsim=289.0; factor=2; } if(strcmp(argv[i],"-m")==0) { i++; memorymode=atoi(argv[i]); if(memorymode!=0 || memorymode!=1) { fprintf(stderr,"Illegal -m switch only -m 0 or -m 1 is allowed!\n\n"); usage(); } } if(strcmp(argv[i],"-A")==0) { force_compare_all=1; } if(strcmp(argv[i],"-only_pairwise")==0) { only_pairwise=1; } if(strcmp(argv[i],"-casp")==0) { caspoutput=1; } if(strcmp(argv[i],"-lgscore")==0) { lgscoreoutput=1; } if(strcmp(argv[i],"-v")==0) { printf("Verbose mode on start talking...\n"); verbose=1; } if(strcmp(argv[i],"-L0")==0) { i++; L=atoi(argv[i]); } if(strcmp(argv[i],"-step")==0) { i++; step=atoi(argv[i]); } if(strcmp(argv[i],"-factor")==0) { i++; factor=atof(argv[i]); } if(strcmp(argv[i],"-minsim")==0) { i++; minsim=atof(argv[i]); } i++; } // Done readin input arguments: // printf("%s\n",matrix_outfile); st_time=times(&st_cpu); // Read all CA coordinates to memory if(memorymode) { // printf("%d",files); dm=malloc(sizeof(dyn_molecule)*files); //for(i=0;imaxlen) { maxlen=dm[i].atm[dm[i].residues-1].resnum; } //if(dm[i].method!=NULL) // { // printf("%s\n",dm[i].method); // } } else { fprintf(stderr,"error reading %s\n",dm[i].filename); exit(0); } //else // { // fprintf(stderr,"%s\n",m[i].filename); //} } } // exit(0); //int residues=100; //int totfiles=; // //dm=malloc(sizeof(dyn_molecule)*totfiles); //for(i=0;imaxlen) { maxlen=userdef_len; } // printf("%d\n",maxlen); //exit(1); for(i=0;i1) { fp=fopen(matrix_outfile,"w"); if(fp==NULL) { fprintf(stderr,"Couldn't open matrix otufile for writing (%s)\n",matrix_outfile); exit(0); } } for(i=0;imaxlen) // { // maxlen=LG[0].residues; // } // In the case of CASP I can assume that the resnum is the correct index for(k=0;k1) { if(userdef_len >= 0 && userdef_len0) { printf("REMARK Input directory: %s\n",dir); } if(strlen(infilename)>0) { printf("REMARK Input file: %s\n",infilename); } printf("REMARK Total number of pdbs: %d\n",files); printf("REMARK Total number of comparisons: %d\n",total_number_of_comparisons); printf("REMARK Total CPU time used: %8.2lf seconds\n",(double)(((double)en_cpu.tms_utime - (double)st_cpu.tms_utime)/100)); printf("REMARK LGscore parameters L=%d minsim=%lf factor=%lf step=%d\n",L,minsim,factor,step); if(userdef_len>0) { printf("REMARK Target sequence length (user defined): %d\n",userdef_len); if(userdef_len0) { maxlen=userdef_len; } for(i=0;i0.0001) { //Sstr[i][j]/=number_of_comparisons[i]; // Sstr[i][j]/=number_of_comparisons[i]; if(caspoutput==1) { rmsd=d0*sqrt(1/Sstr[i][j]-1); if(rmsd>15) { rmsd=15; } printf(" %4.3lf",rmsd); } else { printf(" %4.3lf",Sstr[i][j]); } } else { printf(" X"); } } printf("\n"); } printf("END\n"); } else { printf("Too few pdbfiles (%d) found\n",files); } //Freeing up some allocated memory bytes_freed=0; for(i=0;i\n"); fprintf(stderr,"\t\t-i \n"); fprintf(stderr,"\t\t-casp will output local rmsd as local quality (default is average local S-score)\n"); fprintf(stderr,"\t\t-lgscore will output average LGscore as global quality measure (default is average S-score)\n"); fprintf(stderr,"\t\t-L \n"); fprintf(stderr,"\t\t-t \n"); fprintf(stderr,"\t\t-m \n"); fprintf(stderr,"\t\t-A \n"); fprintf(stderr,"\t\t-only_pairwise \n"); //fprintf(stderr,"\t\t-f \n"); //fprintf(stderr,"\t\t-m \n"); exit(0); } int SameMethod(char *file1, char *file2) { int ident=0; int diffpos=0; int i=0,j=0; int len1,len2; int dot_counter=0; char method1[100]; char method2[100]; len1=strlen(file1); len2=strlen(file2); if(file1[len1-3]=='T' && file1[len1-2]=='S' && file2[len2-3]=='T' && file2[len2-2]=='S' || file1[len1-7]=='A' && file1[len1-6]=='L' & file1[len1-3]=='p' && file1[len1-2]=='d' && file1[len1-1]=='b' && file2[len2-7]=='A' && file2[len2-6]=='L' & file2[len2-3]=='p' && file2[len2-2]=='d' && file2[len2-1]=='b' ) { //CASP naming convension ending with TS{rank} //OR AL{rank}.pdb if(len1 != len2) { return 0; } //All files here will have the same length //Check for identical matches for(i=0;i0) { return 1; } } return 0; } char *basename(char *fname) { int sl; char *s; // if (strlen(fname) > STRING_BUFFER-4) // error(NIL,"filename is too long",NIL); s = strrchr(fname,'/'); if (s) fname = s+1; /* Process suffix */ return fname; }