modseqalign.c

马尔科夫模型的java版本实现

    free((*(msa_seq_infop_template)).msa_seq_1);
    if(hmm.nr_alphabets > 1) {
      free((*(msa_seq_infop_template)).msa_seq_2);
    }
    if(hmm.nr_alphabets > 2) {
      free((*(msa_seq_infop_template)).msa_seq_3);
    }
    if(hmm.nr_alphabets > 3) {
      free((*(msa_seq_infop_template)).msa_seq_4);
    }
    free((*(msa_seq_infop_template)).gaps);
    free((*(msa_seq_infop_template)).gap_shares);
    free((*(msa_seq_infop_template)).lead_columns_start);
    free(msa_seq_infop_template);
    
    free((*(msa_seq_infop_target)).msa_seq_1);
    if(hmm.nr_alphabets > 1) {
      free((*(msa_seq_infop_target)).msa_seq_2);
    }
    if(hmm.nr_alphabets > 2) {
      free((*(msa_seq_infop_target)).msa_seq_3);
    }
    if(hmm.nr_alphabets > 3) {
      free((*(msa_seq_infop_target)).msa_seq_4);
    }
    free((*(msa_seq_infop_target)).gaps);
    free((*(msa_seq_infop_target)).gap_shares);
    free((*(msa_seq_infop_target)).lead_columns_start);
    free(msa_seq_infop_target);
  }
  else if(seq_format == FASTA || seq_format == STANDARD) {
    free(((seq_info_template.seqs))->seq_1);
    if(hmm.nr_alphabets > 1) {
      free((seq_info_template.seqs)->seq_2);
    }
    if(hmm.nr_alphabets > 2) {
      free((seq_info_template.seqs)->seq_3);
      }
    if(hmm.nr_alphabets > 3) {
      free((seq_info_template.seqs)->seq_4);
    }
    free(seq_info_template.seqs);
    
    free(((seq_info_target.seqs))->seq_1);
    if(hmm.nr_alphabets > 1) {
      free((seq_info_target.seqs)->seq_2);
    }
    if(hmm.nr_alphabets > 2) {
      free((seq_info_target.seqs)->seq_3);
      }
    if(hmm.nr_alphabets > 3) {
      free((seq_info_target.seqs)->seq_4);
    }
    free(seq_info_target.seqs);
  }


  if(substmtxfile != NULL) {
    free(subst_mtxp);
    if(hmm.nr_alphabets > 1) {
      free(subst_mtxp_2);
    }
    if(hmm.nr_alphabets > 2) {
      free(subst_mtxp_3);
    }
    if(hmm.nr_alphabets > 3) {
      free(subst_mtxp_4);
    }
    fclose(substmtxfile);
  }
  
  if(freqfile != NULL) {
    free(aa_freqs);
    if(hmm.nr_alphabets > 1) {
      free(aa_freqs_2);
    }
    if(hmm.nr_alphabets > 2) {
      free(aa_freqs_3);
    }
    if(hmm.nr_alphabets > 3) {
      free(aa_freqs_4);
    }
    fclose(freqfile);
  }
  hmm_garbage_collection_multi(hmmfile, &hmm);


  free(path_template);
  free(path_target);
 
}

