slayer.cpp

   1 /*
   2  *  slayer.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/25/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "slayer.h"
  11
  12 /***********************************************************************/
  13 Slayer::Slayer(int win, int increment, int parentThreshold, float div) :
  14                 windowSize(win), windowStep(increment), parentFragmentThreshold(parentThreshold), divRThreshold(div) {}
  15 /***********************************************************************/
  16 void Slayer::getResults(Sequence* query, vector<Sequence*> refSeqs) {
  17         try {
  18
  19                 for (int i = 0; i < refSeqs.size(); i++) {
  20
  21                         for (int j = i+1; j < refSeqs.size(); j++) {
  22
  23                                 //make copies of query and each parent because runBellerophon removes gaps and messes them up
  24                                 Sequence* q = new Sequence(query->getName(), query->getAligned());
  25                                 Sequence* leftParent = new Sequence(refSeqs[i]->getName(), refSeqs[i]->getAligned());
  26                                 Sequence* rightParent = new Sequence(refSeqs[j]->getName(), refSeqs[j]->getAligned());
  27
  28                                 vector<data_struct> divs = runBellerophon(q, leftParent, rightParent);
  29
  30                                 delete q;
  31                                 delete leftParent;
  32                                 delete rightParent;
  33                         }
  34
  35                 }
  36
  37
  38         }
  39         catch(exception& e) {
  40                 errorOut(e, "Slayer", "getResults");
  41                 exit(1);
  42         }
  43 }
  44 /***********************************************************************/
  45 vector<data_struct> Slayer::runBellerophon(Sequence* query, Sequence* parentA, Sequence* parentB) {
  46         try{
  47
  48                 vector<data_struct> data;
  49
  50                 //vertical filter
  51                 vector<Sequence*> temp;
  52                 verticalFilter(temp);
  53
  54                 int alignLength = query->getAligned().length();
  55
  56
  57
  58
  59                 return data;
  60
  61         }
  62         catch(exception& e) {
  63                 errorOut(e, "Slayer", "runBellerophon");
  64                 exit(1);
  65         }
  66 }
  67 /***********************************************************************/
  68 float Slayer::computePercentID(string queryFrag, string parent, int left, int right) {
  69         try {
  70                 int total = 0;
  71                 int matches = 0;
  72
  73                 for (int i = left; i <= right; i++) {
  74                         total++;
  75                         if (queryFrag[i] == parent[i]) {
  76                                 matches++;
  77                         }
  78                 }
  79
  80                 float percentID =( matches/(float)total) * 100;
  81
  82                 return percentID;
  83         }
  84         catch(exception& e) {
  85                 errorOut(e, "Slayer", "computePercentID");
  86                 exit(1);
  87         }
  88 }
  89 /***********************************************************************/
  90 //this is a vertical filter
  91 void Slayer::verticalFilter(vector<Sequence*> seqs) {
  92         try {
  93                 vector<int> gaps;       gaps.resize(seqs[0]->getAligned().length(), 0);
  94
  95                 string filterString = (string(seqs[0]->getAligned().length(), '1'));
  96
  97                 //for each sequence
  98                 for (int i = 0; i < seqs.size(); i++) {
  99
 100                         string seqAligned = seqs[i]->getAligned();
 101
 102                         for (int j = 0; j < seqAligned.length(); j++) {
 103                                 //if this spot is a gap
 104                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
 105                         }
 106                 }
 107
 108                 //zero out spot where all sequences have blanks
 109                 int numColRemoved = 0;
 110                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 111                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
 112                 }
 113
 114                 //for each sequence
 115                 for (int i = 0; i < seqs.size(); i++) {
 116
 117                         string seqAligned = seqs[i]->getAligned();
 118                         string newAligned = "";
 119
 120                         for (int j = 0; j < seqAligned.length(); j++) {
 121                                 //if this spot is not a gap
 122                                 if (filterString[j] == '1') { newAligned += seqAligned[j]; }
 123                         }
 124
 125                         seqs[i]->setAligned(newAligned);
 126                 }
 127         }
 128         catch(exception& e) {
 129                 errorOut(e, "Slayer", "verticalFilter");
 130                 exit(1);
 131         }
 132 }
 133 /***********************************************************************/