sequence.c

clustalx用来做基因序列分析

	for (;;) {
		while (fgets(line,MAXLINE+1,fin) != NULL)
                	if( *line == '{' ) break;

		while( !keyword(line,"name") )
			if (fgets(line,MAXLINE+1,fin) == NULL) return;
			
    	for(i=5;i<=strlen(line);i++)  /* DES */
			if(line[i] != ' ') break;
		strncpy(sname,line+i,MAXNAMES); /* remember entryname */
    		for(i=0;i<=strlen(sname);i++)
			if(sname[i] == ' ') {
				sname[i]=EOS;
				break;
			}
		sname[MAXNAMES]=EOS;
		rtrim(sname);
    	blank_to_(sname);
		
/* look for secondary structure feature table / gap penalty mask */
		while(fgets(line,MAXLINE+1,fin) != NULL) {
			if (keyword(line,"feature")) {
				if (interactive) {
					strcpy(title,"Found secondary structure in alignment file: ");
					strcat(title,sname);
					(*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
				}
				else (*lin2) = 'y';
				if ((*lin2 != 'n') && (*lin2 != 'N'))  {               	
					struct_penalties = SECST;
					for (i=0;i<length;i++)
						sec_struct_mask[i] = '.';
					do {
						if(keyword(line,"feature"))
							get_rsf_feature(&line[7],length);
						fgets(line,MAXLINE+1,fin);
					} while( !keyword(line,"sequence") );
				}
				else {
					do {
						fgets(line,MAXLINE+1,fin);
					} while( !keyword(line,"sequence") );
				}
				strcpy(ss_name,sname);
			}
			else if (keyword(line,"sequence"))
				break;	

			if (struct_penalties != NONE) break;			
		}
						
	}
						
}

