flowdata.cpp

   1 /*
   2  *  flowdata.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 12/22/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "flowdata.h"
  11
  12 //**********************************************************************************************************************
  13
  14 FlowData::FlowData(){}
  15
  16 //**********************************************************************************************************************
  17
  18 FlowData::~FlowData(){  /*      do nothing      */      }
  19
  20 //**********************************************************************************************************************
  21
  22 FlowData::FlowData(int numFlows, float signal, float noise, int maxHomoP, string baseFlow) :
  23                         numFlows(numFlows), signalIntensity(signal), noiseIntensity(noise), maxHomoP(maxHomoP), baseFlow(baseFlow){
  24
  25         try {
  26                 m = MothurOut::getInstance();
  27
  28                 flowData.assign(numFlows, 0);
  29 //              baseFlow = "TACG";
  30                 seqName = "";
  31                 locationString = "";
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "FlowData", "FlowData");
  35                 exit(1);
  36         }
  37
  38 }
  39
  40 //**********************************************************************************************************************
  41
  42 bool FlowData::getNext(ifstream& flowFile){
  43         try {
  44
  45         seqName = getSequenceName(flowFile);
  46         if (m->debug) {  m->mothurOut("[DEBUG]: flow = " + seqName + " "); }
  47                 flowFile >> endFlow;
  48         if (m->debug) {  m->mothurOut(toString(endFlow) + " "); }
  49         if (!m->control_pressed) {
  50             if (m->debug) {  m->mothurOut(" "); }
  51             for(int i=0;i<numFlows;i++) {
  52                 flowFile >> flowData[i];
  53                 if (m->debug) {  m->mothurOut(toString(flowData[i]) + " "); }
  54             }
  55             if (m->debug) {  m->mothurOut("\n"); }
  56             updateEndFlow();
  57             translateFlow();
  58             m->gobble(flowFile);
  59                 }
  60
  61                 if(flowFile){   return 1;       }
  62                 else            {       return 0;       }
  63         }
  64         catch(exception& e) {
  65                 m->errorOut(e, "FlowData", "getNext");
  66                 exit(1);
  67         }
  68
  69 }
  70 //********************************************************************************************************************
  71 string FlowData::getSequenceName(ifstream& flowFile) {
  72         try {
  73                 string name = "";
  74
  75         flowFile >> name;
  76
  77                 if (name.length() != 0) {
  78             m->checkName(name);
  79         }else{ m->mothurOut("Error in reading your flowfile, at position " + toString(flowFile.tellg()) + ". Blank name."); m->mothurOutEndLine(); m->control_pressed = true;  }
  80
  81                 return name;
  82         }
  83         catch(exception& e) {
  84                 m->errorOut(e, "FlowData", "getSequenceName");
  85                 exit(1);
  86         }
  87 }
  88
  89 //**********************************************************************************************************************
  90
  91 void FlowData::updateEndFlow(){
  92         try{
  93
  94         if (baseFlow.length() > 4) { return; }
  95
  96                 //int currLength = 0;
  97                 float maxIntensity = (float) maxHomoP + 0.49;
  98
  99                 int deadSpot = 0;
 100
 101                 while(deadSpot < endFlow){
 102                         int signal = 0;
 103                         int noise = 0;
 104
 105                         for(int i=0;i<baseFlow.length();i++){
 106                                 float intensity = flowData[i + deadSpot];
 107                                 if(intensity > signalIntensity){
 108                                         signal++;
 109
 110                                         if(intensity  < noiseIntensity || intensity > maxIntensity){
 111                                                 noise++;
 112                                         }
 113                                 }
 114                         }
 115
 116                         if(noise > 0 || signal == 0){
 117                                 break;
 118                         }
 119
 120                         deadSpot += baseFlow.length();
 121                 }
 122                 endFlow = deadSpot;
 123
 124         }
 125         catch(exception& e) {
 126                 m->errorOut(e, "FlowData", "findDeadSpot");
 127                 exit(1);
 128         }
 129 }
 130
 131 //**********************************************************************************************************************
 132 //TATGCT
 133 //1 0 0 0 0 1
 134 //then the second positive flow is for a T, but you saw a T between the last and previous flow adn it wasn't positive, so something is missing
 135 //Becomes TNT
 136 void FlowData::translateFlow(){
 137         try{
 138         sequence = "";
 139         set<char> charInMiddle;
 140         int oldspot = -1;
 141         bool updateOld = false;
 142
 143         for(int i=0;i<endFlow;i++){
 144                         int intensity = (int)(flowData[i] + 0.5);
 145                         char base = baseFlow[i % baseFlow.length()];
 146
 147             if (intensity == 0) { //are we in the middle
 148                 if (oldspot != -1) { charInMiddle.insert(base); }
 149             }else if (intensity >= 1) {
 150                 if (oldspot == -1) { updateOld = true;  }
 151                 else {  //check for bases inbetween two 1's
 152                     if (charInMiddle.count(base) != 0) { //we want to covert to an N
 153                         sequence = sequence.substr(0, oldspot+1);
 154                         sequence += 'N';
 155                     }
 156                     updateOld = true;
 157                     charInMiddle.clear();
 158                 }
 159             }
 160
 161                         for(int j=0;j<intensity;j++){
 162                                 sequence += base;
 163                         }
 164
 165             if (updateOld) { oldspot = sequence.length()-1;  updateOld = false; }
 166                 }
 167
 168                 if(sequence.size() > 4){
 169                         sequence = sequence.substr(4);
 170                 }
 171                 else{
 172                         sequence = "NNNN";
 173                 }
 174     }
 175         catch(exception& e) {
 176                 m->errorOut(e, "FlowData", "translateFlow");
 177                 exit(1);
 178         }
 179 }
 180
 181 //**********************************************************************************************************************
 182
 183 void FlowData::capFlows(int mF){
 184
 185         try{
 186
 187                 maxFlows = mF;
 188                 if(endFlow > maxFlows){ endFlow = maxFlows;     }
 189         translateFlow();
 190
 191         }
 192         catch(exception& e) {
 193                 m->errorOut(e, "FlowData", "capFlows");
 194                 exit(1);
 195         }
 196 }
 197
 198 //**********************************************************************************************************************
 199
 200 bool FlowData::hasGoodHomoP(){
 201
 202         try{
 203
 204         float maxIntensity = (float) maxHomoP + 0.49;
 205
 206         for(int i=0;i<endFlow;i++){
 207             if(flowData[i] > maxIntensity){
 208                 return 0;
 209             }
 210         }
 211                 return 1;
 212         }
 213         catch(exception& e) {
 214                 m->errorOut(e, "FlowData", "hasMinFlows");
 215                 exit(1);
 216         }
 217 }
 218
 219 //**********************************************************************************************************************
 220
 221 bool FlowData::hasMinFlows(int minFlows){
 222
 223         try{
 224                 bool pastMin = 0;
 225                 if(endFlow >= minFlows){        pastMin = 1;    }
 226
 227                 return pastMin;
 228         }
 229         catch(exception& e) {
 230                 m->errorOut(e, "FlowData", "hasMinFlows");
 231                 exit(1);
 232         }
 233 }
 234
 235 //**********************************************************************************************************************
 236
 237 Sequence FlowData::getSequence(){
 238
 239         try{
 240                 return Sequence(seqName, sequence);
 241         }
 242         catch(exception& e) {
 243                 m->errorOut(e, "FlowData", "getSequence");
 244                 exit(1);
 245         }
 246 }
 247
 248 //**********************************************************************************************************************
 249
 250 void FlowData::printFlows(ofstream& outFlowFile){
 251         try{
 252 //      outFlowFile << '>' << seqName << locationString << " length=" << seqLength << " numflows=" << maxFlows << endl;
 253                 outFlowFile << seqName << ' ' << endFlow << ' ' << setprecision(2);
 254
 255                 for(int i=0;i<maxFlows;i++){
 256                         outFlowFile << flowData[i] << ' ';
 257                 }
 258                 outFlowFile << endl;
 259         }
 260         catch(exception& e) {
 261                 m->errorOut(e, "FlowData", "printFlows");
 262                 exit(1);
 263         }
 264 }
 265
 266 //**********************************************************************************************************************
 267
 268 void FlowData::printFlows(ofstream& outFlowFile, string scrapCode){
 269         try{
 270                 outFlowFile << seqName << '|' << scrapCode << ' ' << endFlow << ' ' << setprecision(2);
 271
 272                 for(int i=0;i<numFlows;i++){
 273                         outFlowFile << flowData[i] << ' ';
 274                 }
 275                 outFlowFile << endl;
 276         }
 277         catch(exception& e) {
 278                 m->errorOut(e, "FlowData", "printFlows");
 279                 exit(1);
 280         }
 281 }
 282
 283 //**********************************************************************************************************************
 284
 285 string FlowData::getName(){
 286
 287         try{
 288                 return seqName;
 289         }
 290         catch(exception& e) {
 291                 m->errorOut(e, "FlowData", "getName");
 292                 exit(1);
 293         }
 294 }
 295
 296 //**********************************************************************************************************************