#define DEBUG 0
#define IODEBUG 0
#define SCOREOUT 0
+#define SHISHAGONYU 0 // for debug
#define NODIST -9999
static int stdout_dist;
static int store_localhom;
static int store_dist;
+static int nadd;
+static int laste;
+static int lastm;
+static int lastsubopt;
+static int lastonce;
+static int usenaivescoreinsteadofalignmentscore;
+static int specifictarget;
+
+typedef struct _lastres
+{
+ int score;
+ int start1;
+ int start2;
+ char *aln1;
+ char *aln2;
+} Lastres;
+typedef struct _reg
+{
+ int start;
+ int end;
+} Reg;
+
+typedef struct _aln
+{
+ int nreg;
+ Reg *reg1;
+ Reg *reg2;
+} Aln;
+
+typedef struct _lastresx
+{
+ int score;
+ int naln;
+ Aln *aln;
+} Lastresx;
+
+#ifdef enablemultithread
typedef struct _jobtable
{
int i;
int j;
} Jobtable;
-#ifdef enablemultithread
typedef struct _thread_arg
{
int thread_no;
Jobtable *jobpospt;
char **name;
char **seq;
+ char **dseq;
+ int *thereisxineachseq;
LocalHom **localhomtable;
double **distancemtx;
double *selfscore;
char ***bpp;
+ Lastresx **lastresx;
int alloclen;
+ int *targetmap;
pthread_mutex_t *mutex_counter;
pthread_mutex_t *mutex_stdout;
} thread_arg_t;
#endif
+typedef struct _lastcallthread_arg
+{
+ int nq, nd;
+ char **dseq;
+ char **qseq;
+ Lastresx **lastresx;
+#ifdef enablemultithread
+ int thread_no;
+ int *kshare;
+ pthread_mutex_t *mutex;
+#endif
+} lastcallthread_arg_t;
+
static void t2u( char *seq )
{
while( *seq )
}
}
-static float recallpairfoldalign( char **mseq1, char **mseq2, int m1, int m2, int *of1pt, int *of2pt, int alloclen )
+static int removex( char *d, char *m )
+{
+ int val = 0;
+ while( *m != 0 )
+ {
+ if( *m == 'X' || *m == 'x' )
+ {
+ m++;
+ val++;
+ }
+ else
+ {
+ *d++ = *m++;
+ }
+ }
+ *d = 0;
+ return( val );
+}
+
+static void putlocalhom_last( char *s1, char *s2, LocalHom *localhompt, Lastresx *lastresx, char korh )
+{
+ char *pt1, *pt2;
+ int naln, nreg;
+ int iscore;
+ int isumscore;
+ int sumoverlap;
+ LocalHom *tmppt = localhompt;
+ LocalHom *tmppt2;
+ LocalHom *localhompt0;
+ Reg *rpt1, *rpt2;
+ Aln *apt;
+ int nlocalhom = 0;
+ int len;
+
+// fprintf( stderr, "s1=%s\n", s1 );
+// fprintf( stderr, "s2=%s\n", s2 );
+
+
+ naln = lastresx->naln;
+ apt = lastresx->aln;
+
+ if( naln == 0 ) return;
+ while( naln-- )
+ {
+ rpt1 = apt->reg1;
+ rpt2 = apt->reg2;
+ nreg = apt->nreg;
+ isumscore = 0;
+ sumoverlap = 0;
+ while( nreg-- )
+ {
+ if( nlocalhom++ > 0 )
+ {
+// fprintf( stderr, "reallocating ...\n" );
+ tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
+// fprintf( stderr, "done\n" );
+ tmppt = tmppt->next;
+ tmppt->next = NULL;
+ }
+ tmppt->start1 = rpt1->start;
+ tmppt->start2 = rpt2->start;
+ tmppt->end1 = rpt1->end;
+ tmppt->end2 = rpt2->end;
+ tmppt->korh = 'h';
+ if( rpt1 == apt->reg1 ) localhompt0 = tmppt; // ?
+
+// fprintf( stderr, "in putlocalhom, reg1: %d-%d (nreg=%d)\n", rpt1->start, rpt1->end, lastresx->nreg );
+// fprintf( stderr, "in putlocalhom, reg2: %d-%d (nreg=%d)\n", rpt2->start, rpt2->end, lastresx->nreg );
+
+ len = tmppt->end1 - tmppt->start1 + 1;
+
+// fprintf( stderr, "tmppt->start1=%d\n", tmppt->start1 );
+// fprintf( stderr, "tmppt->start2=%d\n", tmppt->start2 );
+
+// fprintf( stderr, "s1+tmppt->start1=%*.*s\n", len, len, s1+tmppt->start1 );
+// fprintf( stderr, "s2+tmppt->start2=%*.*s\n", len, len, s2+tmppt->start2 );
+
+ pt1 = s1 + tmppt->start1;
+ pt2 = s2 + tmppt->start2;
+ iscore = 0;
+ while( len-- )
+ {
+ iscore += n_dis[(int)amino_n[(unsigned char)*pt1++]][(int)amino_n[(unsigned char)*pt2++]]; // - offset \e$B$O$$$i$J$$$+$b\e(B
+// fprintf( stderr, "len=%d, %c-%c, iscore(0) = %d\n", len, *(pt1-1), *(pt2-1), iscore );
+ }
+
+ if( divpairscore )
+ {
+ tmppt->overlapaa = tmppt->end2-tmppt->start2+1;
+ tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
+ }
+ else
+ {
+ isumscore += iscore;
+ sumoverlap += tmppt->end2-tmppt->start2+1;
+ }
+ rpt1++;
+ rpt2++;
+ }
+#if 0
+ fprintf( stderr, "iscore (1)= %d\n", iscore );
+ fprintf( stderr, "al1: %d - %d\n", start1, end1 );
+ fprintf( stderr, "al2: %d - %d\n", start2, end2 );
+#endif
+
+ if( !divpairscore )
+ {
+ for( tmppt2=localhompt0; tmppt2; tmppt2=tmppt2->next )
+ {
+ tmppt2->overlapaa = sumoverlap;
+ tmppt2->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
+// fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
+ }
+ }
+ apt++;
+ }
+}
+
+static int countcomma( char *s )
+{
+ int v = 0;
+ while( *s ) if( *s++ == ',' ) v++;
+ return( v );
+}
+
+static double recallpairfoldalign( char **mseq1, char **mseq2, int m1, int m2, int *of1pt, int *of2pt, int alloclen )
+{
+ static FILE *fp = NULL;
+ double value;
+ char *aln1;
+ char *aln2;
+ int of1tmp, of2tmp;
+
+ if( fp == NULL )
+ {
+ fp = fopen( "_foldalignout", "r" );
+ if( fp == NULL )
+ {
+ fprintf( stderr, "Cannot open _foldalignout\n" );
+ exit( 1 );
+ }
+ }
+
+ aln1 = calloc( alloclen, sizeof( char ) );
+ aln2 = calloc( alloclen, sizeof( char ) );
+
+ readpairfoldalign( fp, *mseq1, *mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
+
+ if( strstr( foldalignopt, "-global") )
+ {
+ fprintf( stderr, "Calling G__align11\n" );
+ value = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
+ *of1pt = 0;
+ *of2pt = 0;
+ }
+ else
+ {
+ fprintf( stderr, "Calling L__align11\n" );
+ value = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, of1pt, of2pt );
+ }
+
+// value = (double)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
+
+ if( aln1[0] == 0 )
+ {
+ fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d. Sequence alignment is used instead.\n", m1+1, m2+1 );
+ }
+ else
+ {
+ strcpy( *mseq1, aln1 );
+ strcpy( *mseq2, aln2 );
+ *of1pt = of1tmp;
+ *of2pt = of2tmp;
+ }
+
+// value = naivepairscore11( *mseq1, *mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
+
+// fclose( fp ); // saigo dake yatta houga yoi.
+
+// fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
+// fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
+
+
+ free( aln1 );
+ free( aln2 );
+
+ return( value );
+}
+
+static void block2reg( char *block, Reg *reg1, Reg *reg2, int start1, int start2 )
+{
+ Reg *rpt1, *rpt2;
+ char *tpt, *npt;
+ int pos1, pos2;
+ int len, glen1, glen2;
+ pos1 = start1;
+ pos2 = start2;
+ rpt1 = reg1;
+ rpt2 = reg2;
+ while( block )
+ {
+ block++;
+// fprintf( stderr, "block = %s\n", block );
+ tpt = strchr( block, ':' );
+ npt = strchr( block, ',' );
+ if( !tpt || tpt > npt )
+ {
+ len = atoi( block );
+ reg1->start = pos1;
+ reg2->start = pos2;
+ pos1 += len - 1;
+ pos2 += len - 1;
+ reg1->end = pos1;
+ reg2->end = pos2;
+// fprintf( stderr, "in loop reg1: %d-%d\n", reg1->start, reg1->end );
+// fprintf( stderr, "in loop reg2: %d-%d\n", reg2->start, reg2->end );
+ reg1++;
+ reg2++;
+ }
+ else
+ {
+ sscanf( block, "%d:%d", &glen1, &glen2 );
+ pos1 += glen1 + 1;
+ pos2 += glen2 + 1;
+ }
+ block = npt;
+
+ }
+ reg1->start = reg1->end = reg2->start = reg2->end = -1;
+
+ while( rpt1->start != -1 )
+ {
+// fprintf( stderr, "reg1: %d-%d\n", rpt1->start, rpt1->end );
+// fprintf( stderr, "reg2: %d-%d\n", rpt2->start, rpt2->end );
+ rpt1++;
+ rpt2++;
+ }
+// *apt1 = *apt2 = 0;
+// fprintf( stderr, "aln1 = %s\n", aln1 );
+// fprintf( stderr, "aln2 = %s\n", aln2 );
+}
+
+
+static void readlastresx_singleq( FILE *fp, int n1, int nameq, Lastresx **lastresx )
+{
+ char *gett;
+ Aln *tmpaln;
+ int prevnaln, naln, nreg;
+#if 0
+ int i, pstart, pend, end1, end2;
+#endif
+ int score, name1, start1, alnSize1, seqSize1;
+ int name2, start2, alnSize2, seqSize2;
+ char strand1, strand2;
+ int includeintoscore;
+ gett = calloc( 10000, sizeof( char ) );
+
+// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
+// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
+
+ while( 1 )
+ {
+ fgets( gett, 9999, fp );
+ if( feof( fp ) ) break;
+ if( gett[0] == '#' ) continue;
+// fprintf( stdout, "gett = %s\n", gett );
+ if( gett[strlen(gett)-1] != '\n' )
+ {
+ fprintf( stderr, "Too long line?\n" );
+ exit( 1 );
+ }
+
+ sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
+ &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
+ &name2, &start2, &alnSize2, &strand2, &seqSize2 );
+
+ if( alg == 'R' && name2 <= name1 ) continue;
+ if( name2 != nameq )
+ {
+ fprintf( stderr, "BUG!!!\n" );
+ exit( 1 );
+ }
+
+// if( lastresx[name1][name2].score ) continue; // dame!!!!
+
+
+ prevnaln = lastresx[name1][name2].naln;
+#if 0
+ for( i=0; i<prevnaln; i++ )
+ {
+ nreg = lastresx[name1][name2].aln[i].nreg;
+
+ pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
+ end1 = start1 + alnSize1;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start1 && start1 <= pend && pend - start1 > 1 ) break;
+ if( pstart <= end1 && end1 <= pend && end1 - pstart > 1 ) break;
+
+ pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
+ end2 = start2 + alnSize2;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start2 && start2 <= pend && pend - start2 > 1 ) break;
+ if( pstart <= end2 && end2 <= pend && end2 - pstart > 1 ) break;
+ }
+ includeintoscore = ( i == prevnaln );
+#else
+ if( prevnaln ) includeintoscore = 0;
+ else includeintoscore = 1;
+#endif
+ if( !includeintoscore && !lastsubopt )
+ continue;
+
+ naln = prevnaln + 1;
+ lastresx[name1][name2].naln = naln;
+// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
+
+ if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
+ exit( 1 );
+ }
+ else
+ lastresx[name1][name2].aln = tmpaln;
+
+ nreg = countcomma( gett )/2 + 1;
+ lastresx[name1][name2].aln[prevnaln].nreg = nreg;
+// lastresx[name1][name2].aln[naln].nreg = -1;
+// lastresx[name1][name2].aln[naln].reg1 = NULL;
+// lastresx[name1][name2].aln[naln].reg2 = NULL;
+// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
+// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
+ block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
+
+ if( includeintoscore )
+ {
+ if( lastresx[name1][name2].score ) score += penalty;
+ lastresx[name1][name2].score += score;
+ }
+
+// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
+ }
+ free( gett );
+}
+
+#ifdef enablemultithread
+#if 0
+static void readlastresx_group( FILE *fp, Lastresx **lastresx )
+{
+ char *gett;
+ Aln *tmpaln;
+ int prevnaln, naln, nreg;
+#if 0
+ int i, pstart, pend, end1, end2;
+#endif
+ int score, name1, start1, alnSize1, seqSize1;
+ int name2, start2, alnSize2, seqSize2;
+ char strand1, strand2;
+ int includeintoscore;
+ gett = calloc( 10000, sizeof( char ) );
+
+// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
+// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
+
+ while( 1 )
+ {
+ fgets( gett, 9999, fp );
+ if( feof( fp ) ) break;
+ if( gett[0] == '#' ) continue;
+// fprintf( stdout, "gett = %s\n", gett );
+ if( gett[strlen(gett)-1] != '\n' )
+ {
+ fprintf( stderr, "Too long line?\n" );
+ exit( 1 );
+ }
+
+ sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
+ &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
+ &name2, &start2, &alnSize2, &strand2, &seqSize2 );
+
+ if( alg == 'R' && name2 <= name1 ) continue;
+
+// if( lastresx[name1][name2].score ) continue; // dame!!!!
+
+ prevnaln = lastresx[name1][name2].naln;
+#if 0
+ for( i=0; i<prevnaln; i++ )
+ {
+ nreg = lastresx[name1][name2].aln[i].nreg;
+
+ pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
+ end1 = start1 + alnSize1;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
+ if( pstart <= end1 && end1 <= pend && end1 - pstart > 3 ) break;
+
+ pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
+ end2 = start2 + alnSize2;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
+ if( pstart <= end2 && end2 <= pend && end2 - pstart > 3 ) break;
+ }
+ includeintoscore = ( i == prevnaln );
+#else
+ if( prevnaln ) includeintoscore = 0;
+ else includeintoscore = 1;
+#endif
+ if( !includeintoscore && !lastsubopt )
+ continue;
+
+ naln = prevnaln + 1;
+ lastresx[name1][name2].naln = naln;
+// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
+
+ if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
+ exit( 1 );
+ }
+ else
+ lastresx[name1][name2].aln = tmpaln;
+
+
+
+ nreg = countcomma( gett )/2 + 1;
+ lastresx[name1][name2].aln[prevnaln].nreg = nreg;
+// lastresx[name1][name2].aln[naln].nreg = -1;
+// lastresx[name1][name2].aln[naln].reg1 = NULL;
+// lastresx[name1][name2].aln[naln].reg2 = NULL;
+// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
+// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
+ block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
+
+ if( includeintoscore )
+ {
+ if( lastresx[name1][name2].score ) score += penalty;
+ lastresx[name1][name2].score += score;
+ }
+
+// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
+ }
+ free( gett );
+}
+#endif
+#endif
+
+static void readlastresx( FILE *fp, int n1, int n2, Lastresx **lastresx, char **seq1, char **seq2 )
+{
+ char *gett;
+ Aln *tmpaln;
+ int prevnaln, naln, nreg;
+#if 0
+ int i, pstart, pend, end1, end2;
+#endif
+ int score, name1, start1, alnSize1, seqSize1;
+ int name2, start2, alnSize2, seqSize2;
+ char strand1, strand2;
+ int includeintoscore;
+ gett = calloc( 10000, sizeof( char ) );
+
+// fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
+// fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
+
+ while( 1 )
+ {
+ fgets( gett, 9999, fp );
+ if( feof( fp ) ) break;
+ if( gett[0] == '#' ) continue;
+// fprintf( stdout, "gett = %s\n", gett );
+ if( gett[strlen(gett)-1] != '\n' )
+ {
+ fprintf( stderr, "Too long line?\n" );
+ exit( 1 );
+ }
+
+ sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
+ &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
+ &name2, &start2, &alnSize2, &strand2, &seqSize2 );
+
+ if( alg == 'R' && name2 <= name1 ) continue;
+
+// if( lastresx[name1][name2].score ) continue; // dame!!!!
+
+ prevnaln = lastresx[name1][name2].naln;
+#if 0
+ for( i=0; i<prevnaln; i++ )
+ {
+ nreg = lastresx[name1][name2].aln[i].nreg;
+
+ pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
+ end1 = start1 + alnSize1;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
+ if( pstart <= end1 && end1 <= pend && end1 - pstart > 3 ) break;
+
+ pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
+ pend = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
+ end2 = start2 + alnSize2;
+// fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
+ if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
+ if( pstart <= end2 && end2 <= pend && end2 - pstart > 3 ) break;
+ }
+ includeintoscore = ( i == prevnaln );
+#else
+ if( prevnaln ) includeintoscore = 0;
+ else includeintoscore = 1;
+#endif
+ if( !includeintoscore && !lastsubopt )
+ continue;
+
+ naln = prevnaln + 1;
+ lastresx[name1][name2].naln = naln;
+// fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
+
+ if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
+ exit( 1 );
+ }
+ else
+ lastresx[name1][name2].aln = tmpaln;
+
+
+
+ nreg = countcomma( gett )/2 + 1;
+ lastresx[name1][name2].aln[prevnaln].nreg = nreg;
+// lastresx[name1][name2].aln[naln].nreg = -1;
+// lastresx[name1][name2].aln[naln].reg1 = NULL;
+// lastresx[name1][name2].aln[naln].reg2 = NULL;
+// fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+ if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
+ {
+ fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
+ exit( 1 );
+ }
+
+// lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
+// lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
+ block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
+
+ if( includeintoscore )
+ {
+ if( lastresx[name1][name2].score ) score += penalty;
+ lastresx[name1][name2].score += score;
+ }
+
+// fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
+ }
+ free( gett );
+}
+
+#ifdef enablemultithread
+#if 0
+static void *lastcallthread_group( void *arg )
+{
+ lastcallthread_arg_t *targ = (lastcallthread_arg_t *)arg;
+ int k, i;
+ int nq = targ->nq;
+ int nd = targ->nd;
+#ifdef enablemultithread
+ int thread_no = targ->thread_no;
+ int *kshare = targ->kshare;
+#endif
+ Lastresx **lastresx = targ->lastresx;
+ char **dseq = targ->dseq;
+ char **qseq = targ->qseq;
+ char command[5000];
+ FILE *lfp;
+ int msize;
+ int klim;
+ int qstart, qend, shou, amari;
+ char kd[1000];
+
+ if( nthread )
+ {
+ shou = nq / nthread;
+ amari = nq - shou * nthread;
+ fprintf( stderr, "shou: %d, amari: %d\n", shou, amari );
+
+ qstart = thread_no * shou;
+ if( thread_no - 1 < amari ) qstart += thread_no;
+ else qstart += amari;
+
+ qend = qstart + shou - 1;
+ if( thread_no < amari ) qend += 1;
+ fprintf( stderr, "%d: %d-%d\n", thread_no, qstart, qend );
+ }
+ k = -1;
+ while( 1 )
+ {
+ if( nthread )
+ {
+ if( qstart > qend ) break;
+ if( k == thread_no ) break;
+ fprintf( stderr, "\n%d-%d / %d (thread %d) \n", qstart, qend, nq, thread_no );
+ k = thread_no;
+ }
+ else
+ {
+ k++;
+ if( k == nq ) break;
+ fprintf( stderr, "\r%d / %d \r", k, nq );
+ }
+
+ if( alg == 'R' ) // if 'r' -> calllast_fast
+ {
+ fprintf( stderr, "Not supported\n" );
+ exit( 1 );
+ }
+ else // 'r'
+ {
+ kd[0] = 0;
+ }
+
+ sprintf( command, "_q%d", k );
+ lfp = fopen( command, "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open %s", command );
+ exit( 1 );
+ }
+ for( i=qstart; i<=qend; i++ )
+ fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
+ fclose( lfp );
+
+// if( alg == 'R' ) msize = MAX(10,k+nq);
+// else msize = MAX(10,nd+nq);
+ if( alg == 'R' ) msize = MAX(10,k*lastm);
+ else msize = MAX(10,nd*lastm);
+
+// fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
+// sprintf( command, "grep '>' _db%sd", kd );
+// system( command );
+ sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
+ if( system( command ) ) exit( 1 );
+
+ sprintf( command, "_lastres%d", k );
+ lfp = fopen( command, "r" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot read _lastres%d", k );
+ exit( 1 );
+ }
+// readlastres( lfp, nd, nq, lastres, dseq, qseq );
+// fprintf( stderr, "Reading lastres\n" );
+ readlastresx_group( lfp, lastresx );
+ fclose( lfp );
+ }
+ return( NULL );
+}
+#endif
+#endif
+
+static void *lastcallthread( void *arg )
+{
+ lastcallthread_arg_t *targ = (lastcallthread_arg_t *)arg;
+ int k, i;
+ int nq = targ->nq;
+ int nd = targ->nd;
+#ifdef enablemultithread
+ int thread_no = targ->thread_no;
+ int *kshare = targ->kshare;
+#endif
+ Lastresx **lastresx = targ->lastresx;
+ char **dseq = targ->dseq;
+ char **qseq = targ->qseq;
+ char command[5000];
+ FILE *lfp;
+ int msize;
+ int klim;
+ char kd[1000];
+
+ k = -1;
+ while( 1 )
+ {
+
+#ifdef enablemultithread
+ if( nthread )
+ {
+ pthread_mutex_lock( targ->mutex );
+ k = *kshare;
+ if( k == nq )
+ {
+ pthread_mutex_unlock( targ->mutex );
+ break;
+ }
+ fprintf( stderr, "\r%d / %d (thread %d) \r", k, nq, thread_no );
+ ++(*kshare);
+ pthread_mutex_unlock( targ->mutex );
+ }
+ else
+#endif
+ {
+ k++;
+ if( k == nq ) break;
+ fprintf( stderr, "\r%d / %d \r", k, nq );
+ }
+
+ if( alg == 'R' ) // if 'r' -> calllast_fast
+ {
+ klim = MIN( k, njob-nadd );
+// klim = k; // dochira demo yoi
+ if( klim == k )
+ {
+ sprintf( command, "_db%dd", k );
+ lfp = fopen( command, "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open _db." );
+ exit( 1 );
+ }
+ for( i=0; i<klim; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
+ fclose( lfp );
+
+// sprintf( command, "md5sum _db%dd > /dev/tty", k );
+// system( command );
+
+ if( dorp == 'd' )
+ sprintf( command, "%s/lastdb _db%dd _db%dd", whereispairalign, k, k );
+ else
+ sprintf( command, "%s/lastdb -p _db%dd _db%dd", whereispairalign, k, k );
+ system( command );
+ sprintf( kd, "%d", k );
+ }
+ else // calllast_fast de tsukutta nowo riyou
+ {
+ kd[0] = 0;
+// fprintf( stderr, "klim=%d, njob=%d, nadd=%d, skip!\n", klim, njob, nadd );
+ }
+ }
+ else // 'r'
+ {
+ kd[0] = 0;
+ }
+
+ sprintf( command, "_q%d", k );
+ lfp = fopen( command, "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open %s", command );
+ exit( 1 );
+ }
+ fprintf( lfp, ">%d\n%s\n", k, qseq[k] );
+ fclose( lfp );
+
+// if( alg == 'R' ) msize = MAX(10,k+nq);
+// else msize = MAX(10,nd+nq);
+ if( alg == 'R' ) msize = MAX(10,k*lastm);
+ else msize = MAX(10,nd*lastm);
+
+// fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
+// sprintf( command, "grep '>' _db%sd", kd );
+// system( command );
+ sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
+ if( system( command ) ) exit( 1 );
+
+ sprintf( command, "_lastres%d", k );
+ lfp = fopen( command, "r" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot read _lastres%d", k );
+ exit( 1 );
+ }
+// readlastres( lfp, nd, nq, lastres, dseq, qseq );
+// fprintf( stderr, "Reading lastres\n" );
+ readlastresx_singleq( lfp, nd, k, lastresx );
+ fclose( lfp );
+ }
+ return( NULL );
+}
+
+
+static void calllast_fast( int nd, char **dseq, int nq, char **qseq, Lastresx **lastresx )
{
- static FILE *fp = NULL;
- float value;
- char *aln1;
- char *aln2;
- int of1tmp, of2tmp;
+ int i, j;
+ FILE *lfp;
+ char command[1000];
- if( fp == NULL )
+ lfp = fopen( "_scoringmatrixforlast", "w" );
+ if( !lfp )
{
- fp = fopen( "_foldalignout", "r" );
- if( fp == NULL )
+ fprintf( stderr, "Cannot open _scoringmatrixforlast" );
+ exit( 1 );
+ }
+ if( dorp == 'd' )
+ {
+ fprintf( lfp, " " );
+ for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
+ fprintf( lfp, "\n" );
+ for( i=0; i<4; i++ )
{
- fprintf( stderr, "Cannot open _foldalignout\n" );
- exit( 1 );
+ fprintf( lfp, "%c ", amino[i] );
+ for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
+ fprintf( lfp, "\n" );
+ }
+ }
+ else
+ {
+ fprintf( lfp, " " );
+ for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
+ fprintf( lfp, "\n" );
+ for( i=0; i<20; i++ )
+ {
+ fprintf( lfp, "%c ", amino[i] );
+ for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
+ fprintf( lfp, "\n" );
}
}
+ fclose( lfp );
- aln1 = calloc( alloclen, sizeof( char ) );
- aln2 = calloc( alloclen, sizeof( char ) );
+// if( alg == 'r' ) // if 'R' -> lastcallthread, kokonoha nadd>0 no toki nomi shiyou
+ {
+ sprintf( command, "_dbd" );
+ lfp = fopen( command, "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open _dbd" );
+ exit( 1 );
+ }
+ if( alg == 'R' )
+ j = njob-nadd;
+ else
+ j = nd;
+ for( i=0; i<j; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
- readpairfoldalign( fp, *mseq1, *mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
+ fclose( lfp );
+ if( dorp == 'd' )
+ sprintf( command, "%s/lastdb _dbd _dbd", whereispairalign );
+ else
+ sprintf( command, "%s/lastdb -p _dbd _dbd", whereispairalign );
+ system( command );
+ }
- if( strstr( foldalignopt, "-global") )
+#ifdef enablemultithread
+ if( nthread )
{
- fprintf( stderr, "Calling G__align11\n" );
- value = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
- *of1pt = 0;
- *of2pt = 0;
+ pthread_t *handle;
+ pthread_mutex_t mutex;
+ lastcallthread_arg_t *targ;
+ int *ksharept;
+ targ = (lastcallthread_arg_t *)calloc( nthread, sizeof( lastcallthread_arg_t ) );
+ handle = calloc( nthread, sizeof( pthread_t ) );
+ ksharept = calloc( 1, sizeof(int) );
+ *ksharept = 0;
+ pthread_mutex_init( &mutex, NULL );
+ for( i=0; i<nthread; i++ )
+ {
+ targ[i].thread_no = i;
+ targ[i].kshare = ksharept;
+ targ[i].nq = nq;
+ targ[i].nd = nd;
+ targ[i].dseq = dseq;
+ targ[i].qseq = qseq;
+ targ[i].lastresx = lastresx;
+ targ[i].mutex = &mutex;
+ pthread_create( handle+i, NULL, lastcallthread, (void *)(targ+i) );
+ }
+
+ for( i=0; i<nthread; i++ )
+ {
+ pthread_join( handle[i], NULL );
+ }
+ pthread_mutex_destroy( &mutex );
+ free( handle );
+ free( targ );
+ free( ksharept );
}
else
+#endif
{
- fprintf( stderr, "Calling L__align11\n" );
- value = L__align11( mseq1, mseq2, alloclen, of1pt, of2pt );
+ lastcallthread_arg_t *targ;
+ targ = (lastcallthread_arg_t *)calloc( 1, sizeof( lastcallthread_arg_t ) );
+ targ[0].nq = nq;
+ targ[0].nd = nd;
+ targ[0].dseq = dseq;
+ targ[0].qseq = qseq;
+ targ[0].lastresx = lastresx;
+ lastcallthread( targ );
+ free( targ );
}
-// value = (float)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
+}
- if( aln1[0] == 0 )
+static void calllast_once( int nd, char **dseq, int nq, char **qseq, Lastresx **lastresx )
+{
+ int i, j;
+ char command[5000];
+ FILE *lfp;
+ int msize;
+ int res;
+
+ fprintf( stderr, "nq=%d\n", nq );
+
+ lfp = fopen( "_db", "w" );
+ if( !lfp )
{
- fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d. Sequence alignment is used instead.\n", m1+1, m2+1 );
+ fprintf( stderr, "Cannot open _db" );
+ exit( 1 );
+ }
+ for( i=0; i<nd; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
+ fclose( lfp );
+
+ if( dorp == 'd' )
+ {
+ sprintf( command, "%s/lastdb _db _db", whereispairalign );
+ system( command );
+ lfp = fopen( "_scoringmatrixforlast", "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open _scoringmatrixforlast" );
+ exit( 1 );
+ }
+ fprintf( lfp, " " );
+ for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
+ fprintf( lfp, "\n" );
+ for( i=0; i<4; i++ )
+ {
+ fprintf( lfp, "%c ", amino[i] );
+ for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
+ fprintf( lfp, "\n" );
+ }
+ fclose( lfp );
+#if 0
+ sprintf( command, "lastex -s 2 -a %d -b %d -p _scoringmatrixforlast -E 10000 _db.prj _db.prj > _lastex", -penalty, -penalty_ex );
+ system( command );
+ lfp = fopen( "_lastex", "r" );
+ fgets( command, 4999, lfp );
+ fgets( command, 4999, lfp );
+ fgets( command, 4999, lfp );
+ fgets( command, 4999, lfp );
+ laste = atoi( command );
+ fclose( lfp );
+ fprintf( stderr, "laste = %d\n", laste );
+ sleep( 10 );
+#else
+// laste = 5000;
+#endif
}
else
{
- strcpy( *mseq1, aln1 );
- strcpy( *mseq2, aln2 );
- *of1pt = of1tmp;
- *of2pt = of2tmp;
+ sprintf( command, "%s/lastdb -p _db _db", whereispairalign );
+ system( command );
+ lfp = fopen( "_scoringmatrixforlast", "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open _scoringmatrixforlast" );
+ exit( 1 );
+ }
+ fprintf( lfp, " " );
+ for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
+ fprintf( lfp, "\n" );
+ for( i=0; i<20; i++ )
+ {
+ fprintf( lfp, "%c ", amino[i] );
+ for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
+ fprintf( lfp, "\n" );
+ }
+ fclose( lfp );
+// fprintf( stderr, "Not written yet\n" );
}
-// value = naivepairscore11( *mseq1, *mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
-
-// fclose( fp ); // saigo dake yatta houga yoi.
-
-// fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
-// fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
+ lfp = fopen( "_q", "w" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot open _q" );
+ exit( 1 );
+ }
+ for( i=0; i<nq; i++ )
+ {
+ fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
+ }
+ fclose( lfp );
+ msize = MAX(10,nd*lastm);
- free( aln1 );
- free( aln2 );
+// fprintf( stderr, "Calling lastal from calllast_once, msize=%d\n", msize );
+ sprintf( command, "%s/lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", whereispairalign, msize, laste, -penalty, -penalty_ex );
+// sprintf( command, "lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", 1, laste, -penalty, -penalty_ex );
+// sprintf( command, "lastal -v -e 40 -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", -penalty, -penalty_ex );
+ res = system( command );
+ if( res )
+ {
+ fprintf( stderr, "LAST aborted\n" );
+ exit( 1 );
+ }
- return( value );
+ lfp = fopen( "_lastres", "r" );
+ if( !lfp )
+ {
+ fprintf( stderr, "Cannot read _lastres" );
+ exit( 1 );
+ }
+// readlastres( lfp, nd, nq, lastres, dseq, qseq );
+ fprintf( stderr, "Reading lastres\n" );
+ readlastresx( lfp, nd, nq, lastresx, dseq, qseq );
+ fclose( lfp );
}
static void callfoldalign( int nseq, char **mseq )
}
}
-static float recalllara( char **mseq1, char **mseq2, int alloclen )
+static double recalllara( char **mseq1, char **mseq2, int alloclen )
{
static FILE *fp = NULL;
static char *ungap1;
static char *ori2;
// int res;
static char com[10000];
- float value;
+ double value;
if( fp == NULL )
exit( 1 );
}
- value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
+ value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
// fclose( fp ); // saigo dake yatta houga yoi.
}
-static float callmxscarna_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
+static double calldafs_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
{
FILE *fp;
int res;
char *com;
- float value;
+ double value;
char *dirname;
t2u( *mseq1 );
t2u( *mseq2 );
- sprintf( com, "%s/_mxscarnainorg", dirname );
+ sprintf( com, "%s/_dafsinorg", dirname );
fp = fopen( com, "w" );
if( !fp )
{
- fprintf( stderr, "Cannot open %s/_mxscarnainorg\n", dirname );
+ fprintf( stderr, "Cannot open %s/_dafsinorg\n", dirname );
exit( 1 );
}
fprintf( fp, ">1\n" );
write1seq( fp, *mseq2 );
fclose( fp );
- sprintf( com, "tr -d '\\r' < %s/_mxscarnainorg > %s/_mxscarnain", dirname, dirname );
+ sprintf( com, "tr -d '\\r' < %s/_dafsinorg > %s/_dafsin", dirname, dirname );
system( com ); // for cygwin, wakaran
-#if 0
- sprintf( com, "cd %s; %s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum", dirname, whereispairalign );
-#else
- sprintf( com, "_mxscarnash%s", dirname );
+ sprintf( com, "_dafssh%s", dirname );
fp = fopen( com, "w" );
fprintf( fp, "cd %s\n", dirname );
- fprintf( fp, "%s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum\n", whereispairalign );
+ fprintf( fp, "%s/dafs --mafft-in _bpp _dafsin > _dafsout 2>_dum\n", whereispairalign );
fprintf( fp, "exit $tatus\n" );
fclose( fp );
- sprintf( com, "tr -d '\\r' < _mxscarnash%s > _mxscarnash%s.unix", dirname, dirname );
+ sprintf( com, "tr -d '\\r' < _dafssh%s > _dafssh%s.unix", dirname, dirname );
system( com ); // for cygwin, wakaran
- sprintf( com, "sh _mxscarnash%s.unix 2>_dum%s", dirname, dirname );
-#endif
+ sprintf( com, "sh _dafssh%s.unix 2>_dum%s", dirname, dirname );
res = system( com );
if( res )
{
- fprintf( stderr, "Error in mxscarna\n" );
+ fprintf( stderr, "Error in dafs\n" );
exit( 1 );
}
- sprintf( com, "%s/_mxscarnaout", dirname );
+ sprintf( com, "%s/_dafsout", dirname );
fp = fopen( com, "r" );
if( !fp )
{
- fprintf( stderr, "Cannot open %s/_mxscarnaout\n", dirname );
+ fprintf( stderr, "Cannot open %s/_dafsout\n", dirname );
exit( 1 );
}
+ myfgets( com, 999, fp ); // nagai kanousei ga arunode
+ fgets( com, 999, fp );
+ myfgets( com, 999, fp ); // nagai kanousei ga arunode
fgets( com, 999, fp );
load1SeqWithoutName_new( fp, *mseq1 );
fgets( com, 999, fp );
// fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
// fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
- value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
+ value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
#if 0
sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
return( value );
}
-#if 0
-static float callmxscarna_slow( char **mseq1, char **mseq2, int alloclen )
+static double callmxscarna_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
{
FILE *fp;
int res;
- static char com[10000];
- float value;
+ char *com;
+ double value;
+ char *dirname;
+
+
+ dirname = calloc( 100, sizeof( char ) );
+ com = calloc( 1000, sizeof( char ) );
+ sprintf( dirname, "_%d-%d", i, j );
+ sprintf( com, "rm -rf %s", dirname );
+ system( com );
+ sprintf( com, "mkdir %s", dirname );
+ system( com );
+
+
+ sprintf( com, "%s/_bpporg", dirname );
+ fp = fopen( com, "w" );
+ if( !fp )
+ {
+ fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
+ exit( 1 );
+ }
+ fprintf( fp, ">a\n" );
+ while( *bpp1 )
+ fprintf( fp, "%s", *bpp1++ );
+
+ fprintf( fp, ">b\n" );
+ while( *bpp2 )
+ fprintf( fp, "%s", *bpp2++ );
+ fclose( fp );
+ sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
+ system( com ); // for cygwin, wakaran
t2u( *mseq1 );
t2u( *mseq2 );
- fp = fopen( "_mxscarnain", "w" );
+
+ sprintf( com, "%s/_mxscarnainorg", dirname );
+ fp = fopen( com, "w" );
if( !fp )
{
- fprintf( stderr, "Cannot open _mxscarnain\n" );
+ fprintf( stderr, "Cannot open %s/_mxscarnainorg\n", dirname );
exit( 1 );
}
fprintf( fp, ">1\n" );
- fprintf( fp, "%s\n", *mseq1 );
+// fprintf( fp, "%s\n", *mseq1 );
+ write1seq( fp, *mseq1 );
fprintf( fp, ">2\n" );
- fprintf( fp, "%s\n", *mseq2 );
+// fprintf( fp, "%s\n", *mseq2 );
+ write1seq( fp, *mseq2 );
+ fclose( fp );
+
+ sprintf( com, "tr -d '\\r' < %s/_mxscarnainorg > %s/_mxscarnain", dirname, dirname );
+ system( com ); // for cygwin, wakaran
+
+#if 0
+ sprintf( com, "cd %s; %s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum", dirname, whereispairalign );
+#else
+ sprintf( com, "_mxscarnash%s", dirname );
+ fp = fopen( com, "w" );
+ fprintf( fp, "cd %s\n", dirname );
+ fprintf( fp, "%s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum\n", whereispairalign );
+ fprintf( fp, "exit $tatus\n" );
fclose( fp );
+//sleep( 10000 );
+
+ sprintf( com, "tr -d '\\r' < _mxscarnash%s > _mxscarnash%s.unix", dirname, dirname );
+ system( com ); // for cygwin, wakaran
- sprintf( com, "env PATH=%s mxscarnamod _mxscarnain > _mxscarnaout 2>_dum", whereispairalign );
+ sprintf( com, "sh _mxscarnash%s.unix 2>_dum%s", dirname, dirname );
+#endif
res = system( com );
if( res )
{
exit( 1 );
}
- fp = fopen( "_mxscarnaout", "r" );
+ sprintf( com, "%s/_mxscarnaout", dirname );
+
+ fp = fopen( com, "r" );
if( !fp )
{
- fprintf( stderr, "Cannot open _mxscarnaout\n" );
+ fprintf( stderr, "Cannot open %s/_mxscarnaout\n", dirname );
exit( 1 );
}
// fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
// fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
- value = (float)naivepairscore11( *mseq1, *mseq2, penalty );
+ value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
+
+#if 0
+ sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
+ if( system( com ) )
+ {
+ fprintf( stderr, "retrying to rmdir\n" );
+ usleep( 2000 );
+ system( com );
+ }
+#endif
+
+ free( dirname );
+ free( com );
+
return( value );
}
-#endif
static void readhat4( FILE *fp, char ***bpp )
{
}
fgets( oneline, 999, fp );
// fprintf( stderr, "oneline=%s\n", oneline );
-// sscanf( oneline, "%d %d %f", &posi, &posj, &prob );
+// sscanf( oneline, "%d %d %lf", &posi, &posj, &prob );
// fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
*bpp = realloc( *bpp, (bppsize+2) * sizeof( char * ) );
(*bpp)[bppsize] = calloc( 100, sizeof( char ) );
fclose( fp );
}
-void arguments( int argc, char *argv[] )
+static void arguments( int argc, char *argv[] )
{
int c;
nthread = 1;
+ laste = 5000;
+ lastm = 3;
+ nadd = 0;
+ lastsubopt = 0;
+ lastonce = 0;
foldalignopt[0] = 0;
laraparams = NULL;
inputfile = NULL;
stdout_dist = 0;
store_dist = 1;
store_localhom = 1;
- dorp = NOTSPECIFIED;
+// dorp = NOTSPECIFIED;
ppenalty = NOTSPECIFIED;
ppenalty_OP = NOTSPECIFIED;
ppenalty_ex = NOTSPECIFIED;
ppenalty_EX = NOTSPECIFIED;
+ penalty_shift_factor = 1000.0;
poffset = NOTSPECIFIED;
kimuraR = NOTSPECIFIED;
pamN = NOTSPECIFIED;
fftThreshold = NOTSPECIFIED;
RNAppenalty = NOTSPECIFIED;
RNApthr = NOTSPECIFIED;
-
+ specificityconsideration = 0.0;
+ usenaivescoreinsteadofalignmentscore = 0;
+ specifictarget = 0;
+ nwildcard = 0;
+
+// reporterr( "argc=%d\n", argc );
+// reporterr( "*argv=%s\n", *argv );
+// reporterr( "(*argv)[0]=%c\n", (*argv)[0] );
while( --argc > 0 && (*++argv)[0] == '-' )
{
+// reporterr( "(*argv)[0] in while loop = %s\n", (*argv) );
while ( ( c = *++argv[0] ) )
{
switch( c )
{
case 'i':
inputfile = *++argv;
- fprintf( stderr, "inputfile = %s\n", inputfile );
+// fprintf( stderr, "inputfile = %s\n", inputfile );
--argc;
goto nextoption;
case 'f':
ppenalty_EX = (int)( atof( *++argv ) * 1000 - 0.5 );
--argc;
goto nextoption;
+ case 'Q':
+ penalty_shift_factor = atof( *++argv );
+ --argc;
+ goto nextoption;
case 'h':
poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
--argc;
goto nextoption;
case 'k':
- kimuraR = atoi( *++argv );
+ kimuraR = myatoi( *++argv );
// fprintf( stderr, "kimuraR = %d\n", kimuraR );
--argc;
goto nextoption;
case 'b':
- nblosum = atoi( *++argv );
+ nblosum = myatoi( *++argv );
scoremtx = 1;
// fprintf( stderr, "blosum %d\n", nblosum );
--argc;
goto nextoption;
case 'j':
- pamN = atoi( *++argv );
+ pamN = myatoi( *++argv );
scoremtx = 0;
TMorJTT = JTT;
- fprintf( stderr, "jtt %d\n", pamN );
+// fprintf( stderr, "jtt %d\n", pamN );
--argc;
goto nextoption;
case 'm':
- pamN = atoi( *++argv );
+ pamN = myatoi( *++argv );
scoremtx = 0;
TMorJTT = TM;
- fprintf( stderr, "TM %d\n", pamN );
+// fprintf( stderr, "TM %d\n", pamN );
--argc;
goto nextoption;
+#if 0
case 'l':
ppslocal = (int)( atof( *++argv ) * 1000 + 0.5 );
pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
// fprintf( stderr, "pslocal = %d\n", pslocal );
--argc;
goto nextoption;
+#else
+ case 'l':
+ if( atof( *++argv ) < 0.00001 ) store_localhom = 0;
+ --argc;
+ goto nextoption;
+#endif
case 'd':
whereispairalign = *++argv;
fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
--argc;
goto nextoption;
case 'C':
- nthread = atoi( *++argv );
- fprintf( stderr, "nthread = %d\n", nthread );
+ nthread = myatoi( *++argv );
+// fprintf( stderr, "nthread = %d\n", nthread );
+ --argc;
+#ifndef enablemultithread
+ nthread = 0;
+#endif
+ goto nextoption;
+ case 'I':
+ nadd = myatoi( *++argv );
+// fprintf( stderr, "nadd = %d\n", nadd );
+ --argc;
+ goto nextoption;
+ case 'w':
+ lastm = myatoi( *++argv );
+ fprintf( stderr, "lastm = %d\n", lastm );
+ --argc;
+ goto nextoption;
+ case 'e':
+ laste = myatoi( *++argv );
+ fprintf( stderr, "laste = %d\n", laste );
+ --argc;
+ goto nextoption;
+ case 'u':
+ specificityconsideration = (double)myatof( *++argv );
+// fprintf( stderr, "specificityconsideration = %f\n", specificityconsideration );
--argc;
goto nextoption;
+ case 'K': // Hontou ha iranai. disttbfast.c, tbfast.c to awaserutame.
+ break;
case 'c':
stdout_dist = 1;
break;
fmodel = 1;
break;
#endif
+#if 0
case 'r':
fmodel = -1;
break;
+#endif
case 'D':
dorp = 'd';
break;
case 'P':
dorp = 'p';
break;
+#if 0
case 'e':
fftscore = 0;
break;
-#if 0
case 'O':
fftNoAnchStop = 1;
break;
#endif
+#if 0
case 'Q':
calledByXced = 1;
break;
-#if 0
case 'x':
disp = 1;
break;
case 'a':
alg = 'a';
break;
-#endif
case 'S':
alg = 'S';
break;
+#endif
+ case 'U':
+ lastonce = 1;
+ break;
+ case 'S':
+ lastsubopt = 1;
+ break;
case 't':
alg = 't';
store_localhom = 0;
case 'L':
alg = 'L';
break;
+ case 'Y':
+ alg = 'Y'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> L;
+ break;
+ case 'Z':
+ usenaivescoreinsteadofalignmentscore = 1;
+ break;
case 's':
alg = 's';
break;
+ case 'G':
+ alg = 'G';
+ break;
case 'B':
alg = 'B';
break;
case 'R':
alg = 'R';
break;
+ case 'r':
+ alg = 'r'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> R, last
+ break;
case 'N':
alg = 'N';
break;
- case 'K':
- alg = 'K';
- break;
case 'A':
alg = 'A';
break;
case 'y':
divpairscore = 1;
break;
+ case '=':
+ specifictarget = 1;
+ break;
+ case ':':
+ nwildcard = 1;
+ break;
/* Modified 01/08/27, default: user tree */
case 'J':
tbutree = 0;
fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
--argc;
goto nextoption;
+#if 0
case 'z':
- fftThreshold = atoi( *++argv );
+ fftThreshold = myatoi( *++argv );
--argc;
goto nextoption;
case 'w':
- fftWinSize = atoi( *++argv );
+ fftWinSize = myatoi( *++argv );
--argc;
goto nextoption;
case 'Z':
checkC = 1;
break;
+#endif
default:
fprintf( stderr, "illegal option %c\n", c );
argc = 0;
}
if( argc != 0 )
{
- fprintf( stderr, "options: Check source file !\n" );
+ fprintf( stderr, "pairlocalalign options: Check source file !\n" );
exit( 1 );
}
if( tbitr == 1 && outgap == 0 )
return( val );
}
+static double score2dist( double pscore, double selfscore1, double selfscore2)
+{
+ double val;
+ double bunbo;
+// fprintf( stderr, "In score2dist\n" );
+
+ if( (bunbo=MIN( selfscore1, selfscore2 )) == 0.0 )
+ val = 2.0;
+ else if( bunbo < pscore ) // mondai ari
+ val = 0.0;
+ else
+ val = ( 1.0 - pscore / bunbo ) * 2.0;
+ return( val );
+}
+
#if enablemultithread
-static void *athread( void *arg )
+static void *athread( void *arg ) // alg='R', alg='r' -> tsukawarenai.
{
thread_arg_t *targ = (thread_arg_t *)arg;
- int i, j;
- int clus1, clus2;
- int off1, off2;
+ int i, ilim, j, jst;
+ int off1, off2, dum1, dum2, thereisx;
int intdum;
- double bunbo;
- float pscore = 0.0; // by D.Mathog
- double *effarr;
+ double pscore = 0.0; // by D.Mathog
double *effarr1;
double *effarr2;
- char **pair;
- char *indication1, *indication2;
- char **mseq1, **mseq2;
+ char **mseq1, **mseq2, **distseq1, **distseq2, **dumseq1, **dumseq2;
char **aseq;
+ double **dynamicmtx = NULL;
+ double dist;
+ double scoreoffset;
// thread_arg
int thread_no = targ->thread_no;
Jobtable *jobpospt = targ->jobpospt;
char **name = targ->name;
char **seq = targ->seq;
+ char **dseq = targ->dseq;
+ int *thereisxineachseq = targ->thereisxineachseq;
LocalHom **localhomtable = targ->localhomtable;
double **distancemtx = targ->distancemtx;
double *selfscore = targ->selfscore;
char ***bpp = targ->bpp;
+ Lastresx **lastresx = targ->lastresx;
int alloclen = targ->alloclen;
+ int *targetmap = targ->targetmap;
// fprintf( stderr, "thread %d start!\n", thread_no );
- effarr = AllocateDoubleVec( njob );
- for( i=0; i<njob; i++ ) effarr[i] = 1.0;
- effarr1 = AllocateDoubleVec( njob );
- effarr2 = AllocateDoubleVec( njob );
- indication1 = AllocateCharVec( 150 );
- indication2 = AllocateCharVec( 150 );
- pair = AllocateCharMtx( njob, njob );
- for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
- for( i=0; i<njob; i++ ) pair[i][i] = 1;
+ effarr1 = AllocateDoubleVec( 1 );
+ effarr2 = AllocateDoubleVec( 1 );
mseq1 = AllocateCharMtx( njob, 0 );
mseq2 = AllocateCharMtx( njob, 0 );
- aseq = AllocateCharMtx( njob, alloclen+10 );
+ if( alg == 'N' )
+ {
+ dumseq1 = AllocateCharMtx( 1, alloclen+10 );
+ dumseq2 = AllocateCharMtx( 1, alloclen+10 );
+ }
+ distseq1 = AllocateCharMtx( 1, 0 );
+ distseq2 = AllocateCharMtx( 1, 0 );
+ aseq = AllocateCharMtx( 2, alloclen+10 );
+ if( specificityconsideration > 0.0 ) dynamicmtx = AllocateDoubleMtx( nalphabets, nalphabets );
+
+ if( alg == 'Y' || alg == 'r' ) ilim = njob - nadd;
+ else ilim = njob - 1;
+
while( 1 )
{
if( j == njob )
{
i++;
- j = i + 1;
- if( i == njob-1 )
+
+ if( alg == 'Y' || alg == 'r' ) jst = njob - nadd;
+ else jst = i + 1;
+ j = jst;
+
+ if( i == ilim )
{
// fprintf( stderr, "thread %d end!\n", thread_no );
pthread_mutex_unlock( targ->mutex_counter );
commonIP = NULL;
if( commonJP ) FreeIntMtx( commonJP );
commonJP = NULL;
- Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
+ Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
+ G__align11( NULL, NULL, NULL, 0, 0, 0 ); // 20130603
G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
- L__align11( NULL, NULL, 0, NULL, NULL );
- genL__align11( NULL, NULL, 0, NULL, NULL );
- free( effarr );
+ L__align11( NULL, 0.0, NULL, NULL, 0, NULL, NULL );
+ L__align11_noalign( NULL, NULL, NULL );
+ genL__align11( NULL, NULL, NULL, 0, NULL, NULL );
free( effarr1 );
free( effarr2 );
- free( indication1 );
- free( indication2 );
- FreeCharMtx( pair );
free( mseq1 );
free( mseq2 );
+ if( alg == 'N' )
+ {
+ FreeCharMtx( dumseq1 );
+ FreeCharMtx( dumseq2 );
+ }
+ free( distseq1 );
+ free( distseq2 );
FreeCharMtx( aseq );
+ if( dynamicmtx ) FreeDoubleMtx( dynamicmtx );
return( NULL );
}
}
pthread_mutex_unlock( targ->mutex_counter );
- if( j == i+1 || j % 100 == 0 )
+// if( j == i+1 || j % 100 == 0 )
+ if( j == i+1 && i % 10 == 0 )
{
- fprintf( stderr, "% 5d / %d (by thread %3d) \r", i, njob, thread_no );
+ fprintf( stderr, "% 5d / %d (by thread %3d) \r", i, njob-nadd, thread_no );
// fprintf( stderr, "% 5d - %5d / %d (thread %d)\n", i, j, njob, thread_no );
}
if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
{
- if( store_dist ) distancemtx[i][j] = 2.0;
+ if( store_dist )
+ {
+ if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
+ else distancemtx[i][j-i] = 3.0;
+ }
if( stdout_dist)
{
pthread_mutex_lock( targ->mutex_stdout );
- fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 2.0 );
+ fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
pthread_mutex_unlock( targ->mutex_stdout );
}
continue;
}
- strcpy( aseq[i], seq[i] );
- strcpy( aseq[j], seq[j] );
- clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
- clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
+ strcpy( aseq[0], seq[i] );
+ strcpy( aseq[1], seq[j] );
+// clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
+// clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
+// fprintf( stderr, "Skipping conjuction..\n" );
+
+ effarr1[0] = 1.0;
+ effarr2[0] = 1.0;
+ mseq1[0] = aseq[0];
+ mseq2[0] = aseq[1];
+
+ thereisx = thereisxineachseq[i] + thereisxineachseq[j];
+// strcpy( distseq1[0], dseq[i] ); // nen no tame
+// strcpy( distseq2[0], dseq[j] ); // nen no tame
+ distseq1[0] = dseq[i];
+ distseq2[0] = dseq[j];
+
// fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
// fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
if( use_fft )
{
- pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
+ pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, 1, 1, alloclen, &intdum, NULL, 0, NULL );
// fprintf( stderr, "pscore (fft) = %f\n", pscore );
off1 = off2 = 0;
}
{
switch( alg )
{
+ case( 'R' ):
+ if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
+ pscore = 0.0;
+ else
+ pscore = (double)lastresx[i][j].score; // all pair
+ break;
+ case( 'r' ):
+ if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) )
+ pscore = (double)lastresx[i][j-(njob-nadd)].score;
+ else
+ pscore = 0.0;
+ break;
case( 'L' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
- L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
+ if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
+ pscore = 0.0;
+ else
+ {
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+// if( store_localhom )
+ if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
+ {
+ pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+ if( ( scoreoffset = dist2offset( dist ) ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ L__align11( dynamicmtx, scoreoffset, mseq1, mseq2, alloclen, &off1, &off2 );
+ }
+ }
+#endif
+ }
+ else
+ pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
+ }
+ }
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
+ break;
+ case( 'Y' ):
+ if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
+ {
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+ if( store_localhom )
+ {
+ pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
+ }
+ else
+ pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
+ }
+ }
+ else
+ pscore = 0.0;
break;
case( 'A' ):
- pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+// if( store_localhom )
+ if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
+ {
+ pscore = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
+ if( thereisx ) pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+// dist = score2dist( L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ), selfscore[i], selfscore[j] ); // 2014/Feb/20
+ if( dist2offset( dist ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ G__align11( dynamicmtx, mseq1, mseq2, alloclen, outgap, outgap );
+
+ }
+// pscore = (double)naivepairscore11( *mseq1, *mseq2, 0.0 );
+ }
+#endif
+ }
+ else
+ pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
+ }
off1 = off2 = 0;
break;
case( 'N' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
- genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
+ pscore = genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
+ if( thereisx )
+ {
+ strcpy( dumseq1[0], distseq1[0] );
+ strcpy( dumseq2[0], distseq2[0] );
+ pscore = genL__align11( n_dis_consweight_multi, dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
+ }
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+ if( dist2offset( dist ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ genL__align11( dynamicmtx, mseq1, mseq2, alloclen, &off1, &off2 );
+ }
+ }
+#endif
+ }
break;
case( 't' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
off1 = off2 = 0;
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
+ pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
break;
case( 's' ):
pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
off1 = off2 = 0;
break;
+ case( 'G' ):
+ pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
+ off1 = off2 = 0;
+ break;
#if 0
case( 'a' ):
pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
off1 = off2 = 0;
break;
case( 'K' ):
- pscore = genG__align11( mseq1, mseq2, alloclen );
+ pscore = genG__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen );
off1 = off2 = 0;
break;
case( 'H' ):
#if SCOREOUT
fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
#endif
- if( !store_localhom )
- ;
- else if( alg == 'H' )
- putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
- else if( alg != 'S' && alg != 'V' )
+// if( pscore > 0.0 && ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) ) x-ins-i de seido teika
+ if( ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) )
{
- putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
+ if( !store_localhom )
+ ;
+ else if( specifictarget && targetmap[i] == -1 && targetmap[j] == -1)
+ ;
+ else if( alg == 'R' )
+ putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j, 'h' );
+ else if( alg == 'r' )
+ putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd), 'h' );// ?????
+ else if( alg == 'H' )
+ putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else if( alg == 'Y' )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else if( !specifictarget && alg != 'S' && alg != 'V' )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-i, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else
+// putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
+ {
+ if( targetmap[i] != -1 && targetmap[j] != -1 )
+ {
+ putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' ); // sukoshi muda.
+ }
+ else if( targetmap[j] != -1 )
+ putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
+ else if( targetmap[i] != -1 )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+#if 0
+ if( targetmap[i] != -1 )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+
+ else if( targetmap[j] != -1 )
+ putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
+#endif
+ else
+ {
+ reporterr( "okashii\n" );
+ exit( 1 );
+ }
+ }
}
+ pscore = score2dist( pscore, selfscore[i], selfscore[j] );
+
+// pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 );
+// pscore = score2dist( pscore, selfscore[i], selfscore[j] );
+// reporterr( "->pscore = %f\n", pscore );
- if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 || bunbo < pscore )
- pscore = 2.0;
- else
- pscore = ( 1.0 - pscore / bunbo ) * 2.0;
}
else
{
pthread_mutex_unlock( targ->mutex_stdout );
}
#endif // mutex
- if( store_dist) distancemtx[i][j] = pscore;
+ if( store_dist )
+ {
+ if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
+ else distancemtx[i][j-i] = pscore;
+ }
}
}
#endif
-static void pairalign( char **name, int nlen[M], char **seq, char **aseq, char **mseq1, char **mseq2, double *effarr, int alloclen )
+static void pairalign( char **name, int *nlen, char **seq, char **aseq, char **dseq, int *thereisxineachseq, char **mseq1, char **mseq2, int alloclen, Lastresx **lastresx, double **distancemtx, LocalHom **localhomtable, int ngui )
{
- int i, j, ilim;
- int clus1, clus2;
- int off1, off2;
- float pscore = 0.0; // by D.Mathog
- static char *indication1, *indication2;
+ int i, j, ilim, jst, jj;
+ int off1, off2, dum1, dum2, thereisx;
+ double pscore = 0.0; // by D.Mathog
FILE *hat2p, *hat3p;
- double **distancemtx;
+// double **distancemtx;
double *selfscore;
double *effarr1;
double *effarr2;
char *pt;
char *hat2file = "hat2";
- LocalHom **localhomtable = NULL, *tmpptr;
- static char **pair;
+// LocalHom **localhomtable = NULL,
+ LocalHom *tmpptr;
int intdum;
- double bunbo;
char ***bpp = NULL; // mxscarna no toki dake
+ char **distseq1, **distseq2;
+ char **dumseq1, **dumseq2;
+ double dist;
+ double scoreoffset;
+ int ntarget;
+ int *targetmap, *targetmapr;
+
+
+ if( specifictarget )
+ {
+ targetmap = calloc( njob, sizeof( int ) );
+ ntarget = 0;
+ for( i=0; i<njob; i++ )
+ {
+ targetmap[i] = -1;
+ if( !strncmp( name[i]+1, "_focus_", 7 ) )
+ targetmap[i] = ntarget++;
+ }
+ targetmapr = calloc( ntarget, sizeof( int ) );
+ for( i=0; i<njob; i++ )
+ if( targetmap[i] != -1 ) targetmapr[targetmap[i]] = i;
+
+ if( ntarget == 0 )
+ {
+ reporterr( "\n\nAdd '>_focus_' to the title lines of the sequences to be focused on.\n\n" );
+ exit( 1 );
+ }
+ else
+ {
+ reporterr( "nfocus = %d \n", ntarget );
+ }
+ }
+ else
+ {
+ ntarget = njob;
+ targetmap = calloc( njob, sizeof( int ) );
+ targetmapr = calloc( njob, sizeof( int ) );
+ for( i=0; i<njob; i++ )
+ targetmap[i] = targetmapr[i] = i;
+ }
+
+#if 0
+ for( i=0; i<njob; i++ )
+ reporterr( "targetmap[%d] = %d\n", i, targetmap[i] );
+ for( i=0; i<ntarget; i++ )
+ reporterr( "targetmapr[%d] = %d\n", i, targetmapr[i] );
+#endif
- if( store_localhom )
+ if( store_localhom && localhomtable == NULL )
{
- localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
- for( i=0; i<njob; i++)
+ if( alg == 'Y' || alg == 'r' )
+ {
+ ilim = njob - nadd;
+ jst = nadd;
+ }
+ else
{
- localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
- for( j=0; j<njob; j++)
+ ilim = ntarget;
+ jst = njob;
+ }
+ localhomtable = (LocalHom **)calloc( ilim, sizeof( LocalHom *) );
+ for( i=0; i<ilim; i++)
+ {
+ localhomtable[i] = (LocalHom *)calloc( jst, sizeof( LocalHom ) );
+ for( j=0; j<jst; j++)
{
localhomtable[i][j].start1 = -1;
localhomtable[i][j].end1 = -1;
localhomtable[i][j].opt = -1.0;
localhomtable[i][j].next = NULL;
localhomtable[i][j].nokori = 0;
+ localhomtable[i][j].extended = -1;
+ localhomtable[i][j].last = localhomtable[i]+j;
+ localhomtable[i][j].korh = 'h';
}
+ if( !specifictarget && alg != 'Y' && alg != 'r' ) jst--;
}
}
- if( store_dist ) distancemtx = AllocateDoubleMtx( njob, njob );
+ if( store_dist )
+ {
+ if( ngui == 0 )
+ {
+ if( alg == 'Y' || alg == 'r' )
+ distancemtx = AllocateDoubleMtx( njob, nadd );
+ else
+ distancemtx = AllocateDoubleHalfMtx( njob );
+// distancemtx = AllocateDoubleMtx( njob, njob );
+ }
+ }
else distancemtx = NULL;
+
+ if( alg == 'N' )
+ {
+ dumseq1 = AllocateCharMtx( 1, alloclen+10 );
+ dumseq2 = AllocateCharMtx( 1, alloclen+10 );
+ }
+ distseq1 = AllocateCharMtx( 1, 0 ); // muda
+ distseq2 = AllocateCharMtx( 1, 0 ); // muda
+
selfscore = AllocateDoubleVec( njob );
effarr1 = AllocateDoubleVec( njob );
effarr2 = AllocateDoubleVec( njob );
- indication1 = AllocateCharVec( 150 );
- indication2 = AllocateCharVec( 150 );
-#if 0
-#else
- pair = AllocateCharMtx( njob, njob );
-#endif
#if 0
fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
// writePre( njob, name, nlen, aseq, 0 );
- for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) pair[i][j] = 0;
- for( i=0; i<njob; i++ ) pair[i][i] = 1;
+ reporterr( "All-to-all alignment.\n" );
+ if( alg == 'R' )
+ {
+ fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
+ if( lastonce )
+ calllast_once( njob, seq, njob, seq, lastresx );
+ else
+ calllast_fast( njob, seq, njob, seq, lastresx );
+ fprintf( stderr, "done.\n" );
+// nthread = 0; // igo multithread nashi
+ }
+ if( alg == 'r' )
+ {
+ fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
+ fprintf( stderr, "nadd=%d\n", nadd );
+#if 1 // last_fast ha, lastdb ga muda
+ if( lastonce )
+ calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
+ else
+ calllast_fast( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
+#else
+ calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
+#endif
+
+ fprintf( stderr, "nadd=%d\n", nadd );
+ fprintf( stderr, "done.\n" );
+// nthread = 0; // igo multithread nashi
+ }
if( alg == 'H' )
{
fprintf( stderr, "done.\n" );
fprintf( stderr, "Running MXSCARNA (Tabei et al. http://www.ncrna.org/software/mxscarna)\n" );
}
+ if( alg == 'G' )
+ {
+ fprintf( stderr, "Preparing bpp\n" );
+// bpp = AllocateCharCub( njob, nlenmax, 0 );
+ bpp = calloc( njob, sizeof( char ** ) );
+ preparebpp( njob, bpp );
+ fprintf( stderr, "done.\n" );
+ fprintf( stderr, "Running DAFS (Sato et al. http://www.ncrna.org/)\n" );
+ }
for( i=0; i<njob; i++ )
{
pscore = 0.0;
for( pt=seq[i]; *pt; pt++ )
- pscore += amino_dis[(int)*pt][(int)*pt];
+ pscore += amino_dis[(unsigned char)*pt][(unsigned char)*pt];
selfscore[i] = pscore;
-
+// fprintf( stderr, "selfscore[%d] = %f\n", i, selfscore[i] );
}
#if enablemultithread
- if( nthread > 0 )
+ if( nthread > 0 ) // alg=='r' || alg=='R' -> nthread:=0 (sukoshi ue)
{
Jobtable jobpos;
pthread_t *handle;
pthread_mutex_t mutex_stdout;
thread_arg_t *targ;
+ if( alg == 'Y' || alg == 'r' ) jobpos.j = njob - nadd - 1;
+ else jobpos.j = 0;
jobpos.i = 0;
- jobpos.j = 0;
targ = calloc( nthread, sizeof( thread_arg_t ) );
handle = calloc( nthread, sizeof( pthread_t ) );
targ[i].jobpospt = &jobpos;
targ[i].name = name;
targ[i].seq = seq;
+ targ[i].dseq = dseq;
+ targ[i].thereisxineachseq = thereisxineachseq;
targ[i].localhomtable = localhomtable;
targ[i].distancemtx = distancemtx;
targ[i].selfscore = selfscore;
targ[i].bpp = bpp;
+ targ[i].lastresx = lastresx;
targ[i].alloclen = alloclen;
+ targ[i].targetmap = targetmap;
targ[i].mutex_counter = &mutex_counter;
targ[i].mutex_stdout = &mutex_stdout;
else
#endif
{
- ilim = njob - 1;
+ double **dynamicmtx = NULL;
+ if( specificityconsideration > 0.0 ) dynamicmtx = AllocateDoubleMtx( nalphabets, nalphabets );
+
+ if( alg == 'Y' || alg == 'r' ) ilim = njob - nadd;
+ else ilim = njob - 1;
for( i=0; i<ilim; i++ )
{
if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
- fprintf( stderr, "% 5d / %d\r", i, njob );
+ fprintf( stderr, "% 5d / %d\r", i, njob-nadd );
fflush( stderr );
- for( j=i+1; j<njob; j++ )
+
+ if( alg == 'Y' || alg == 'r' ) jst = njob - nadd;
+ else jst = i + 1;
+ for( j=jst; j<njob; j++ )
{
if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
{
- if( store_dist ) distancemtx[i][j] = 2.0;
- if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 2.0 );
+ if( store_dist )
+ {
+ if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
+ else distancemtx[i][j-i] = 3.0;
+ }
+ if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
continue;
}
- strcpy( aseq[i], seq[i] );
- strcpy( aseq[j], seq[j] );
- clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
- clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
+ strcpy( aseq[0], seq[i] );
+ strcpy( aseq[1], seq[j] );
+// clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
+// clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
+// fprintf( stderr, "Skipping conjuction..\n" );
+
+ effarr1[0] = 1.0;
+ effarr2[0] = 1.0;
+ mseq1[0] = aseq[0];
+ mseq2[0] = aseq[1];
+
+ thereisx = thereisxineachseq[i] + thereisxineachseq[j];
+// strcpy( distseq1[0], dseq[i] ); // nen no tame
+// strcpy( distseq2[0], dseq[j] ); // nen no tame
+ distseq1[0] = dseq[i];
+ distseq2[0] = dseq[j];
+
// fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
// fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
if( use_fft )
{
- pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
+ pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, 1, 1, alloclen, &intdum, NULL, 0, NULL );
// fprintf( stderr, "pscore (fft) = %f\n", pscore );
off1 = off2 = 0;
}
{
switch( alg )
{
- case( 'a' ):
- pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
- off1 = off2 = 0;
- break;
case( 't' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
+ pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
off1 = off2 = 0;
break;
case( 'A' ):
- pscore = G__align11( mseq1, mseq2, alloclen, outgap, outgap );
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+// if( store_localhom )
+ if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
+ {
+ pscore = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
+ if( thereisx ) pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+// dist = score2dist( L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ), selfscore[i], selfscore[j] ); // 2014/Feb/20
+ if( dist2offset( dist ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ G__align11( dynamicmtx, mseq1, mseq2, alloclen, outgap, outgap );
+ }
+// pscore = (double)naivepairscore11( *mseq1, *mseq2, 0.0 );
+ }
+#endif
+ }
+ else
+ pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
+ }
off1 = off2 = 0;
break;
case( 'N' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
- genL__align11( mseq1, mseq2, alloclen, &off1, &off2 );
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+ pscore = genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
+ if( thereisx )
+ {
+ strcpy( dumseq1[0], distseq1[0] );
+ strcpy( dumseq2[0], distseq2[0] );
+ pscore = genL__align11( n_dis_consweight_multi, dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
+ }
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+// fprintf( stderr, "dist = %f\n", score2dist( pscore, selfscore[i], selfscore[j] ) );
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+ if( dist2offset( dist ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ genL__align11( dynamicmtx, mseq1, mseq2, alloclen, &off1, &off2 );
+ }
+ }
+#endif
+ }
+ break;
+ case( 'R' ):
+ if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
+ pscore = 0.0;
+ else
+ pscore = (double)lastresx[i][j].score; // all pair
+ break;
+ case( 'r' ):
+ if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) )
+ pscore = (double)lastresx[i][j-(njob-nadd)].score;
+ else
+ pscore = 0.0;
+ break;
+ case( 'L' ):
+ if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
+ pscore = 0.0;
+ else
+ {
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+// if( store_localhom )
+ if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
+ {
+ pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 ); // all pair
+ if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // all pair
+#if 1
+ if( specificityconsideration > 0.0 )
+ {
+ dist = score2dist( pscore, selfscore[i], selfscore[j] );
+ if( ( scoreoffset = dist2offset( dist ) ) < 0.0 )
+ {
+ makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
+ strcpy( mseq1[0], seq[i] );
+ strcpy( mseq2[0], seq[j] );
+ L__align11( dynamicmtx, scoreoffset, mseq1, mseq2, alloclen, &off1, &off2 );
+ }
+ }
+#endif
+ }
+ else
+ pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // all pair
+ }
+ }
+// pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
break;
+ case( 'Y' ):
+ if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
+ {
+ if( usenaivescoreinsteadofalignmentscore )
+ {
+ L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
+ }
+ else
+ {
+ if( store_localhom )
+ {
+ pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
+ if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
+ }
+ else
+ pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
+ }
+ }
+ else
+ pscore = 0.0;
+ break;
+ case( 'a' ):
+ pscore = Aalign( mseq1, mseq2, effarr1, effarr2, 1, 1, alloclen );
+ off1 = off2 = 0;
+ break;
+#if 0
case( 'K' ):
pscore = genG__align11( mseq1, mseq2, alloclen );
off1 = off2 = 0;
break;
- case( 'L' ):
- pscore = G__align11_noalign( amino_dis, penalty, penalty_ex, mseq1, mseq2, alloclen );
- L__align11( mseq1, mseq2, alloclen, &off1, &off2 );
- break;
+#endif
case( 'H' ):
pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
break;
pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
off1 = off2 = 0;
break;
+ case( 'G' ):
+ pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
+ off1 = off2 = 0;
+ break;
case( 'M' ):
pscore = MSalign11( mseq1, mseq2, alloclen );
break;
#if SCOREOUT
fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
#endif
- if( !store_localhom )
- ;
- else if( alg == 'H' )
- putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
- else if( alg != 'S' && alg != 'V' )
- putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
-
-
- if( (bunbo=MIN( selfscore[i], selfscore[j] )) == 0.0 || bunbo < pscore )
- pscore = 2.0;
- else
- pscore = ( 1.0 - pscore / bunbo ) * 2.0;
+// if( pscore > 0.0 && ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) ) // x-ins-i de seido teika
+ if( ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) )
+ {
+ if( !store_localhom )
+ ;
+ else if( specifictarget && targetmap[i] == -1 && targetmap[j] == -1)
+ ;
+ else if( alg == 'R' )
+ putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j, 'h' );
+ else if( alg == 'r' )
+ putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd), 'h' );// ?????
+ else if( alg == 'H' )
+ putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else if( alg == 'Y' )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else if( !specifictarget && alg != 'S' && alg != 'V' )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-i, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else
+ {
+ if( targetmap[i] != -1 && targetmap[j] != -1 )
+ {
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' ); // sukoshi muda.
+ }
+ else if( targetmap[j] != -1 )
+ putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
+ else if( targetmap[i] != -1 )
+ putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
+ else
+ {
+ reporterr( "okashii\n" );
+ exit( 1 );
+ }
+ }
+ }
+
+ pscore = score2dist( pscore, selfscore[i], selfscore[j] );
}
else
{
}
}
if( stdout_dist ) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
- if( store_dist) distancemtx[i][j] = pscore;
+ if( store_dist)
+ {
+ if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
+ else distancemtx[i][j-i] = pscore;
+ }
}
}
+ if( dynamicmtx ) FreeDoubleMtx( dynamicmtx );
}
- if( store_dist )
+ if( store_dist && ngui == 0 )
{
hat2p = fopen( hat2file, "w" );
if( !hat2p ) ErrorExit( "Cannot open hat2." );
- WriteHat2_pointer( hat2p, njob, name, distancemtx );
+ if( alg == 'Y' || alg == 'r' )
+ WriteHat2_part_pointer( hat2p, njob, nadd, name, distancemtx );
+ else
+// WriteHat2_pointer( hat2p, njob, name, distancemtx );
+ WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, distancemtx ); // jissiha double
fclose( hat2p );
}
hat3p = fopen( "hat3", "w" );
if( !hat3p ) ErrorExit( "Cannot open hat3." );
- if( store_localhom )
+ if( store_localhom && ngui == 0 )
{
+
fprintf( stderr, "\n\n##### writing hat3\n" );
- ilim = njob-1;
+ if( alg == 'Y' || alg == 'r' )
+ ilim = njob-nadd;
+ else if( specifictarget )
+ ilim = ntarget;
+ else
+ ilim = njob-1;
for( i=0; i<ilim; i++ )
{
- for( j=i+1; j<njob; j++ )
+ if( alg == 'Y' || alg == 'r' )
+ {
+ jst = njob-nadd;
+ jj = 0;
+ }
+ else if( specifictarget )
{
- for( tmpptr=localhomtable[i]+j; tmpptr; tmpptr=tmpptr->next )
+ jst = 0;
+ jj = 0;
+ }
+ else
+ {
+ jst = i;
+ jj = 0;
+ }
+ for( j=jst; j<njob; j++, jj++ )
+ {
+ for( tmpptr=localhomtable[i]+jj; tmpptr; tmpptr=tmpptr->next )
{
+// fprintf( stderr, "j=%d, jj=%d\n", j, jj );
if( tmpptr->opt == -1.0 ) continue;
// tmptmptmptmptmp
// if( alg == 'B' || alg == 'T' )
// fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, 1.0, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
// else
- fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
+ if( targetmap[j] == -1 || targetmap[i] < targetmap[j] )
+ fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", targetmapr[i], j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
+// fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2+1, tmpptr->end2+1 ); // zettai dame!!!!
}
}
}
+// if( ngui == 0 )
+// {
#if DEBUG
- fprintf( stderr, "calling FreeLocalHomTable\n" );
+ fprintf( stderr, "calling FreeLocalHomTable\n" );
#endif
- FreeLocalHomTable( localhomtable, njob );
+ if( alg == 'Y' || alg == 'r' )
+ FreeLocalHomTable_part( localhomtable, (njob-nadd), nadd );
+ else if( specifictarget )
+ FreeLocalHomTable_part( localhomtable, ntarget, njob );
+ else
+ FreeLocalHomTable_half( localhomtable, njob );
#if DEBUG
- fprintf( stderr, "done. FreeLocalHomTable\n" );
+ fprintf( stderr, "done. FreeLocalHomTable\n" );
#endif
+// }
}
fclose( hat3p );
free( selfscore );
free( effarr1 );
free( effarr2 );
- free( indication1 );
- free( indication2 );
- if( store_dist ) FreeDoubleMtx( distancemtx );
-}
-
-static void WriteOptions( FILE *fp )
-{
-
- if( dorp == 'd' ) fprintf( fp, "DNA\n" );
- else
- {
- if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
- else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
- else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
- }
- fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
- if( use_fft ) fprintf( fp, "FFT on\n" );
-
- fprintf( fp, "tree-base method\n" );
- if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
- else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
- if( tbitr || tbweight )
+ if( alg == 'N' )
{
- fprintf( fp, "iterate at each step\n" );
- if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
- if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
- if( tbweight ) fprintf( fp, " weighted\n" );
- fprintf( fp, "\n" );
+ FreeCharMtx( dumseq1 );
+ FreeCharMtx( dumseq2 );
}
+ free( distseq1 );
+ free( distseq2 );
+ if( store_dist && ngui == 0 ) FreeDoubleHalfMtx( distancemtx, njob );
- fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
-
- if( alg == 'a' )
- fprintf( fp, "Algorithm A\n" );
- else if( alg == 'A' )
- fprintf( fp, "Algorithm A+\n" );
- else if( alg == 'S' )
- fprintf( fp, "Apgorithm S\n" );
- else
- fprintf( fp, "Unknown algorithm\n" );
-
- if( use_fft )
- {
- fprintf( fp, "FFT on\n" );
- if( dorp == 'd' )
- fprintf( fp, "Basis : 4 nucleotides\n" );
- else
- {
- if( fftscore )
- fprintf( fp, "Basis : Polarity and Volume\n" );
- else
- fprintf( fp, "Basis : 20 amino acids\n" );
- }
- fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
- fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
- }
- else
- fprintf( fp, "FFT off\n" );
- fflush( fp );
+ free( targetmap );
+ free( targetmapr );
}
-
-int main( int argc, char *argv[] )
+
+int pairlocalalign( int ngui, int lgui, char **namegui, char **seqgui, double **distancemtx, LocalHom **localhomtable, int argc, char **argv )
{
- int nlen[M];
+ int *nlen, *thereisxineachseq;
char **name, **seq;
char **mseq1, **mseq2;
char **aseq;
char **bseq;
- double *eff;
- int i;
+ char **dseq;
+ int i, j, k;
FILE *infp;
char c;
int alloclen;
+ Lastresx **lastresx;
+
+// reporterr( "argc=%d, argv[0]=%s\n", argc, argv[0] );
arguments( argc, argv );
-#ifndef enablemultithread
- nthread = 0;
-#endif
- if( inputfile )
+
+ if( !ngui )
{
- infp = fopen( inputfile, "r" );
- if( !infp )
+ if( inputfile )
+ {
+ infp = fopen( inputfile, "r" );
+ if( !infp )
+ {
+ fprintf( stderr, "Cannot open %s\n", inputfile );
+ exit( 1 );
+ }
+ }
+ else
+ infp = stdin;
+
+ getnumlen( infp );
+ rewind( infp );
+
+ if( njob < 2 )
+ {
+ fprintf( stderr, "At least 2 sequences should be input!\n"
+ "Only %d sequence found.\n", njob );
+ exit( 1 );
+ }
+ if( njob > M )
{
- fprintf( stderr, "Cannot open %s\n", inputfile );
+ fprintf( stderr, "The number of sequences must be < %d\n", M );
+ fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
exit( 1 );
}
}
- else
- infp = stdin;
-
- getnumlen( infp );
- rewind( infp );
- if( njob < 2 )
+ if( ( alg == 'r' || alg == 'R' ) && dorp == 'p' )
{
- fprintf( stderr, "At least 2 sequences should be input!\n"
- "Only %d sequence found.\n", njob );
+ fprintf( stderr, "Not yet supported\n" );
exit( 1 );
}
- if( njob > M )
+
+ alloclen = nlenmax*2;
+ if( ngui )
+ {
+ seq = seqgui;
+ name = namegui;
+ }
+ else
{
- fprintf( stderr, "The number of sequences must be < %d\n", M );
- fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
- exit( 1 );
+ seq = AllocateCharMtx( njob, alloclen+10 );
+ name = AllocateCharMtx( njob, B );
}
- alloclen = nlenmax*2;
- seq = AllocateCharMtx( njob, alloclen+10 );
- aseq = AllocateCharMtx( njob, alloclen+10 );
+ aseq = AllocateCharMtx( 2, alloclen+10 );
bseq = AllocateCharMtx( njob, alloclen+10 );
+ dseq = AllocateCharMtx( njob, alloclen+10 );
mseq1 = AllocateCharMtx( njob, 0 );
mseq2 = AllocateCharMtx( njob, 0 );
- name = AllocateCharMtx( njob, B );
+ nlen = AllocateIntVec( njob );
+ thereisxineachseq = AllocateIntVec( njob );
+
- eff = AllocateDoubleVec( njob );
+ if( alg == 'R' )
+ {
+ lastresx = calloc( njob+1, sizeof( Lastresx * ) );
+ for( i=0; i<njob; i++ )
+ {
+ lastresx[i] = calloc( njob+1, sizeof( Lastresx ) ); // muda
+ for( j=0; j<njob; j++ )
+ {
+ lastresx[i][j].score = 0;
+ lastresx[i][j].naln = 0;
+ lastresx[i][j].aln = NULL;
+ }
+ lastresx[i][njob].naln = -1;
+ }
+ lastresx[njob] = NULL;
+ }
+ else if( alg == 'r' )
+ {
+// fprintf( stderr, "Allocating lastresx (%d), njob=%d, nadd=%d\n", njob-nadd+1, njob, nadd );
+ lastresx = calloc( njob-nadd+1, sizeof( Lastresx * ) );
+ for( i=0; i<njob-nadd; i++ )
+ {
+// fprintf( stderr, "Allocating lastresx[%d]\n", i );
+ lastresx[i] = calloc( nadd+1, sizeof( Lastresx ) );
+ for( j=0; j<nadd; j++ )
+ {
+// fprintf( stderr, "Initializing lastresx[%d][%d]\n", i, j );
+ lastresx[i][j].score = 0;
+ lastresx[i][j].naln = 0;
+ lastresx[i][j].aln = NULL;
+ }
+ lastresx[i][nadd].naln = -1;
+ }
+ lastresx[njob-nadd] = NULL;
+ }
+ else
+ lastresx = NULL;
#if 0
Read( name, nlen, seq );
#else
- readData_pointer( infp, name, nlen, seq );
+ if( !ngui )
+ {
+ readData_pointer( infp, name, nlen, seq );
+ fclose( infp );
+ }
#endif
- fclose( infp );
constants( njob, seq );
initFiles();
- WriteOptions( trap_g );
+// WriteOptions( trap_g );
c = seqcheck( seq );
if( c )
// writePre( njob, name, nlen, seq, 0 );
- for( i=0; i<njob; i++ ) eff[i] = 1.0;
- for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
- pairalign( name, nlen, bseq, aseq, mseq1, mseq2, eff, alloclen );
+ for( i=0; i<njob; i++ )
+ {
+ gappick0( bseq[i], seq[i] );
+ thereisxineachseq[i] = removex( dseq[i], bseq[i] );
+ }
+
+ pairalign( name, nlen, bseq, aseq, dseq, thereisxineachseq, mseq1, mseq2, alloclen, lastresx, distancemtx, localhomtable, ngui );
fprintf( trap_g, "done.\n" );
#if DEBUG
fprintf( stderr, "closing trap_g\n" );
#endif
fclose( trap_g );
+ fclose( prep_g );
// writePre( njob, name, nlen, aseq, !contin );
#if 0
{
fprintf( stderr, "\nThe order of pairwise alignments is not identical to that in the input file, because of the parallel calculation. Reorder them by yourself.\n" );
}
- FreeCharMtx( seq );
+
+#if 1
+ if( lastresx )
+ {
+ for( i=0; lastresx[i]; i++ )
+ {
+ for( j=0; lastresx[i][j].naln!=-1; j++ )
+ {
+ for( k=0; k<lastresx[i][j].naln; k++ )
+ {
+ free( lastresx[i][j].aln[k].reg1 );
+ free( lastresx[i][j].aln[k].reg2 );
+ }
+ free( lastresx[i][j].aln );
+ }
+ free( lastresx[i] );
+ }
+ free( lastresx );
+ }
+#endif
+ if( ngui == 0 )
+ {
+ FreeCharMtx( seq );
+ FreeCharMtx( name );
+ }
FreeCharMtx( aseq );
FreeCharMtx( bseq );
- FreeCharMtx( name );
+ FreeCharMtx( dseq );
free( mseq1 );
free( mseq2 );
- free( eff );
+ free( nlen );
+ free( thereisxineachseq );
+ freeconstants();
+
+ if( !ngui )
+ {
+ FreeCommonIP();
+ }
+ Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
+ G__align11( NULL, NULL, NULL, 0, 0, 0 ); // 20130603
+ G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
+ L__align11( NULL, 0.0, NULL, NULL, 0, NULL, NULL );
+ L__align11_noalign( NULL, NULL, NULL );
+ genL__align11( NULL, NULL, NULL, 0, NULL, NULL );
+
+#if SHISHAGONYU
+ if( ngui )
+ {
+ char buf[100];
+ for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
+ {
+ sprintf( buf, "%5.3f", distancemtx[i][j-i] );
+ distancemtx[i][j-i] = 0.0;
+ sscanf( buf, "%lf", distancemtx[i]+j-i );
+// distancemtx[i][j-i] = 0.001 * (int)(distancemtx[i][j-i] * 1000 + 0.5);
+ }
+
+ }
+#endif
+
return( 0 );
}
+
git://source.jalview.org / jabaws.git / blobdiff