void align_paths(int **path_templatep, int **path_targetp, int seq_len_target, int seq_len_template, FILE *outfile,
		 int seq_format)
{
  int *path_template, *path_target;
  struct align_mtx_element_s *align_mtx;
  int i,j;
  int score_up, score_left, score_left_up, score;
  char last_up, last_left, last_left_up, last;
  int alignment_length, cur_alignment_pos;
  int target_pos, template_pos;
  struct alignment_s *alignment;
  

  alignment_length = seq_len_target + seq_len_template;

  path_template = *path_templatep;
  path_target = *path_targetp;
  
  align_mtx = (struct align_mtx_element_s*)(malloc_or_die((seq_len_target+1) * (seq_len_template+1) *
							  sizeof(struct align_mtx_element_s)));
  
  (align_mtx + get_mtx_index(0,0,seq_len_target+1))->score = 0;
  (align_mtx + get_mtx_index(0,0,seq_len_target+1))->last = 'e';
  for(i = 1; i <= seq_len_target; i++) {
    (align_mtx + get_mtx_index(0,i,seq_len_target+1))->score = 0;
    (align_mtx + get_mtx_index(0,i,seq_len_target+1))->last = 'l';
  }

  for(j = 1; j <= seq_len_template; j++) {
    (align_mtx + get_mtx_index(j,0,seq_len_target+1))->score = 0;
    (align_mtx + get_mtx_index(j,0,seq_len_target+1))->last = 'u';
  }
  

  for(j = 1; j <= seq_len_template; j++) {
    for(i = 1; i <= seq_len_target; i++) {
      score_up = (align_mtx + get_mtx_index(j-1,i,seq_len_target+1))->score;
      score_left = (align_mtx + get_mtx_index(j,i-1,seq_len_target+1))->score;
      score_left_up = (align_mtx + get_mtx_index(j-1,i-1,seq_len_target+1))->score;
      last_up = (align_mtx + get_mtx_index(j-1,i,seq_len_target+1))->last;
      last_left = (align_mtx + get_mtx_index(j,i-1,seq_len_target+1))->last;
      last_left_up = (align_mtx + get_mtx_index(j-1,i-1,seq_len_target+1))->last;
      if(*(path_target + i) == *(path_template + j)) {
	score_left_up += 1;
      }

      if(score_left_up > score_left && score_left_up > score_up) {
	/* left_up */
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
      }
      else if(score_up > score_left && score_up > score_left_up) {
	/* up */
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
      }
      else if(score_left > score_left_up && score_left > score_up) {
	/* left */
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	(align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
      }
      else if(score_left_up > score_left) {
	/* left_up + up */
	if(last_left_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}
	else if(last_left == -1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
	}
	else if(last_left == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
	}
	else {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}
      }
      else if(score_left_up > score_up) {
	/* left_up + left */
	if(last_left_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a' ;
	}
	else if(last_up == 1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else if(last_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}	
      }
      else if(score_up > score_left_up) {
	/* up + left */
	if(last_up == 1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else if(last_left == 1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
	}
	else if(last_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else if(last_left == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}
	else {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
      }
      else {
	/* all scores equal */
	if(last_left_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}
	else if(last_up == 1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else if(last_left == 1) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
	}
	else if(last_up == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'u';
	}
	else if(last_left == 0) {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'l';
	}
	else {
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->score = score_left_up;
	  (align_mtx + get_mtx_index(j,i,seq_len_target+1))->last = 'a';
	}
      }
    }
  }
#ifdef DEBUG_ALIGN
  dump_align_matrix(seq_len_template+1, seq_len_target+1, align_mtx);
#endif
  
  target_pos = seq_len_target;
  template_pos = seq_len_template;
  
  alignment = (struct alignment_s*)(malloc_or_die(alignment_length * sizeof(struct alignment_s)));
  cur_alignment_pos = 0;
  while((align_mtx + get_mtx_index(template_pos, target_pos, seq_len_target+1))->last != 'e') {
    /* get letter combination */
    score = (align_mtx + get_mtx_index(template_pos, target_pos, seq_len_target+1))->score;
    last = (align_mtx + get_mtx_index(template_pos, target_pos, seq_len_target+1))->last;
    if(last == 'a') {
      (alignment + cur_alignment_pos)->target_pos = target_pos;
      (alignment + cur_alignment_pos)->template_pos = template_pos;
      if(seq_format == FASTA || seq_format == STANDARD) {
	strncpy((alignment + cur_alignment_pos)->target_letter, seq_info_target.seqs->seq_1 + target_pos-1, 5);
	strncpy((alignment + cur_alignment_pos)->template_letter, seq_info_template.seqs->seq_1 + template_pos-1, 5);
      }
      else if(seq_format == PROFILE || seq_format == MSA_STANDARD) {
	strncpy((alignment + cur_alignment_pos)->target_letter,
		(msa_seq_infop_target->msa_seq_1 + ((target_pos-1) * (hmm.a_size+1)))->query_letter, 5);
	strncpy((alignment + cur_alignment_pos)->template_letter,
		(msa_seq_infop_template->msa_seq_1 + ((template_pos-1) * (hmm.a_size+1)))->query_letter, 5);
      }
      target_pos--;
      template_pos--;
    }
    else if(last == 'u') {
      (alignment + cur_alignment_pos)->target_pos = -1;
      (alignment + cur_alignment_pos)->template_pos = template_pos;
      if(seq_format == FASTA || seq_format == STANDARD) {
	strncpy((alignment + cur_alignment_pos)->template_letter,
		(msa_seq_infop_target->msa_seq_1 + ((target_pos-1) * (hmm.a_size+1)))->query_letter, 5);
	(alignment + cur_alignment_pos)->target_letter[0] = '-';
	(alignment + cur_alignment_pos)->target_letter[1] = '\0';
      }
      else if(seq_format == PROFILE || seq_format == MSA_STANDARD) {
	strncpy((alignment + cur_alignment_pos)->template_letter,
		(msa_seq_infop_template->msa_seq_1 + ((template_pos-1) * (hmm.a_size+1)))->query_letter, 5);
	(alignment + cur_alignment_pos)->target_letter[0] = '-';
	(alignment + cur_alignment_pos)->target_letter[1] = '\0';
      }
      template_pos--;
    }
    else if(last == 'l') {
      (alignment + cur_alignment_pos)->target_pos = target_pos;
      (alignment + cur_alignment_pos)->template_pos = -1;
      if(seq_format == FASTA || seq_format == STANDARD) {
	strncpy((alignment + cur_alignment_pos)->target_letter, seq_info_target.seqs->seq_1 + target_pos-1, 5);
	(alignment + cur_alignment_pos)->template_letter[0] = '-';
	(alignment + cur_alignment_pos)->template_letter[1] = '\0';
      }
      else if(seq_format == PROFILE || seq_format == MSA_STANDARD) {
	strncpy((alignment + cur_alignment_pos)->target_letter,
		(msa_seq_infop_target->msa_seq_1 + ((target_pos-1) * (hmm.a_size+1)))->query_letter, 5);
	(alignment + cur_alignment_pos)->template_letter[0] = '-';
	(alignment + cur_alignment_pos)->template_letter[1] = '\0';
      }
      target_pos--;
    }
    else {

    }
    cur_alignment_pos++;
    
  }
#ifdef DEBUG_ALIGN
  for(i = cur_alignment_pos - 1; i >= 0; i--) {
    printf("%d\t%d\n", (alignment + i)->target_pos, (alignment + i)->template_pos);
  }
  for(i = cur_alignment_pos - 1; i >= 0; i--) {
    printf("%s\t%s\n", (alignment + i)->target_letter, (alignment + i)->template_letter);
  }
#endif
  
  
  fprintf(outfile, "target\ttemplate\n");
  for(i = cur_alignment_pos - 1; i >= 0; i--) {
    fprintf(outfile, "%s\t%s\n", (alignment + i)->target_letter, (alignment + i)->template_letter);
  }
  fclose(outfile);
  

  free(alignment);
  free(align_mtx);
}