+++ /dev/null
-/**************** Statistical Significance Parameter Subroutine ****************
-
- $Name: fa_34_26_5 $ - $Id: karlin.c,v 1.18 2006/06/01 16:05:30 wrp Exp $
-
- Version 1.0 February 2, 1990
- Version 2.0 March 18, 1993
-
- Program by: Stephen Altschul
-
- Address: National Center for Biotechnology Information
- National Library of Medicine
- National Institutes of Health
- Bethesda, MD 20894
-
- Internet: altschul@ncbi.nlm.nih.gov
-
- See: Karlin, S. & Altschul, S.F. "Methods for Assessing the Statistical
- Significance of Molecular Sequence Features by Using General Scoring
- Schemes," Proc. Natl. Acad. Sci. USA 87 (1990), 2264-2268.
-
- Computes the parameters lambda and K for use in calculating the
- statistical significance of high-scoring segments or subalignments.
-
- The scoring scheme must be integer valued. A positive score must be
- possible, but the expected (mean) score must be negative.
-
- A program that calls this routine must provide the value of the lowest
- possible score, the value of the greatest possible score, and a pointer
- to an array of probabilities for the occurence of all scores between
- these two extreme scores. For example, if score -2 occurs with
- probability 0.7, score 0 occurs with probability 0.1, and score 3
- occurs with probability 0.2, then the subroutine must be called with
- low = -2, high = 3, and pr pointing to the array of values
- { 0.7, 0.0, 0.1, 0.0, 0.0, 0.2 }. The calling program must also provide
- pointers to lambda and K; the subroutine will then calculate the values
- of these two parameters. In this example, lambda=0.330 and K=0.154.
-
- The parameters lambda and K can be used as follows. Suppose we are
- given a length N random sequence of independent letters. Associated
- with each letter is a score, and the probabilities of the letters
- determine the probability for each score. Let S be the aggregate score
- of the highest scoring contiguous segment of this sequence. Then if N
- is sufficiently large (greater than 100), the following bound on the
- probability that S is greater than or equal to x applies:
-
- P( S >= x ) <= 1 - exp [ - KN exp ( - lambda * x ) ].
-
- In other words, the p-value for this segment can be written as
- 1-exp[-KN*exp(-lambda*S)].
-
- This formula can be applied to pairwise sequence comparison by assigning
- scores to pairs of letters (e.g. amino acids), and by replacing N in the
- formula with N*M, where N and M are the lengths of the two sequences
- being compared.
-
- In addition, letting y = KN*exp(-lambda*S), the p-value for finding m
- distinct segments all with score >= S is given by:
-
- 2 m-1 -y
- 1 - [ 1 + y + y /2! + ... + y /(m-1)! ] e
-
- Notice that for m=1 this formula reduces to 1-exp(-y), which is the same
- as the previous formula.
-
-*******************************************************************************/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-#define MAXIT 25 /* Maximum number of iterations used in calculating lambda */
-#define NMAP_X 23
-#define NMAP 33
-
-#define TINY 1e-6
-
-/* first build a residue map to automatically put residues in score bins */
-
-#include "defs.h"
-#include "param.h"
-
-/* initialize the Karlin frequency, probability arrays using
- a specific query sequence */
-
-int karlin(int , int, double *, double *, double *);
-static int karlin_k(int , int , double *, double *, double *, double *);
-
-void init_karlin(const unsigned char *aa0, int n0, struct pstruct *ppst,
- double *aa0_f, double **kp)
-{
- int kar_nsq, kar_range, kar_min, kar_max;
-
- const unsigned char *aa0p;
- int i;
- int r_cnt[NMAP+1];
- double fn0, *kar_p;
-
- kar_range = ppst->pam_h - ppst->pam_l + 1;
- if (*kp == NULL) {
- if ((kar_p=(double *)calloc(kar_range+1,sizeof(double)))==NULL) {
- fprintf(stderr," cannot allocate kar_p array: %d\n",kar_range+1);
- exit(1);
- }
- *kp = kar_p;
- }
- kar_nsq = ppst->nsq; /* alphabet size */
- kar_min = ppst->pam_l; /* low pam value */
- kar_max = ppst->pam_h; /* high pam value */
-
- /* must have at least 1 residue of each type */
- r_cnt[NMAP]=0;
- for (i=1; i<=kar_nsq; i++) r_cnt[i]=1;
-
- fn0 = 100.0/(double)(n0+kar_nsq); /* weight of each residue */
-
- aa0p = aa0;
- /* increment residue count for each residue in query sequence */
- while (*aa0p) r_cnt[ppst->hsqx[*aa0p++]]++;
-
- /* map all unmapped residues to 'X' */
- r_cnt[NMAP_X] += r_cnt[NMAP];
-
- for (i=1; i<=kar_nsq; i++) aa0_f[i] = fn0*(double)r_cnt[i];
-}
-
-double nt_f[] = {0.0, 0.25, 0.25, 0.25, 0.25 };
-
-/* Robinson and Robinson frequencies */
-double aa_f[] = {
-/* NULL */ 0.00,
-/* A */ 0.0780474700897585,
-/* R */ 0.0512953149316987,
-/* N */ 0.0448725775979007,
-/* D */ 0.0536397361638076,
-/* C */ 0.0192460110427568,
-/* Q */ 0.0426436013507063,
-/* E */ 0.0629485981204668,
-/* G */ 0.0737715654561964,
-/* H */ 0.0219922696262025,
-/* I */ 0.0514196403000682,
-/* L */ 0.090191394464413,
-/* K */ 0.0574383201866657,
-/* M */ 0.0224251883196316,
-/* F */ 0.0385564048655621,
-/* P */ 0.0520279465667327,
-/* S */ 0.0711984743501224,
-/* T */ 0.0584129422708473,
-/* W */ 0.013298374223799,
-/* Y */ 0.0321647488738564,
-/* V */ 0.0644094211988074};
-
-/* initialize the Karlin frequency, probability arrays using
- an "average" composition (average length if n0 <=0) */
-
-void
-init_karlin_a(struct pstruct *ppst, double *aa0_f, double **kp)
-{
- int kar_nsq, kar_range;
-
- int i;
- double fn0, *kar_p;
-
- kar_range = ppst->pam_h - ppst->pam_l + 1;
- if (*kp == NULL) {
- if ((kar_p=(double *)calloc(kar_range+1,sizeof(double)))==NULL) {
- fprintf(stderr," cannot allocate kar_p array: %d\n",kar_range+1);
- exit(1);
- }
- *kp = kar_p;
- }
-
- if (ppst->nt_align) {
- kar_nsq = 4;
- for (i=1; i<=kar_nsq; i++) aa0_f[i] = nt_f[i];
- }
- else if (ppst->dnaseq==SEQT_PROT || ppst->dnaseq == SEQT_UNK) {
- kar_nsq = 20;
- for (i=1; i<=kar_nsq; i++) aa0_f[i] = aa_f[i];
- }
- else {
- kar_nsq = ppst->nsq;
- fn0 = 1.0/(double)(kar_nsq-1);
- for (i=1; i< kar_nsq; i++) aa0_f[i] = fn0;
- aa0_f[kar_nsq]=0.0;
- }
-
-}
-
-/* calculate set up karlin() to calculate Lambda, K, by calculating
- aa1 frequencies */
-int
-do_karlin(const unsigned char *aa1, int n1,
- int **pam2, struct pstruct *ppst,
- double *aa0_f, double *kar_p, double *lambda, double *H)
-{
- register unsigned const char *aap;
- int kar_range, kar_min, kar_max, kar_nsq;
- int r_cnt[NMAP+1];
- double aa1_f[NMAP];
- double fn1, kar_tot;
- int i, j;
-
- kar_nsq = ppst->nsq;
- kar_min = ppst->pam_l;
- kar_max = ppst->pam_h;
- kar_range = kar_max - kar_min + 1;
-
- r_cnt[NMAP]=0;
- for (i=1; i<=kar_nsq; i++) r_cnt[i]=1;
-
- /* residue counts */
-
- aap=aa1;
- while (*aap) r_cnt[ppst->hsqx[*aap++]]++;
-
- r_cnt[NMAP_X] += r_cnt[NMAP];
-
- /* residue frequencies */
- fn1 = 100.0/(double)(n1+kar_nsq);
- for (i=1; i<=kar_nsq; i++) aa1_f[i]= fn1*(double)r_cnt[i];
-
- for (i=0; i<=kar_range; i++) kar_p[i] = 0.0;
-
- for (i=1; i<=kar_nsq; i++) {
- for (j=1; j<=kar_nsq; j++)
- kar_p[pam2[i][j]-kar_min] += aa0_f[i]*aa1_f[j];
- }
-
- kar_tot = 0.0;
- for (i=0; i<=kar_range; i++) kar_tot += kar_p[i];
- if (kar_tot <= 0.00001) return 0;
-
- for (i=0; i<=kar_range; i++) kar_p[i] /= kar_tot;
-
- return karlin(kar_min, kar_max, kar_p, lambda, H);
-}
-
-int
-do_karlin_a(int **pam2, struct pstruct *ppst,
- double *aa0_f, double *kar_p, double *lambda, double *K, double *H)
-{
- double *aa1fp;
- int kar_range, kar_min, kar_max, kar_nsq;
- double aa1_f[NMAP];
- double fn1, kar_tot;
- int i, j;
-
- kar_min = ppst->pam_l;
- kar_max = ppst->pam_h;
- kar_range = kar_max - kar_min + 1;
-
- kar_tot = 0.0;
- if (ppst->nt_align ) {
- kar_nsq = 4;
- aa1fp = nt_f;
- for (i=1; i<=kar_nsq; i++) {kar_tot += aa1fp[i];}
- for (i=1; i<=kar_nsq; i++) {aa1_f[i]= aa1fp[i]/kar_tot;}
- }
- else if (!ppst->nt_align) {
- kar_nsq = 20;
- aa1fp = aa_f;
- for (i=1; i<=kar_nsq; i++) {kar_tot += aa1fp[i];}
- for (i=1; i<=kar_nsq; i++) {aa1_f[i]= aa1fp[i]/kar_tot;}
- }
- else {
- kar_nsq = ppst->nsq;
- fn1 = 1.0/(double)(kar_nsq-1);
- for (i=1; i< kar_nsq; i++) aa1_f[i] = fn1;
- aa1_f[kar_nsq]=0.0;
- }
-
- for (i=0; i<=kar_range; i++) kar_p[i] = 0.0;
-
- for (i=1; i<=kar_nsq; i++) {
- for (j=1; j<kar_nsq; j++)
- kar_p[pam2[i][j]-kar_min] += aa0_f[i]*aa1_f[j];
- }
-
- kar_tot = 0.0;
- for (i=0; i<=kar_range; i++) kar_tot += kar_p[i];
- if (kar_tot <= 0.00001) return 0;
-
- for (i=0; i<=kar_range; i++) kar_p[i] /= kar_tot;
-
- return karlin_k(kar_min, kar_max, kar_p, lambda, K, H);
-}
-
-/* take a array of letters and pam information and get *lambda, *H */
-int
-karlin(int low, /* Lowest score (must be negative) */
- int high, /* Highest score (must be positive) */
- double *pr, /* Probabilities for various scores */
- double *lambda_p, /* Pointer to parameter lambda */
- double *H_p) /* Pointer to parameter H */
-{
- int i,range, nit;
- double up,new,sum,av,beta,ftemp;
- double lambda;
- double *p,*ptr1;
-
- /* Calculate the parameter lambda */
-
- p = pr;
- range = high-low;
-
- /* check for E() < 0.0 */
- sum = 0;
- ptr1 = pr;
- for (i=low; i <= high ; i++) sum += i* (*ptr1++);
- if (sum >= 0.0) {
-#ifdef DEBUG
- fprintf(stderr," (karlin lambda) non-negative expected score: %.4lg\n",
- sum);
-#endif
- return 0;
- }
-
- /* up is upper bound on lambda */
- up=0.5;
- do {
- up *= 2.0;
- ptr1=p;
-
- beta=exp(up);
-
- ftemp=exp(up*(low-1));
- sum = 0.0;
- for (i=0; i<=range; ++i) sum+= *ptr1++ * (ftemp*=beta);
- }
- while (sum<1.0);
-
- /* avoid overflow from very large lambda*S */
-/*
- do {
- up /= 2.0;
- ptr1=p;
- beta=exp(up);
-
- ftemp=exp(up*(low-1));
- sum = 0.0;
- for (i=0; i<=range; ++i) sum+= *ptr1++ * (ftemp*=beta);
- } while (sum > 2.0);
-
- up *= 2.0;
-*/ /* we moved past, now back up */
-
- /* for (lambda=j=0;j<25;++j) { */
- lambda = 0.0;
- nit = 0;
- while ( nit++ < MAXIT ) {
- new = (lambda+up)/2.0;
- beta = exp(new);
- ftemp = exp(new*(low-1));
- ptr1=p;
- sum = 0.0;
- /* multiply by exp(new) for each score */
- for (i=0;i<=range;++i) sum+= *ptr1++ * (ftemp*=beta);
-
- if (sum > 1.0 + TINY) up=new;
- else {
- if ( fabs(lambda - new) < TINY ) goto done;
- lambda = new;
- }
- }
-
- if (lambda <= 1e-10) {
- lambda = -1.0;
- return 0;
- }
-
- done:
- *lambda_p = lambda;
-
- /* Calculate the parameter K */
-
- ptr1=p;
- ftemp=exp(lambda*(low-1));
- for (av=0.0, i=low; i<=high; ++i)
- av+= *ptr1++ *i*(ftemp*=beta);
- *H_p= lambda*av;
-
- return 1; /* Parameters calculated successfully */
-}
-
-static int a_gcd (int, int);
-
-/* take a array of letters and pam information and get *lambda, *K, *H */
-static int
-karlin_k(int low, /* Lowest score (must be negative) */
- int high, /* Highest score (must be positive) */
- double *pr, /* Probabilities for various scores */
- double *lambda_p, /* Pointer to parameter lambda */
- double *K_p,
- double *H_p) /* Pointer to parameter H */
-{
- int i,j,range,lo,hi,first,last, nit;
- double up,new,sum,Sum,av,beta,oldsum,ratio,ftemp;
- double lambda;
- double *p,*P,*ptrP,*ptr1,*ptr2;
-
- /* Calculate the parameter lambda */
-
- p = pr;
- range = high-low;
-
- /* check for E() < 0.0 */
- sum = 0;
- ptr1 = pr;
- for (i=low; i <= high ; i++) sum += i* (*ptr1++);
- if (sum >= 0.0) {
-#ifdef DEBUG
- fprintf(stderr," (karlin lambda) non-negative expected score: %.4lg\n",
- sum);
-#endif
- return 0;
- }
-
- /* up is upper bound on lambda */
- up=0.5;
- do {
- up *= 2.0;
- ptr1=p;
-
- beta=exp(up);
-
- ftemp=exp(up*(low-1));
- sum = 0.0;
- for (i=0; i<=range; ++i) sum+= *ptr1++ * (ftemp*=beta);
- }
- while (sum<1.0);
-
- /* avoid overflow from very large lambda*S */
- /*
- do {
- up /= 2.0;
- ptr1=p;
- beta=exp(up);
-
- ftemp=exp(up*(low-1));
- sum = 0.0;
- for (i=0; i<=range; ++i) sum+= *ptr1++ * (ftemp*=beta);
- } while (sum > 2.0);
-
- up *= 2.0;
- */
- /* we moved past, now back up */
-
- /* for (lambda=j=0;j<25;++j) { */
- lambda = 0.0;
- nit = 0;
- while ( nit++ < MAXIT ) {
- new = (lambda+up)/2.0;
- beta = exp(new);
- ftemp = exp(new*(low-1));
- ptr1=p;
- sum = 0.0;
- /* multiply by exp(new) for each score */
- for (i=0;i<=range;++i) sum+= *ptr1++ * (ftemp*=beta);
-
- if (sum > 1.0 + TINY) up=new;
- else {
- if ( fabs(lambda - new) < TINY ) goto done;
- lambda = new;
- }
- }
-
- if (lambda <= 1e-10) {
- lambda = -1.0;
- return 0;
- }
-
- done:
- *lambda_p = lambda;
-
- /* Calculate the parameter H */
-
- ptr1=p;
- ftemp=exp(lambda*(low-1));
- for (av=0.0, i=low; i<=high; ++i) av+= *ptr1++ *i*(ftemp*=beta);
- *H_p= lambda*av;
-
- /* Calculate the pamameter K */
- Sum=lo=hi=0;
- P= (double *) calloc(MAXIT*range+1,sizeof(double));
- for (*P=sum=oldsum=j=1;j<=MAXIT && sum>0.001;Sum+=sum/=j++) {
- first=last=range;
- for (ptrP=P+(hi+=high)-(lo+=low); ptrP>=P; *ptrP-- =sum) {
- ptr1=ptrP-first;
- ptr2=p+first;
- for (sum=0,i=first; i<=last; ++i) sum += *ptr1-- * *ptr2++;
- if (first) --first;
- if (ptrP-P<=range) --last;
- }
- ftemp=exp(lambda*(lo-1));
- for (sum=0,i=lo;i;++i) sum+= *++ptrP * (ftemp*=beta);
- for (;i<=hi;++i) sum+= *++ptrP;
- ratio=sum/oldsum;
- oldsum=sum;
- }
- for (;j<=200;Sum+=oldsum/j++) oldsum*=ratio;
- for (i=low;!p[i-low];++i);
- for (j= -i;i<high && j>1;) if (p[++i-low]) j=a_gcd(j,i);
- *K_p = (j*exp(-2*Sum))/(av*(1.0-exp(- lambda*j)));
- free(P);
-
- return 1; /* Parameters calculated successfully */
-}
-
-int
-a_gcd(int a, int b)
-{
- int c;
-
- if (b<0) b= -b;
- if (b>a) { c=a; a=b; b=c; }
- for (;b;b=c) { c=a%b; a=b; }
- return a;
-}
-
git://source.jalview.org / jabaws.git / blobdiff