palign.c

在任务级并行平台P2HP上

/* Change int h to int gh everywhere  DES June 1994 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "clustalw.h"
#include <jni.h>

#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))

#define gap(k)  ((k) <= 0 ? 0 : g + gh * (k))
#define tbgap(k)  ((k) <= 0 ? 0 : tb + gh * (k))
#define tegap(k)  ((k) <= 0 ? 0 : te + gh * (k))

/*
*	Prototypes
*/
static void add(sint v);
static sint calc_score(sint iat, sint jat, sint v1, sint v2);
static float tracepath(sint tsb1,sint tsb2);
static void forward_pass(char *ia, char *ib, sint n, sint m);
static void reverse_pass(char *ia, char *ib);
static sint diff(sint A, sint B, sint M, sint N, sint tb, sint te);
static void del(sint k);


static void  read_paramfile(char * p);
static void writetofile(char * filename, int a, int b);

/*
*   Global variables
*/
#ifdef MAC
#define pwint   short
#else
#define pwint   int
#endif
static sint		int_scale;

extern double   **tmat;
extern float    pw_go_penalty;
extern float    pw_ge_penalty;
extern float	transition_weight;
extern sint 	nseqs;
extern sint 	max_aa;
extern sint 	gap_pos1,gap_pos2;
extern sint  	max_aln_length;
extern sint 	*seqlen_array;
extern sint 	debug;
extern sint  	mat_avscore;
extern short 	blosum30mt[],pam350mt[],idmat[],pw_usermat[],pw_userdnamat[];
extern short    clustalvdnamt[],swgapdnamt[];
extern short    gon250mt[];
extern short 	def_dna_xref[],def_aa_xref[],pw_dna_xref[],pw_aa_xref[];
extern Boolean  dnaflag;
extern char 	**seq_array;
extern char 	*amino_acid_codes;
extern char 	pw_mtrxname[];
extern char 	pw_dnamtrxname[];

static float 	mm_score;
static sint 	print_ptr,last_print;
static sint 	*displ;
static pwint 	*HH, *DD, *RR, *SS;
static sint 	g, gh;
static sint   	seq1, seq2;
static sint     matrix[NUMRES][NUMRES];
static pwint    maxscore;
static sint    	sb1, sb2, se1, se2;

//added by jf.
static sint       Asequence,Bsequence,Esequence,Astart,Bend;
static sint       align_model;

/******************************
Modified by jf
2005.1.18
*******************************/