static void get_gde_ss(sint length)
{
	static char title[MAXLINE+1];
	static char line[MAXLINE+1];
	static char lin2[MAXLINE+1];
	static char sname[MAXNAMES+1];
	sint i, len, offset = 0;
        unsigned char c;

	for (;;) {
		line[0] = '\0';
/* search for the next comment line */
		while(*line != '"')
			if (fgets(line,MAXLINE+1,fin) == NULL) return;

/* is it a secondary structure entry? */
		if (strncmp(&line[1],"SS_",3) == 0) {
			for (i=1;i<=MAXNAMES-3;i++) {
				if (line[i+3] == '(' || line[i+3] == '\n')
						break;
				sname[i-1] = line[i+3];
			}
			i--;
			sname[i]=EOS;
			if (sname[i-1] == '(') sscanf(&line[i+3],"%d",&offset);
			else offset = 0;
			for(i--;i > 0;i--) 
				if(isspace(sname[i])) {
					sname[i]=EOS;	
				}
				else break;		
			blank_to_(sname);

			if (interactive) {
				strcpy(title,"Found secondary structure in alignment file: ");
				strcat(title,sname);
				(*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
			}
			else (*lin2) = 'y';
			if ((*lin2 != 'n') && (*lin2 != 'N'))  {               	
				struct_penalties = SECST;
				for (i=0;i<length;i++)
					sec_struct_mask[i] = '.';
				len = 0;
				while(fgets(line,MAXLINE+1,fin)) {
					if(*line == '%' || *line == '#' || *line == '"') break;
					for(i=offset;i < length;i++) {
						c=line[i];
						if(c == '\n' || c == EOS) 
							break;			/* EOL */
						sec_struct_mask[len++]=c;
					}
					if (len > length) break;
				}
				strcpy(ss_name,sname);
			}
		}
/* or is it a gap penalty mask entry? */
		else if (strncmp(&line[1],"GM_",3) == 0) {
			for (i=1;i<=MAXNAMES-3;i++) {
				if (line[i+3] == '(' || line[i+3] == '\n')
						break;
				sname[i-1] = line[i+3];
			}
			i--;
			sname[i]=EOS;
			if (sname[i-1] == '(') sscanf(&line[i+3],"%d",&offset);
			else offset = 0;
			for(i--;i > 0;i--) 
				if(isspace(sname[i])) {
					sname[i]=EOS;	
				}
				else break;		
			blank_to_(sname);

			if (interactive) {
				strcpy(title,"Found gap penalty mask in alignment file: ");
				strcat(title,sname);
				(*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
			}
			else (*lin2) = 'y';
			if ((*lin2 != 'n') && (*lin2 != 'N'))  {               	
				struct_penalties = GMASK;
				for (i=0;i<length;i++)
					gap_penalty_mask[i] = '1';
				len = 0;
				while(fgets(line,MAXLINE+1,fin)) {
					if(*line == '%' || *line == '#' || *line == '"') break;
					for(i=offset;i < length;i++) {
						c=line[i];
						if(c == '\n' || c == EOS) 
							break;			/* EOL */
						gap_penalty_mask[len++]=c;
					}
					if (len > length) break;
				}
				strcpy(ss_name,sname);
			}
		}
		if (struct_penalties != NONE) break;			
	}			
			
}

static void get_swiss_feature(char *line, sint len)
{
	char c, s, feature[MAXLINE+1];
	int  i, start_pos, end_pos;
	
	if (sscanf(line,"%s%d%d",feature,&start_pos,&end_pos) != 3) {
		return;
	}

	if (strcmp(feature,"HELIX") == 0) {
		c = 'A';
		s = '$';
	}
	else if (strcmp(feature,"STRAND") == 0) {
		c = 'B';
		s = '%';
	}
	else
		return;
			
	if(start_pos >=len || end_pos>=len) return;

	sec_struct_mask[start_pos-1] = s;
	for (i=start_pos;i<end_pos-1;i++)
		sec_struct_mask[i] = c;
	sec_struct_mask[end_pos-1] = s;
		
}

static void get_rsf_feature(char *line, sint len)
{
	char c, s;
	char str1[MAXLINE+1],str2[MAXLINE+1],feature[MAXLINE+1];
	int  i, tmp,start_pos, end_pos;
	
	if (sscanf(line,"%d%d%d%s%s%s",&start_pos,&end_pos,&tmp,str1,str2,feature) != 6) {
		return;
	}

	if (strcmp(feature,"HELIX") == 0) {
		c = 'A';
		s = '$';
	}
	else if (strcmp(feature,"STRAND") == 0) {
		c = 'B';
		s = '%';
	}
	else
		return;
			
	if(start_pos>=len || end_pos >= len) return;
	sec_struct_mask[start_pos-1] = s;
	for (i=start_pos;i<end_pos-1;i++)
		sec_struct_mask[i] = c;
	sec_struct_mask[end_pos-1] = s;
		
}

static void get_swiss_mask(char *line, sint len)
{
	int  i, value, start_pos, end_pos;
	
	if (sscanf(line,"%d%d%d",&value,&start_pos,&end_pos) != 3) {
		return;
	}

	if (value < 1 || value > 9) return;
	
	if(start_pos>=len || end_pos >= len) return;
	for (i=start_pos-1;i<end_pos;i++)
		gap_penalty_mask[i] = value+'0';
		
}

static char * get_seq(char *sname,sint *len,char *tit)
{
	static char line[MAXLINE+1];
	char *seq = NULL;
	sint i, offset = 0;
        unsigned char c=EOS;
	Boolean got_seq=FALSE;

	switch(seqFormat) {

/************************************/
		case EMBLSWISS:
			while( !linetype(line,"ID") )
				if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
			
                        for(i=5;i<=strlen(line);i++)  /* DES */
				if(line[i] != ' ') break;
			strncpy(sname,line+i,MAXNAMES); /* remember entryname */
                	for(i=0;i<=strlen(sname);i++)
                        	if(sname[i] == ' ') {
                                	sname[i]=EOS;
                                	break;
                        	}

			sname[MAXNAMES]=EOS;
			rtrim(sname);
                        blank_to_(sname);

						
			while( !linetype(line,"SQ") )
				fgets(line,MAXLINE+1,fin);
				
			*len=0;
			while(fgets(line,MAXLINE+1,fin)) {
				if(got_seq && blankline(line)) break;
 				if( strlen(line) > 2 && line[strlen(line)-2]=='.' && line[strlen(line)-3]=='.' ) 
					continue;
				if(seq==NULL)
					seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
				else
					seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
				for(i=0;i<=MAXLINE;i++) {
					c=line[i];
				if(c == '\n' || c == EOS || c == '/')
					break;			/* EOL */
				c=chartab[c];
				if(c) {
					got_seq=TRUE;
					seq[++(*len)]=c;
				}
				}
			if(c == '/') break;
			}
		break;
		
/************************************/
		case PIR:
			while(*line != '>')
				fgets(line,MAXLINE+1,fin);			
                        for(i=4;i<=strlen(line);i++)  /* DES */
				if(line[i] != ' ') break;
			strncpy(sname,line+i,MAXNAMES); /* remember entryname */
			sname[MAXNAMES]=EOS;
			rtrim(sname);
                        blank_to_(sname);

			fgets(line,MAXLINE+1,fin);
			strncpy(tit,line,MAXTITLES);
			tit[MAXTITLES]=EOS;
			i=strlen(tit);
			if(tit[i-1]=='\n') tit[i-1]=EOS;
			
			*len=0;
			while(fgets(line,MAXLINE+1,fin)) {
				if(seq==NULL)
					seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
				else
					seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
				for(i=0;i<=MAXLINE;i++) {
					c=line[i];
				if(c == '\n' || c == EOS || c == '*')
					break;			/* EOL */
			
				c=chartab[c];
				if(c) seq[++(*len)]=c;
				}
			if(c == '*') break;
			}
		break;
/***********************************************/
		case PEARSON:
			while(*line != '>')
				fgets(line,MAXLINE+1,fin);
			
                        for(i=1;i<=strlen(line);i++)  /* DES */
				if(line[i] != ' ') break;
			strncpy(sname,line+i,MAXNAMES); /* remember entryname */
                        for(i=1;i<=strlen(sname);i++)  /* DES */
				if(sname[i] == ' ') break;
			sname[i]=EOS;
			rtrim(sname);
                        blank_to_(sname);

			*tit=EOS;
			
			*len=0;
			while(fgets(line,MAXLINE+1,fin)) {
				if(seq==NULL)
					seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
				else
					seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
				for(i=0;i<=MAXLINE;i++) {
					c=line[i];
				if(c == '\n' || c == EOS || c == '>')
					break;			/* EOL */
			
				c=chartab[c];
				if(c) seq[++(*len)]=c;
			}
			if(c == '>') break;
			}
		break;
/**********************************************/
		case GDE:
			if (dnaflag) {
				while(*line != '#')
					fgets(line,MAXLINE+1,fin);
			}
			else {
				while(*line != '%')
					fgets(line,MAXLINE+1,fin);
			}
			
			for (i=1;i<=MAXNAMES;i++) {
				if (line[i] == '(' || line[i] == '\n')
                                    break;
				sname[i-1] = line[i];
			}
			i--;
			sname[i]=EOS;
			if (sname[i-1] == '(') sscanf(&line[i],"%d",&offset);
			else offset = 0;
			for(i--;i > 0;i--) 
				if(isspace(sname[i])) {
					sname[i]=EOS;	
				}
				else break;		
                        blank_to_(sname);

			*tit=EOS;
			
			*len=0;
			for (i=0;i<offset;i++) seq[++(*len)] = '-';
			while(fgets(line,MAXLINE+1,fin)) {
			if(*line == '%' || *line == '#' || *line == '"') break;
				if(seq==NULL)
					seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
				else
					seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
				for(i=0;i<=MAXLINE;i++) {
					c=line[i];
				if(c == '\n' || c == EOS) 
					break;			/* EOL */
			
				c=chartab[c];
				if(c) seq[++(*len)]=c;
				}
			}
		break;
/***********************************************/
		case RSF:
			while(*line != '{')
				if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
			
			while( !keyword(line,"name") )
				if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
			
                        for(i=5;i<=strlen(line);i++)  /* DES */
				if(line[i] != ' ') break;
			strncpy(sname,line+i,MAXNAMES); /* remember entryname */
                	for(i=0;i<=strlen(sname);i++)
                        	if(sname[i] == ' ') {
                                	sname[i]=EOS;
                                	break;
                        	}

			sname[MAXNAMES]=EOS;
			rtrim(sname);
                        blank_to_(sname);

						
			while( !keyword(line,"sequence") )
				if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
				
			*len=0;
			while(fgets(line,MAXLINE+1,fin)) {
				if(seq==NULL)
					seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
				else
					seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
				for(i=0;i<=MAXLINE;i++) {
					c=line[i];
					if(c == EOS || c == '}')
						break;			/* EOL */
					if( c=='.')
						seq[++(*len)]='-';
					c=chartab[c];
					if(c) seq[++(*len)]=c;
				}
				if(c == '}') break;
			}
		break;
/***********************************************/
	}
	
	seq[*len+1]=EOS;

	return seq;
}


sint readseqs(sint first_seq) /*first_seq is the #no. of the first seq. to read */
{
	char line[FILENAMELEN+1];
	char fileName[FILENAMELEN+1];

	static char *seq1,sname1[MAXNAMES+1],title[MAXTITLES+1];