| 1 | #include "mltaln.h" |
|---|
| 2 | #include <sys/types.h> |
|---|
| 3 | #include <unistd.h> |
|---|
| 4 | #define DEBUG 0 |
|---|
| 5 | #define TEST 0 |
|---|
| 6 | |
|---|
| 7 | |
|---|
| 8 | int howmanyx( char *s ) |
|---|
| 9 | { |
|---|
| 10 | int val = 0; |
|---|
| 11 | if( scoremtx == -1 ) |
|---|
| 12 | { |
|---|
| 13 | do |
|---|
| 14 | { |
|---|
| 15 | if( !strchr( "atgcuATGCU", *s ) ) val++; |
|---|
| 16 | } while( *++s ); |
|---|
| 17 | } |
|---|
| 18 | else |
|---|
| 19 | { |
|---|
| 20 | do |
|---|
| 21 | { |
|---|
| 22 | if( !strchr( "ARNDCQEGHILKMFPSTWYV", *s ) ) val++; |
|---|
| 23 | } while( *++s ); |
|---|
| 24 | } |
|---|
| 25 | return( val ); |
|---|
| 26 | } |
|---|
| 27 | |
|---|
| 28 | void arguments( int argc, char *argv[] ) |
|---|
| 29 | { |
|---|
| 30 | int c; |
|---|
| 31 | |
|---|
| 32 | inputfile = NULL; |
|---|
| 33 | disopt = 0; |
|---|
| 34 | divpairscore = 0; |
|---|
| 35 | swopt = ""; |
|---|
| 36 | |
|---|
| 37 | while( --argc > 0 && (*++argv)[0] == '-' ) |
|---|
| 38 | { |
|---|
| 39 | while ( (c = *++argv[0]) ) |
|---|
| 40 | { |
|---|
| 41 | switch( c ) |
|---|
| 42 | { |
|---|
| 43 | case 'i': |
|---|
| 44 | inputfile = *++argv; |
|---|
| 45 | fprintf( stderr, "inputfile = %s\n", inputfile ); |
|---|
| 46 | --argc; |
|---|
| 47 | goto nextoption; |
|---|
| 48 | case 'I': |
|---|
| 49 | disopt = 1; |
|---|
| 50 | break; |
|---|
| 51 | case 'A': |
|---|
| 52 | swopt = "-A"; |
|---|
| 53 | break; |
|---|
| 54 | default: |
|---|
| 55 | fprintf( stderr, "illegal option %c\n", c ); |
|---|
| 56 | argc = 0; |
|---|
| 57 | break; |
|---|
| 58 | } |
|---|
| 59 | } |
|---|
| 60 | nextoption: |
|---|
| 61 | ; |
|---|
| 62 | } |
|---|
| 63 | if( argc != 0 ) |
|---|
| 64 | { |
|---|
| 65 | fprintf( stderr, "options: -i\n" ); |
|---|
| 66 | exit( 1 ); |
|---|
| 67 | } |
|---|
| 68 | } |
|---|
| 69 | |
|---|
| 70 | int main( int argc, char *argv[] ) |
|---|
| 71 | { |
|---|
| 72 | int ktuple; |
|---|
| 73 | int i, j; |
|---|
| 74 | FILE *hat2p; |
|---|
| 75 | FILE *hat3p; |
|---|
| 76 | FILE *infp; |
|---|
| 77 | char **seq = NULL; // by D.Mathog |
|---|
| 78 | char **seq1; |
|---|
| 79 | char **name; |
|---|
| 80 | char **name1; |
|---|
| 81 | static int nlen1[M]; |
|---|
| 82 | double **mtx; |
|---|
| 83 | double **mtx2; |
|---|
| 84 | static int nlen[M]; |
|---|
| 85 | static char b[B]; |
|---|
| 86 | double max; |
|---|
| 87 | char com[1000]; |
|---|
| 88 | int opt[M]; |
|---|
| 89 | int res; |
|---|
| 90 | char *home; |
|---|
| 91 | char *fastapath; |
|---|
| 92 | char queryfile[B]; |
|---|
| 93 | char datafile[B]; |
|---|
| 94 | char fastafile[B]; |
|---|
| 95 | char hat2file[B]; |
|---|
| 96 | int pid = (int)getpid(); |
|---|
| 97 | LocalHom **localhomtable, *tmpptr; |
|---|
| 98 | #if 0 |
|---|
| 99 | home = getenv( "HOME" ); |
|---|
| 100 | #else /* $HOME wo tsukau to fasta ni watasu hikisuu ga afureru */ |
|---|
| 101 | home = NULL; |
|---|
| 102 | #endif |
|---|
| 103 | fastapath = getenv( "FASTA_4_MAFFT" ); |
|---|
| 104 | if( !fastapath ) |
|---|
| 105 | fastapath = "fasta34"; |
|---|
| 106 | |
|---|
| 107 | #if DEBUG |
|---|
| 108 | if( home ) fprintf( stderr, "home = %s\n", home ); |
|---|
| 109 | #endif |
|---|
| 110 | if( !home ) home = ""; |
|---|
| 111 | sprintf( queryfile, "%s/tmp/query-%d", home, pid ); |
|---|
| 112 | sprintf( datafile, "%s/tmp/data-%d", home, pid ); |
|---|
| 113 | sprintf( fastafile, "%s/tmp/fasta-%d", home, pid ); |
|---|
| 114 | sprintf( hat2file, "hat2-%d", pid ); |
|---|
| 115 | |
|---|
| 116 | |
|---|
| 117 | arguments( argc, argv ); |
|---|
| 118 | |
|---|
| 119 | if( inputfile ) |
|---|
| 120 | { |
|---|
| 121 | infp = fopen( inputfile, "r" ); |
|---|
| 122 | if( !infp ) |
|---|
| 123 | { |
|---|
| 124 | fprintf( stderr, "Cannot open %s\n", inputfile ); |
|---|
| 125 | exit( 1 ); |
|---|
| 126 | } |
|---|
| 127 | } |
|---|
| 128 | else |
|---|
| 129 | infp = stdin; |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | #if 0 |
|---|
| 134 | PreRead( stdin, &njob, &nlenmax ); |
|---|
| 135 | #else |
|---|
| 136 | dorp = NOTSPECIFIED; |
|---|
| 137 | getnumlen( infp ); |
|---|
| 138 | #endif |
|---|
| 139 | |
|---|
| 140 | if( dorp == 'd' ) |
|---|
| 141 | { |
|---|
| 142 | scoremtx = -1; |
|---|
| 143 | pamN = NOTSPECIFIED; |
|---|
| 144 | } |
|---|
| 145 | else |
|---|
| 146 | { |
|---|
| 147 | nblosum = 62; |
|---|
| 148 | scoremtx = 1; |
|---|
| 149 | } |
|---|
| 150 | constants( njob, seq ); |
|---|
| 151 | |
|---|
| 152 | rewind( infp ); |
|---|
| 153 | |
|---|
| 154 | name = AllocateCharMtx( njob, B+1 ); |
|---|
| 155 | name1 = AllocateCharMtx( njob, B+1 ); |
|---|
| 156 | seq = AllocateCharMtx( njob, nlenmax+1 ); |
|---|
| 157 | seq1 = AllocateCharMtx( 2, nlenmax+1 ); |
|---|
| 158 | mtx = AllocateDoubleMtx( njob, njob ); |
|---|
| 159 | mtx2 = AllocateDoubleMtx( njob, njob ); |
|---|
| 160 | localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) ); |
|---|
| 161 | for( i=0; i<njob; i++) |
|---|
| 162 | { |
|---|
| 163 | localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) ); |
|---|
| 164 | for( j=0; j<njob; j++) |
|---|
| 165 | { |
|---|
| 166 | localhomtable[i][j].start1 = -1; |
|---|
| 167 | localhomtable[i][j].end1 = -1; |
|---|
| 168 | localhomtable[i][j].start2 = -1; |
|---|
| 169 | localhomtable[i][j].end2 = -1; |
|---|
| 170 | localhomtable[i][j].opt = -1.0; |
|---|
| 171 | localhomtable[i][j].next = NULL; |
|---|
| 172 | } |
|---|
| 173 | } |
|---|
| 174 | |
|---|
| 175 | #if 0 |
|---|
| 176 | FRead( infp, name, nlen, seq ); |
|---|
| 177 | #else |
|---|
| 178 | readData_pointer( infp, name, nlen, seq ); |
|---|
| 179 | #endif |
|---|
| 180 | fclose( infp ); |
|---|
| 181 | |
|---|
| 182 | if( scoremtx == -1 ) ktuple = 6; |
|---|
| 183 | else ktuple = 1; |
|---|
| 184 | |
|---|
| 185 | for( i=0; i<njob; i++ ) |
|---|
| 186 | { |
|---|
| 187 | gappick0( seq1[0], seq[i] ); |
|---|
| 188 | strcpy( seq[i], seq1[0] ); |
|---|
| 189 | } |
|---|
| 190 | for( j=0; j<njob; j++ ) |
|---|
| 191 | { |
|---|
| 192 | sprintf( name1[j], "+==========+%d ", j ); |
|---|
| 193 | nlen1[j] = nlen[j]; |
|---|
| 194 | } |
|---|
| 195 | hat2p = fopen( datafile, "w" ); |
|---|
| 196 | if( !hat2p ) ErrorExit( "Cannot open datafile." ); |
|---|
| 197 | WriteForFasta( hat2p, njob, name1, nlen1, seq ); |
|---|
| 198 | fclose( hat2p ); |
|---|
| 199 | |
|---|
| 200 | for( i=0; i<njob; i++ ) |
|---|
| 201 | { |
|---|
| 202 | // fprintf( stderr, "### i = %d\n", i ); |
|---|
| 203 | hat2p = fopen( datafile, "w" ); |
|---|
| 204 | if( !hat2p ) ErrorExit( "Cannot open datafile." ); |
|---|
| 205 | WriteForFasta( hat2p, njob-i, name1+i, nlen1+i, seq+i ); |
|---|
| 206 | fclose( hat2p ); |
|---|
| 207 | |
|---|
| 208 | seq1[0] = seq[i]; |
|---|
| 209 | nlen1[0] = nlen[i]; |
|---|
| 210 | |
|---|
| 211 | hat2p = fopen( queryfile, "w" ); |
|---|
| 212 | if( !hat2p ) ErrorExit( "Cannot open queryfile." ); |
|---|
| 213 | WriteForFasta( hat2p, 1, name1+i, nlen1, seq1 ); |
|---|
| 214 | fclose( hat2p ); |
|---|
| 215 | |
|---|
| 216 | |
|---|
| 217 | if( scoremtx == -1 ) |
|---|
| 218 | sprintf( com, "%s %s -z3 -m10 -n -Q -b%d -E%d -d%d %s %s %d > %s", fastapath, swopt, M, M, M, queryfile, datafile, ktuple, fastafile ); |
|---|
| 219 | else |
|---|
| 220 | sprintf( com, "%s %s -z3 -m10 -Q -b%d -E%d -d%d %s %s %d > %s", fastapath, swopt, M, M, M, queryfile, datafile, ktuple, fastafile ); |
|---|
| 221 | res = system( com ); |
|---|
| 222 | if( res ) ErrorExit( "error in fasta" ); |
|---|
| 223 | |
|---|
| 224 | |
|---|
| 225 | |
|---|
| 226 | hat2p = fopen( fastafile, "r" ); |
|---|
| 227 | if( hat2p == NULL ) |
|---|
| 228 | ErrorExit( "file 'fasta.$$' does not exist\n" ); |
|---|
| 229 | if( scoremtx == -1 ) |
|---|
| 230 | res = ReadFasta34m10_nuc( hat2p, mtx[i], i, name1, localhomtable[i] ); |
|---|
| 231 | else |
|---|
| 232 | res = ReadFasta34m10( hat2p, mtx[i], i, name1, localhomtable[i] ); |
|---|
| 233 | fclose( hat2p ); |
|---|
| 234 | |
|---|
| 235 | if( res < njob - i ) |
|---|
| 236 | { |
|---|
| 237 | fprintf( stderr, "count (fasta34 -z 3) = %d\n", res ); |
|---|
| 238 | exit( 1 ); |
|---|
| 239 | } |
|---|
| 240 | |
|---|
| 241 | |
|---|
| 242 | if( i == 0 ) |
|---|
| 243 | for( j=0; j<njob; j++ ) opt[j] = (int)mtx[0][j]; |
|---|
| 244 | |
|---|
| 245 | |
|---|
| 246 | #if 0 |
|---|
| 247 | { |
|---|
| 248 | int ii, jj; |
|---|
| 249 | if( i < njob-1 ) for( jj=i; jj<i+5; jj++ ) |
|---|
| 250 | fprintf( stdout, "mtx[%d][%d] = %f\n", i+1, jj+1, mtx[i][jj] ); |
|---|
| 251 | } |
|---|
| 252 | #endif |
|---|
| 253 | fprintf( stderr, "query : %4d / %5d\r", i+1, njob ); |
|---|
| 254 | } |
|---|
| 255 | |
|---|
| 256 | for( i=0; i<njob; i++ ) |
|---|
| 257 | { |
|---|
| 258 | max = mtx[i][i]; |
|---|
| 259 | if( max == 0.0 ) |
|---|
| 260 | { |
|---|
| 261 | for( j=0; j<njob; j++ ) |
|---|
| 262 | mtx2[i][j] = 2.0; |
|---|
| 263 | } |
|---|
| 264 | else |
|---|
| 265 | { |
|---|
| 266 | for( j=0; j<njob; j++ ) |
|---|
| 267 | { |
|---|
| 268 | // fprintf( stderr, "##### mtx[%d][%d] = %f\n", i, j, mtx[i][j] ); |
|---|
| 269 | mtx2[i][j] = ( max - mtx[MIN(i,j)][MAX(i,j)] ) / max * 2.0; |
|---|
| 270 | // fprintf( stdout, "max = %f, mtx[%d][%d] = %f -> %f\n", max, i+1, j+1, mtx[i][j], mtx2[i][j] ); |
|---|
| 271 | } |
|---|
| 272 | } |
|---|
| 273 | } |
|---|
| 274 | for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ ) |
|---|
| 275 | { |
|---|
| 276 | // fprintf( stdout, "mtx2[%d][%d] = %f, %f\n", i+1, j+1, mtx2[i][j], mtx2[j][i] ); |
|---|
| 277 | mtx2[i][j] = MIN( mtx2[i][j], mtx2[j][i] ); |
|---|
| 278 | } |
|---|
| 279 | |
|---|
| 280 | #if 0 |
|---|
| 281 | { |
|---|
| 282 | int ii, jj; |
|---|
| 283 | if( i < njob-1 ) for( jj=i+1; jj<njob; jj++ ) |
|---|
| 284 | fprintf( stderr, "mtx2[][] = %f\n", mtx2[i][jj] ); |
|---|
| 285 | } |
|---|
| 286 | #endif |
|---|
| 287 | |
|---|
| 288 | for( i=0; i<njob; i++ ) name[i][0] = '='; |
|---|
| 289 | |
|---|
| 290 | if( disopt ) |
|---|
| 291 | { |
|---|
| 292 | strcpy( b, name[0] ); |
|---|
| 293 | sprintf( name[0], "=query====lgth=%04d-%04d %.*s", nlen[0], howmanyx( seq[0] ), B-30, b ); |
|---|
| 294 | #if 0 |
|---|
| 295 | strins( b, name[0] ); |
|---|
| 296 | #endif |
|---|
| 297 | for( i=1; i<njob; i++ ) |
|---|
| 298 | { |
|---|
| 299 | strcpy( b, name[i] ); |
|---|
| 300 | sprintf( name[i], "=opt=%04d=lgth=%04d-%04d %.*s", opt[i], nlen[i], howmanyx( seq[i] ), B-30, b ); |
|---|
| 301 | #if 0 |
|---|
| 302 | strins( b, name[i] ); |
|---|
| 303 | #endif |
|---|
| 304 | } |
|---|
| 305 | } |
|---|
| 306 | |
|---|
| 307 | hat2p = fopen( hat2file, "w" ); |
|---|
| 308 | if( !hat2p ) ErrorExit( "Cannot open hat2." ); |
|---|
| 309 | WriteHat2_pointer( hat2p, njob, name, mtx2 ); |
|---|
| 310 | fclose( hat2p ); |
|---|
| 311 | |
|---|
| 312 | #if 1 |
|---|
| 313 | fprintf( stderr, "##### writing hat3\n" ); |
|---|
| 314 | hat3p = fopen( "hat3", "w" ); |
|---|
| 315 | if( !hat3p ) ErrorExit( "Cannot open hat3." ); |
|---|
| 316 | for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ ) |
|---|
| 317 | { |
|---|
| 318 | for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next ) |
|---|
| 319 | { |
|---|
| 320 | if( tmpptr->opt == -1.0 ) continue; |
|---|
| 321 | fprintf( hat3p, "%d %d %d %6.3f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next ); |
|---|
| 322 | } |
|---|
| 323 | } |
|---|
| 324 | fclose( hat3p ); |
|---|
| 325 | #endif |
|---|
| 326 | |
|---|
| 327 | sprintf( com, "/bin/rm %s %s %s", queryfile, datafile, fastafile ); |
|---|
| 328 | system( com ); |
|---|
| 329 | |
|---|
| 330 | #if 0 |
|---|
| 331 | sprintf( com, ALNDIR "/supgsdl < %s", hat2file ); |
|---|
| 332 | res = system( com ); |
|---|
| 333 | if( res ) ErrorExit( "error in spgsdl" ); |
|---|
| 334 | #endif |
|---|
| 335 | |
|---|
| 336 | sprintf( com, "mv %s hat2", hat2file ); |
|---|
| 337 | res = system( com ); |
|---|
| 338 | if( res ) ErrorExit( "error in mv" ); |
|---|
| 339 | |
|---|
| 340 | SHOWVERSION; |
|---|
| 341 | exit( 0 ); |
|---|
| 342 | } |
|---|
