--- /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