+++ /dev/null
-/*
- * dist.c
- *
- * $Id$
- *
- *
- *****************************************************************************
- *
- * Copyright (c) 2004, Luke Sheneman
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * + Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * + Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in
- * the documentation and/or other materials provided with the
- * distribution.
- * + The names of its contributors may not be used to endorse or promote
- * products derived from this software without specific prior
- * written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- *****************************************************************************
- *
- * Compute a distance matrix given a set of sequences
- *
- *****************************************************************************
- *
- * AUTHOR:
- *
- * Luke Sheneman
- * sheneman@cs.uidaho.edu
- *
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <string.h>
-#include <ctype.h>
-
-#include "common.h"
-#include "dayhoff.h"
-#include "fasta.h"
-#include "distclearcut.h"
-
-
-
-
-/*
- * NJ_build_distance_matrix() -
- *
- * Given a filename for an alignment, read the alignment
- * into memory and then compute the distance matrix
- * using the appropriate correction model
- */
-DMAT *
-NJ_build_distance_matrix(NJ_ARGS *nj_args) {
-
- DMAT *dmat;
- NJ_alignment *alignment;
-
- /* Read an alignment in FASTA format */
- alignment =
- NJ_read_fasta(nj_args);
-
- if(!alignment) {
- return(NULL);
- }
-
- /*
- * Given a global multiple sequence alignment (MSA) and
- * a specified distance correction model, compute a
- * corrected distance matrix
- *
- * From proteins, we may want to allow users to specify
- * a substitution matrix (feature)
- */
-
- dmat =
- NJ_compute_dmat(nj_args,
- alignment);
-
- // NJ_print_taxanames(dmat);
-
- if(!dmat) {
- fprintf(stderr, "Clearcut: Error computing distance matrix\n");
- }
-
- /* now free the memory associated with the alignment */
- NJ_free_alignment(alignment);
-
- return(dmat);
-}
-
-
-
-
-
-/*
- * NJ_compute_dmat() -
- *
- * Given an alignment and a correction model, compute the
- * distance matrix and return it
- *
- */
-DMAT *
-NJ_compute_dmat(NJ_ARGS *nj_args,
- NJ_alignment *alignment) {
-
- DMAT *dmat;
- long int i;
-
-
- /* allocate distance matrix here */
- dmat = (DMAT *)calloc(1, sizeof(DMAT));
- if(!dmat) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_compute_dmat()\n");
- return(NULL);
- }
-
- dmat->ntaxa = alignment->nseq;
- dmat->size = alignment->nseq;
-
- /* allocate memory to hold the taxa names */
- dmat->taxaname = (char **)calloc(alignment->nseq, sizeof(char *));
- if(!dmat->taxaname) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_compute_dmat()\n");
- return(NULL);
- }
-
- /* copy sequence titles */
- for(i=0;i<alignment->nseq;i++) {
- dmat->taxaname[i] = (char *)calloc(strlen(alignment->titles[i])+1, sizeof(char));
- if(!dmat->taxaname[i]) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_compute_dmat()\n");
- return(NULL);
- }
-
- strncpy(dmat->taxaname[i], alignment->titles[i], strlen(alignment->titles[i]));
- }
-
- /* allocate val matrix in dmat */
- dmat->val = (float *)calloc(dmat->ntaxa*dmat->ntaxa, sizeof(float));
-
- if(!dmat->val) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_compute_dmat()\n");
- return(NULL);
- }
-
- /* now lets allocate space for the r and r2 columns */
- dmat->r = (float *)calloc(dmat->ntaxa, sizeof(float));
- dmat->r2 = (float *)calloc(dmat->ntaxa, sizeof(float));
-
- /* track some memory addresses */
- dmat->rhandle = dmat->r;
- dmat->r2handle = dmat->r2;
- dmat->valhandle = dmat->val;
-
- /* apply model correction to matrix */
- switch(nj_args->correction_model) {
-
- case NJ_MODEL_JUKES:
-
- if(nj_args->dna_flag) {
- NJ_DNA_jc_correction(dmat, alignment);
- } else if(nj_args->protein_flag) {
- NJ_PROTEIN_jc_correction(dmat, alignment);
- } else {
- fprintf(stderr, "Clearcut: Need to know sequence type for Jukes-Cantor model correction.\n");
- return(NULL);
- }
-
- break;
-
- case NJ_MODEL_KIMURA:
-
- if(nj_args->dna_flag) {
- NJ_DNA_k2p_correction(dmat, alignment);
- } else if(nj_args->protein_flag) {
- NJ_PROTEIN_kimura_correction(dmat, alignment);
- } else {
- fprintf(stderr, "Clearcut: Need to know sequence type for Kimura model correction.\n");
- return(NULL);
- }
-
- break;
-
- case NJ_MODEL_NONE:
-
- NJ_no_correction(dmat, alignment);
-
- break;
-
- default:
- fprintf(stderr, "Clearcut: Invalid distance correction model.\n");
- return(NULL);
- }
-
- return(dmat);
-}
-
-
-
-
-
-/*
- * NJ_no_correction() -
- *
- * Compute the distance matrix without correction
- * (straight percent ID)
- *
- * Resolve ambiguities in sequence data by skipping
- * those nucleotides/residues
- *
- */
-void
-NJ_no_correction(DMAT *dmat,
- NJ_alignment *alignment) {
-
- long int i, j;
- float pdiff;
-
- /* compute pairwise percent identity */
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
- pdiff = 1.0 - NJ_pw_percentid(alignment, i, j);
- dmat->val[NJ_MAP(i, j, dmat->size)] = pdiff;
- }
- }
-
- return;
-}
-
-
-
-
-/*
- * NJ_DNA_jc_correction() -
- *
- * Compute the distance matrix with jukes-cantor correction
- * and assign high distance if sequence divergence exceeds
- * 0.75
- *
- * Jukes, T.H. (1969), Evolution of protein molecules. In H.N. Munro (Ed.),
- * Mammalian Protein Metabolism, Volume III, Chapter 24, pp. 21-132.
- * New York: Academic Press
- *
- */
-void
-NJ_DNA_jc_correction(DMAT *dmat,
- NJ_alignment *alignment) {
-
- long int i, j;
- long int k;
- float d, cutoff, dist;
- long int residues;
-
- cutoff = 0.75;
-
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
-
- k = NJ_pw_differences(alignment, i, j, &residues);
- d = 1.0 - NJ_pw_percentid(alignment, i, j);
-
- if(d > cutoff) {
- dist = NJ_BIGDIST;
- } else {
- dist = (-0.75) * log(1.0 - (4.0/3.0)*d);
- }
-
- if(fabs(dist) < FLT_EPSILON) {
- dmat->val[NJ_MAP(i, j, dmat->size)] = 0.0;
- } else {
- dmat->val[NJ_MAP(i, j, dmat->size)] = dist;
- }
- }
- }
-
-
-
- return;
-}
-
-
-
-
-
-
-/*
- * NJ_PROTEIN_jc_correction() -
- *
- * This function performs modified jukes/cantor correction on
- * a protein alignment
- *
- * Jukes, T.H. (1969), Evolution of protein molecules. In H.N. Munro (Ed.),
- * Mammalian Protein Metabolism, Volume III, Chapter 24, pp. 21-132.
- * New York: Academic Press
- *
- */
-void
-NJ_PROTEIN_jc_correction(DMAT *dmat,
- NJ_alignment *alignment) {
-
- long int i, j;
- long int residues;
- long int diff;
- float dist, x;
-
-
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
-
- diff = NJ_pw_differences(alignment, i, j, &residues);
-
- if(!diff || !residues) {
- dist = 0.0;
- } else {
-
- dist = (float)diff/(float)residues;
- x = ((20.0/19.0)*dist);
-
- if(NJ_FLT_GT(x, 1.0)) {
- dist = NJ_BIGDIST;
- } else {
- dist = -(19.0/20.0) * log(1.0 - x);
- }
- }
-
- dmat->val[NJ_MAP(i, j, dmat->size)] = dist;
- }
- }
-
- return;
-}
-
-
-
-
-
-
-/*
- * NJ_DNA_k2p_correction() -
- *
- * Correct a distance matrix using k2p correction using
- * cutoffs to avoid problems with logarithms.
- *
- * dist = -0.5ln(1-2P-Q) - 0.25ln(1-2Q)
- *
- * But due to the logarithms, this is only valid when
- *
- * (2P+Q <= 1) &&
- * (2Q <= 1)
- *
- * So assign arbitary distances when these constraints are
- * not strictly followed.
- *
- * Kimura, M. (1980), A simple method for estimating evolutionary
- * rates of base substitutions through comparative studies of
- * nucleotide sequences. J. Mol. Evol., 16, 111-120
- *
- */
-void
-NJ_DNA_k2p_correction(DMAT *dmat,
- NJ_alignment *alignment) {
-
- long int i, j;
- float P; /* proportion of transitions */
- float Q; /* proportion of transversions */
- long int nucleotides;
- long int transitions, transversions;
- float dist;
- float log_x = 0.0; /* the params for the first log */
- float log_y = 0.0; /* the params for the second log */
-
- int blowup; /* a flag to specify if we have a log blowup */
-
-
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
-
- blowup = 0;
-
- /* count the number of transitions and transversions */
- NJ_DNA_count_tt(alignment, i, j, &transitions, &transversions, &nucleotides);
-
- if(!nucleotides) { /* sequences have no non-ambiguous overlap in alignment */
- P = 0.0;
- Q = 0.0;
- } else {
- P = (float)transitions / (float)nucleotides;
- Q = (float)transversions / (float)nucleotides;
- }
-
- /* the first log blows up if 2*P+Q = 1.0 */
- if(NJ_FLT_EQ((2.0 * P + Q), 1.0)) {
- blowup = 1;
- } else {
- if( NJ_FLT_LT(1.0 - 2.0*P - Q, 0.0) ) {
- blowup = 1;
- } else {
- log_x = log(1.0 - 2.0*P - Q);
- }
- }
-
- /* the second log blows up if 2*Q >= 1.0 */
- if( NJ_FLT_EQ((2.0 * Q), 1.0) ||
- NJ_FLT_GT((2.0 * Q), 1.0) ) {
- blowup = 1;
- } else {
- log_y = log(1.0 - 2.0*Q);
- }
-
- /* if our logarithms blow up, we just set the distance to the max */
- if(blowup) {
- dist = NJ_BIGDIST;
- } else {
- dist = (-0.5)*log_x - 0.25*log_y;
- }
-
- if(fabs(dist) < FLT_EPSILON) {
- dmat->val[NJ_MAP(i, j, dmat->size)] = 0.0;
- } else {
- dmat->val[NJ_MAP(i, j, dmat->size)] = dist;
- }
- }
- }
-
- return;
-}
-
-
-
-
-/*
- * NJ_PROTEIN_kimura_correction() -
- *
- * Perform Kimura correction for distances derived from protein
- * alignments.
- *
- * Kimura, M. (1983), The Neutral Theory of Molecular Evolution.
- * p. 75., Cambridge University Press, Cambridge, England
- *
- */
-void
-NJ_PROTEIN_kimura_correction(DMAT *dmat,
- NJ_alignment *alignment) {
-
- long int i, j;
- long int residues;
- long int diff;
- float dist;
-
-
- printf("NJ_PROTEIN_kimura_correction()\n");
-
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
- diff = NJ_pw_differences(alignment, i, j, &residues);
-
- if(!diff || !residues) {
- dist = 0.0;
- } else {
- dist = (float)diff/(float)residues;
- }
-
- if(NJ_FLT_LT(dist, 0.75)) {
- if(NJ_FLT_GT(dist, 0.0) ) {
- dist = -log(1.0 - dist - (dist * dist/5.0) );
- }
- } else {
- if(NJ_FLT_GT(dist, 0.93) ) {
- dist = 10.0;
- } else {
- dist = (float)NJ_dayhoff[ (int)((dist*1000.0)-750.0) ] / 100.0 ;
- }
- }
-
- dmat->val[NJ_MAP(i, j, dmat->size)] = dist;
- }
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_DNA_count_tt() -
- *
- * Count the number of transitions and transversions
- * between two aligned DNA sequences
- *
- * This routine automatically skips ambiguities when
- * counting transitions and transversions.
- *
- */
-void
-NJ_DNA_count_tt(NJ_alignment *alignment,
- long int x,
- long int y,
- long int *transitions,
- long int *transversions,
- long int *residues) {
-
- long int tmp_transitions = 0;
- long int tmp_transversions = 0;
- long int tmp_residues = 0;
- char a, b;
- long int i;
-
- for(i=0;i<alignment->length;i++) {
-
- a = toupper(alignment->data[x*alignment->length+i]);
- b = toupper(alignment->data[y*alignment->length+i]);
-
- if( (a == 'A' && b == 'T') ||
- (a == 'T' && b == 'A') ||
- (a == 'A' && b == 'C') ||
- (a == 'C' && b == 'A') ||
- (a == 'T' && b == 'G') ||
- (a == 'G' && b == 'T') ||
- (a == 'C' && b == 'G') ||
- (a == 'G' && b == 'C') ) {
- tmp_transversions++;
- }
-
- if( (a == 'C' && b == 'T') ||
- (a == 'T' && b == 'C') ||
- (a == 'G' && b == 'A') ||
- (a == 'A' && b == 'G') ) {
- tmp_transitions++;
- }
-
- /* count the number of residues */
- if(a != NJ_AMBIGUITY_CHAR &&
- b != NJ_AMBIGUITY_CHAR ) {
- tmp_residues++;
- }
-
- }
-
- *transitions = tmp_transitions;
- *transversions = tmp_transversions;
-
- if(residues) {
- *residues = tmp_residues;
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_pw_percentid() -
- *
- * Given an alignment and a specification
- * for two rows, compute the pairwise
- * percent identity between the two
- *
- */
-float
-NJ_pw_percentid(NJ_alignment *alignment,
- long int x,
- long int y) {
-
- float pid;
- long int i;
- long int residues;
- long int same;
- char c1, c2;
-
- residues = 0;
- same = 0;
- for(i=0;i<alignment->length;i++) {
-
- c1 = alignment->data[x*alignment->length+i];
- c2 = alignment->data[y*alignment->length+i];
-
- if( c1 != NJ_AMBIGUITY_CHAR ||
- c2 != NJ_AMBIGUITY_CHAR ) {
-
- residues++;
-
- if(c1 == c2) {
- same++;
- }
- }
-
- }
-
- pid = (float)same/(float)residues;
-
- return(pid);
-}
-
-
-
-/*
- * NJ_pw_differences() -
- *
- * Given an alignment and a specification
- * for two rows in the alignment, compute the
- * number of differences between the two sequences
- *
- * With respect to ambiguity codes, we will want to
- * disregard those sites entirely in our count.
- *
- */
-long int
-NJ_pw_differences(NJ_alignment *alignment,
- long int x,
- long int y,
- long int *residues) {
-
- long int i;
- long int diff;
- char c1, c2;
- long int tmp_residues;
-
- diff = 0;
- tmp_residues = 0;
- for(i=0;i<alignment->length;i++) {
-
- c1 = alignment->data[x*alignment->length+i];
- c2 = alignment->data[y*alignment->length+i];
-
- if( c1 != NJ_AMBIGUITY_CHAR ||
- c2 != NJ_AMBIGUITY_CHAR ) {
-
- tmp_residues++;
-
- if(c1 != c2) {
- diff++;
- }
- }
-
- }
-
- *residues = tmp_residues;
-
- return(diff);
-}
-
-
-
-
-
-
projects / mothur.git / blobdiff