+++ /dev/null
-/**
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <stdio.h>
-#include <cstdlib>
-#include <iostream>
-#include <algorithm>
-#include <string>
-#include <stdio.h>
-#include "../alignment/Alignment.h"
-#ifdef HAVE_MHASH_H
-#include "mhash.h"
-#endif
-
-#include "Stats.h"
-
-
-using namespace std;
-
-namespace clustalw
-{
-
-
-Stats::Stats()
-{
- enabled = false;
-}
-
-
-Stats::~Stats()
-{
-}
-
-
-
-
-/* adopted from Sean Eddy'ssquid:aligneval.c */
-/* Function: PairwiseIdentity()
- *
- * Purpose: Calculate the pairwise fractional identity between
- * two aligned sequences s1 and s2. This is simply
- * (idents / MIN(len1, len2)).
- *
- * Note how many ways there are to calculate pairwise identity,
- * because of the variety of choices for the denominator:
- * idents/(idents+mismat) has the disadvantage that artifactual
- * gappy alignments would have high "identities".
- * idents/(AVG|MAX)(len1,len2) both have the disadvantage that
- * alignments of fragments to longer sequences would have
- * artifactually low "identities".
- *
- * Case sensitive; also, watch out in nucleic acid alignments;
- * U/T RNA/DNA alignments will be counted as mismatches!
- */
-/* float
-PairwiseIdentity(char @s1, char @s2)
-{
- int idents; /@ total identical positions @/
- int len1, len2; /@ lengths of seqs @/
- int x; /@ position in aligned seqs @/
-
- idents = len1 = len2 = 0;
- for (x = 0; s1[x] != '\0' && s2[x] != '\0'; x++)
- {
- if (!isgap(s1[x])) {
- len1++;
- if (s1[x] == s2[x]) idents++;
- }
- if (!isgap(s2[x])) len2++;
- }
- if (len2 < len1) len1 = len2;
- return (len1 == 0 ? 0.0 : (float) idents / (float) len1);
-}
-*/
-
-
-/* s1/s2 are unit-offset as usual
- *
- */
-float
-Stats::pairwiseIdentity(Alignment *alnObj, int s1, int s2)
-{
- int idents; /* total identical positions */
- int len1, len2; /* real lengths of seqs */
- int x; /* position in aligned seqs */
-
- const vector<int>* seq1 = alnObj->getSequence(s1);
- const vector<int>* seq2 = alnObj->getSequence(s2);
-
- idents = len1 = len2 = 0;
- // cerr << "comparing " << alnObj->getName(s1).c_str() << ":" << alnObj->getName(s2).c_str() << " " << s1 << ":" << s2 << "\n";
-
- // sequence length should be identical, but be paranoid
- for (x = 1; x<=alnObj->getSeqLength(s1) && x<=alnObj->getSeqLength(s2); x++)
- {
- if (! alnObj->isGap(s1, x)) {
- len1++;
- //cerr << " pos " << x << ": " << (*seq1)[x] << ":" << (*seq2)[x] << "\n";
- if ((*seq1)[x] == (*seq2)[x])
- idents++;
- }
- //DEBUG
- //else {
- //cerr << " gap at pos " << x << " (" << s1 << ")\n";
- //}
- if (! alnObj->isGap(s2, x))
- len2++;
- //DEBUG
- //else
- // cerr << " gap at pos " << x << " (" << s2 << ")\n";
- }
- if (len2 < len1)
- len1 = len2;
- return (len1 == 0 ? 0.0 : (float) idents / (float) len1);
-}
-
-
-
-#ifdef HAVE_MHASH_H
-
-string
-Stats::ConcatInputHash(Alignment *alnObj)
-{
- vector<string> rawSeqArray;
- string ret;
- char *hash;
-
- // collect all sequences and sort
- const clustalw::SeqArray* seqArray = alnObj->getSeqArray();
- for(int s = 1; s <= alnObj->getNumSeqs(); s++)
- {
- string seq;
- for(int r = 1; r <= alnObj->getSeqLength(s); r++)
- {
- int val = (*seqArray)[s][r];
- seq.append(1, clustalw::userParameters->getAminoAcidCode(val));
- }
- rawSeqArray.push_back(seq);
- }
- std::sort(rawSeqArray.begin(), rawSeqArray.end());
-
- // concatenate sorted seqs
- string concatSeq;
- std::vector<string>::iterator iter;
- for(iter=rawSeqArray.begin(); iter != rawSeqArray.end(); ++iter)
- concatSeq.append(*iter);
-
- // build hash and return
- hash = Md5Hash(concatSeq.c_str());
- if (hash==NULL)
- {
- ret="HASHING_FAILURE";
- } else {
- for (int i=0; i<strlen(hash); i++)
- ret.append(1, hash[i]);
- free(hash);
- }
- return ret;
-}
-
-
-string
-Stats::Md5ForSeq(Alignment *alnObj, int s)
-{
- string seq;
- const clustalw::SeqArray* seqArray = alnObj->getSeqArray();
- char *hash;
- string ret;
-
- for(int l = 1; l <= alnObj->getSeqLength(s); l++)
- {
- int val = (*seqArray)[s][l];
- // continue if gap
- if((val < 0) || (val > userParameters->getMaxAA()))
- continue;
- seq.append(1, clustalw::userParameters->getAminoAcidCode(val));
- }
- std::transform(seq.begin(), seq.end(), seq.begin(), ::toupper);
-
- hash = Md5Hash(seq.c_str());
- if (hash==NULL)
- {
- ret = "HASHING_FAILURE";
- } else {
- ret = hash;
- }
- return ret;
-}
-
-/* create md5 hash from input string
- * returns NULL on failure
- * user must free returned string
- */
-char *
-Stats::Md5Hash(const char *thread)
-{
- MHASH td;
- char *tmpcstr;
- int i;
- unsigned char *hash;
- char *rethash;
-
- td = mhash_init(MHASH_MD5);
- if (td == MHASH_FAILED)
- return NULL;
-
- if (thread==NULL)
- return NULL;
-
- mhash(td, thread, strlen(thread));
- //mhash_deinit(td, hash);
- hash = (unsigned char*) mhash_end(td);
-
- rethash = (char*) calloc(mhash_get_block_size(MHASH_MD5)*2+1, sizeof(char));
- for (i = 0; i < mhash_get_block_size(MHASH_MD5); i++) {
- sprintf(&rethash[i*2], "%.2x", hash[i]);
- }
-
- return rethash;
-}
-
-#endif
-// HAVE_MHASH_H
-
-
-
-
-void
-Stats::logCmdLine(int argc, char **argv)
-{
- FILE *fp = fopen(logfilename.c_str(), "a");
- if (fp == NULL)
- {
- cerr << "couldn't open file " << logfilename << " for logging of stats\n";
- return;
- }
-
- if (argc > 1)
- {
- for (int i=1; i<argc; i++)
- {
- // remove non-interesting stuff
- if (strstr(argv[i], "-infile=")==NULL &&
- strstr(argv[i], "-outfile=")==NULL &&
- strstr(argv[i], "-stats=")==NULL &&
- strstr(argv[i], "-align")==NULL)
- {
- fprintf(fp, "cmdline non-default arg: %s\n", argv[i]);
- }
- }
- }
-}
-
-
-/* log some statistics for input sequences, i.e. alnObj should hold
- * the unaligned sequences
- *
- */
-void
-Stats::logInputSeqStats(Alignment *alnObj)
-{
- int i;
- std::vector<double> lengths;
- time_t t = time(NULL);
- tm s = *localtime(&t);
- int shortest;
- string hash;
-
- FILE *fp = fopen(logfilename.c_str(), "a");
- if (fp == NULL)
- {
- cerr << "couldn't open file " << logfilename << " for logging of stats\n";
- return;
- }
-
- fprintf(fp, "logging job: %s on %s", userParameters->getSeqName().c_str(), asctime(&s));
- fprintf(fp, "clustal version: %s\n", userParameters->getRevisionLevel().c_str());
-
- fprintf(fp, "seq type: ");
- if (userParameters->getDNAFlag())
- fprintf(fp, "DNA");
- else
- fprintf(fp, "protein");
- fprintf(fp, "\n");
-
- fprintf(fp, "numseqs: %d\n", alnObj->getNumSeqs());
-
- // create a vector of seq lengths for later
- // and get shortest seq at the same time
- shortest=alnObj->getLengthLongestSequence();
- for (i = 1; i <= alnObj->getNumSeqs(); i++) {
- int l = alnObj->getSeqLength(i);
- lengths.push_back(l);
- if (l<shortest)
- shortest=l;
- }
-
- fprintf(fp, "seqlen longest: %d\n", alnObj->getLengthLongestSequence());
- fprintf(fp, "seqlen shortest: %d\n", shortest);
-
- fprintf(fp, "seqlen avg: %.2f\n", utilityObject->average(lengths));
- fprintf(fp, "seqlen std-dev: %.2f\n", utilityObject->stdDev(lengths));
- fprintf(fp, "seqlen median: %.2f\n", utilityObject->median(lengths));
-
-
-#ifdef HAVE_MHASH_H
- //hash = concatInputHash(alnObj);
- //fprintf(fp, "seq hash: %s\n", hash.c_str());
-
- for (int s = 1; s <= alnObj->getNumSeqs(); s++) {
- string md5 = Md5ForSeq(alnObj, s);
- fprintf(fp, "md5 for seq %d: %s\n", s, md5.c_str());
-
- }
-#else
- fprintf(fp, "md5: disabled\n");
-#endif
-
- fclose(fp);
-}
-
-
-
-/* log some statistics for aligned sequences, i.e. alnObj should hold
- * the already aligned sequences
- *
- */
-void
-Stats::logAlignedSeqStats(Alignment *alnObj)
-{
- FILE *fp = fopen(logfilename.c_str(), "a");
- if (fp == NULL)
- {
- cerr << "couldn't open file " << logfilename << " for logging of stats\n";
- return;
- }
-
- // alignment length is the length of any sequence
- fprintf(fp, "aln len: %d\n", alnObj->getSeqLength(1));
-
-
- std::vector<double> pwIdents;
- double lowestPwId = 1.0;
- double hightestPwId = 0.0;
-
- // create vector of pairwise identities
- for(int s1 = 1; s1 <= alnObj->getNumSeqs(); s1++)
- {
- for(int s2 = s1+1; s2 <= alnObj->getNumSeqs(); s2++)
- {
- double thisPwId = pairwiseIdentity(alnObj, s1, s2);
- pwIdents.push_back(thisPwId);
- if (thisPwId>hightestPwId)
- hightestPwId=thisPwId;
- if (thisPwId<lowestPwId)
- lowestPwId=thisPwId;
- //fprintf(fp, "ident %s:%s %d:%d=%f\n", alnObj->getName(s1).c_str(), alnObj->getName(s2).c_str(), s1, s2, PairwiseIdentity(alnObj, s1, s2));
- }
- }
- fprintf(fp, "aln pw-id highest: %.2f\n", hightestPwId);
- fprintf(fp, "aln pw-id lowest: %.2f\n", lowestPwId);
- fprintf(fp, "aln pw-id avg: %.2f\n", utilityObject->average(pwIdents));
- fprintf(fp, "aln pw-id std-dev: %.2f\n", utilityObject->stdDev(pwIdents));
- fprintf(fp, "aln pw-id median: %.2f\n", utilityObject->median(pwIdents));
-
- fclose(fp);
-}
-
-}
-
git://source.jalview.org / jabaws.git / blobdiff