//sint pairalign(sint istart, sint jstart, sint jend)
sint pairalign(sint istart, sint iend, sint jstart, sint jend,char *datapath,char *resultpath,char *paramfilename)
{
	short	 *mat_xref;
	static sint    si, sj, i;
	static sint    n,m,len1,len2;
	static sint    maxres;                          
	static short    *matptr;
	static char   c;
	static float gscale,ghscale;
	
	/***
	added by jf.
	2005.3.15.
	***/
	char paramfile[FILENAMELEN+1];
	char resultfile[FILENAMELEN+1];

	Boolean wfFlag=TRUE;
	
	displ = (sint *)ckalloc((2*max_aln_length+1) * sizeof(sint));
	HH = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
	DD = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
	RR = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
	SS = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
	
#ifdef MAC
	int_scale = 10;
#else
	int_scale = 100;
#endif
	gscale=ghscale=1.0;
	if (dnaflag)
    {
		if (debug>1) fprintf(stdout,"matrix %s\n",pw_dnamtrxname);
		if (strcmp(pw_dnamtrxname, "iub") == 0)
		{ 
			matptr = swgapdnamt;
			mat_xref = def_dna_xref;
		}
		else if (strcmp(pw_dnamtrxname, "clustalw") == 0)
		{ 
			matptr = clustalvdnamt;
			mat_xref = def_dna_xref;
			gscale=0.6667;
			ghscale=0.751;
		}
		else
		{
			matptr = pw_userdnamat;
			mat_xref = pw_dna_xref;
		}
		maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
		if (maxres == 0) return((sint)-1);
		
		matrix[0][4]=transition_weight*matrix[0][0];
		matrix[4][0]=transition_weight*matrix[0][0];
		matrix[2][11]=transition_weight*matrix[0][0];
		matrix[11][2]=transition_weight*matrix[0][0];
		matrix[2][12]=transition_weight*matrix[0][0];
		matrix[12][2]=transition_weight*matrix[0][0];
    }
	else
    {
		if (debug>1) fprintf(stdout,"matrix %s\n",pw_mtrxname);
		if (strcmp(pw_mtrxname, "blosum") == 0)
		{
			matptr = blosum30mt;
			mat_xref = def_aa_xref;
		}
		else if (strcmp(pw_mtrxname, "pam") == 0)
		{
			matptr = pam350mt;
			mat_xref = def_aa_xref;
		}
		else if (strcmp(pw_mtrxname, "gonnet") == 0)
		{
			matptr = gon250mt;
			int_scale /= 10;
			mat_xref = def_aa_xref;
		}
		else if (strcmp(pw_mtrxname, "id") == 0)
		{
			matptr = idmat;
			mat_xref = def_aa_xref;
		}
		else
		{
			matptr = pw_usermat;
			mat_xref = pw_aa_xref;
		}
		
		maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
		if (maxres == 0) return((sint)-1);
    }
	
	
	
	
	//begin modify
	strcpy(paramfile,datapath);
	strcat(paramfile,paramfilename);
	
	read_paramfile(paramfile);
	
	strcpy(resultfile,resultpath); 
	align_model==0;
	//added by jf. 
	//2005.3.11
	if(align_model==1)
	{
		si=MAX(0,istart); 
		//***for (si=MAX(0,istart);si<nseqs && si<iend;si++)
		// {
		n = seqlen_array[si+1];
		len1 = 0;
		for (i=1;i<=n;i++) {
			c = seq_array[si+1][i];
			if ((c!=gap_pos1) && (c != gap_pos2)) len1++;
		}
		
		for (sj=MAX(si+1,jstart+1);sj<nseqs && sj<jend;sj++)
		{
			m = seqlen_array[sj+1];
			
			if(n==0 || m==0) {
			/***
			modify the matrix tmat
				***/
				tmat[si+1][sj+1]=1.0;
				tmat[sj+1][si+1]=1.0;
				continue;
			}
			len2 = 0;
			for (i=1;i<=m;i++) {
				c = seq_array[sj+1][i];
				if ((c!=gap_pos1) && (c != gap_pos2)) len2++;
			}
			
			if (dnaflag) {
				g = 2 * (float)pw_go_penalty * int_scale*gscale;
				gh = pw_ge_penalty * int_scale*ghscale;
			}
			else {
				if (mat_avscore <= 0)
					g = 2 * (float)(pw_go_penalty + log((double)(MIN(n,m))))*int_scale;
				else
					g = 2 * mat_avscore * (float)(pw_go_penalty +
					log((double)(MIN(n,m))))*gscale;
				gh = pw_ge_penalty * int_scale;
			}
			
			if (debug>1) fprintf(stdout,"go %d ge %d\n",(pint)g,(pint)gh);
			
			/*
			align the sequences
			*/
			seq1 = si+1;
			seq2 = sj+1;
			
			forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
				n, m);
			
			reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0]);
			
			last_print = 0;
			print_ptr = 1;
			/*
			sb1 = sb2 = 1;
			se1 = n-1;
			se2 = m-1;
			*/
			
			/* use Myers and Miller to align two sequences */
			
			maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
				(sint)0, (sint)0);
			
			/* calculate percentage residue identity */
			
			mm_score = tracepath(sb1,sb2);
			
			if(len1==0 || len2==0) mm_score=0;
			else
				mm_score /= (float)MIN(len1,len2);
			
			tmat[si+1][sj+1] = ((float)100.0 - mm_score)/(float)100.0;
			tmat[sj+1][si+1] = ((float)100.0 - mm_score)/(float)100.0;
			
			if (debug>1)
			{
				fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
					(pint)si+1,(pint)sj+1, 
					(pint)mm_score, 
					(pint)maxscore/(MIN(len1,len2)*100));
			}
			else
			{
				info("Sequences (%d:%d) Aligned. Score:  %d",
					(pint)si+1,(pint)sj+1, 
					(pint)mm_score);
			}
			writetofile(wfFlag,resultfile,si,sj, Bsequence+si, Bsequence+sj);
			wfFlag=FALSE;
		}
	}
	else if(align_model==0)
	{
		for (si=MAX(0,Asequence);si<nseqs && si<iend && si<=Bsequence;si++)
		{
			n = seqlen_array[si+1];
			len1 = 0;
			for (i=1;i<=n;i++) {
				c = seq_array[si+1][i];
				if ((c!=gap_pos1) && (c != gap_pos2)) len1++;
			}
			if (si==Asequence) {
				sj=MAX(si+1,jstart+1);
				sj=MAX(sj,Astart);
			}
			else sj=MAX(si+1,jstart+1);
			for (;sj<nseqs && sj<jend;sj++)
			{
				if (si==Bsequence&&sj>Bend) {
					break;
				}
				m = seqlen_array[sj+1];
				
				if(n==0 || m==0) {
					/***
						modify the matrix tmat
						***/
						tmat[si+1][sj+1]=1.0;
					tmat[sj+1][si+1]=1.0;
					continue;
				}
				len2 = 0;
				for (i=1;i<=m;i++) {
					c = seq_array[sj+1][i];
					if ((c!=gap_pos1) && (c != gap_pos2)) len2++;
				}
				
				if (dnaflag) {
					g = 2 * (float)pw_go_penalty * int_scale*gscale;
					gh = pw_ge_penalty * int_scale*ghscale;
				}
				else {
					if (mat_avscore <= 0)
						g = 2 * (float)(pw_go_penalty + log((double)(MIN(n,m))))*int_scale;
					else
						g = 2 * mat_avscore * (float)(pw_go_penalty +
						log((double)(MIN(n,m))))*gscale;
					gh = pw_ge_penalty * int_scale;
				}
				
				if (debug>1) fprintf(stdout,"go %d ge %d\n",(pint)g,(pint)gh);
				
				/*align the sequences*/
					seq1 = si+1;
				seq2 = sj+1;
				
				forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
					n, m);
				
				reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0]);
				
				last_print = 0;
				print_ptr = 1;
				/*
					sb1 = sb2 = 1;
				se1 = n-1;
				se2 = m-1;
				*/
					
					/* use Myers and Miller to align two sequences */
					
					maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
					(sint)0, (sint)0);
				
				/* calculate percentage residue identity */
					
					mm_score = tracepath(sb1,sb2);
				
				if(len1==0 || len2==0) mm_score=0;
				else
					mm_score /= (float)MIN(len1,len2);
				
				tmat[si+1][sj+1] = ((float)100.0 - mm_score)/(float)100.0;
				tmat[sj+1][si+1] = ((float)100.0 - mm_score)/(float)100.0;
				
				if (debug>1)
				{
					fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
						(pint)si+1,(pint)sj+1, 
						(pint)mm_score, 
						(pint)maxscore/(MIN(len1,len2)*100));
				}
				else
				{
					info("Sequences (%d:%d) Aligned. Score:  %d",
						(pint)si+1,(pint)sj+1, 
						(pint)mm_score);
				}
				//writetofile(resultfile,si,sj, Bsequence+si+1,Bsequence+sj+1);
				writetofile(wfFlag,resultfile,si,sj, si+1, sj+1);
				wfFlag=FALSE;
			}
		}
	}
	/*
	{
			for (si=MAX(0,istart);si<nseqs && si<iend;si++)
			{
				n = seqlen_array[si+1];
				len1 = 0;
				for (i=1;i<=n;i++) {
					c = seq_array[si+1][i];
					if ((c!=gap_pos1) && (c != gap_pos2)) len1++;
				}
				
				for (sj=MAX(si+1,jstart+1);sj<nseqs && sj<jend;sj++)
				{
					m = seqlen_array[sj+1];
					
					if(n==0 || m==0) {
					/ ***
					modify the matrix tmat
						*** /
						tmat[si+1][sj+1]=1.0;
						tmat[sj+1][si+1]=1.0;
						continue;
					}
					len2 = 0;
					for (i=1;i<=m;i++) {
						c = seq_array[sj+1][i];
						if ((c!=gap_pos1) && (c != gap_pos2)) len2++;
					}
					
					if (dnaflag) {
						g = 2 * (float)pw_go_penalty * int_scale*gscale;
						gh = pw_ge_penalty * int_scale*ghscale;
					}
					else {
						if (mat_avscore <= 0)
							g = 2 * (float)(pw_go_penalty + log((double)(MIN(n,m))))*int_scale;
						else
							g = 2 * mat_avscore * (float)(pw_go_penalty +
							log((double)(MIN(n,m))))*gscale;
						gh = pw_ge_penalty * int_scale;
					}
					
					if (debug>1) fprintf(stdout,"go %d ge %d\n",(pint)g,(pint)gh);
					
					/ *
					align the sequences
					* /
					seq1 = si+1;
					seq2 = sj+1;
					
					forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
						n, m);
					
					reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0]);
					
					last_print = 0;
					print_ptr = 1;
					/ *
					sb1 = sb2 = 1;
					se1 = n-1;
					se2 = m-1;
					* /
					
					/ * use Myers and Miller to align two sequences * /
					
					maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
						(sint)0, (sint)0);
					
					/ * calculate percentage residue identity * /
					
					mm_score = tracepath(sb1,sb2);
					
					if(len1==0 || len2==0) mm_score=0;
					else
						mm_score /= (float)MIN(len1,len2);
					
					tmat[si+1][sj+1] = ((float)100.0 - mm_score)/(float)100.0;
					tmat[sj+1][si+1] = ((float)100.0 - mm_score)/(float)100.0;
					
					if (debug>1)
					{
						fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
							(pint)si+1,(pint)sj+1, 
							(pint)mm_score, 
							(pint)maxscore/(MIN(len1,len2)*100));
					}
					else
					{
						info("Sequences (%d:%d) Aligned. Score:  %d",
							(pint)si+1,(pint)sj+1, 
							(pint)mm_score);
					}
					writetofile(resultfile,si,sj, Bsequence+si+1,Bsequence+sj+1);
				}
			}
		}*/
	
	displ=ckfree((void *)displ);
	HH=ckfree((void *)HH);
	DD=ckfree((void *)DD);
	RR=ckfree((void *)RR);
	SS=ckfree((void *)SS);
	
	
	return((sint)1);
}

/***
added by jf.