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
  44         try {
  45         seqName = getSequenceName(flowFile);
  46                 flowFile >> endFlow;
  47         if (!m->control_pressed) {
  48             for(int i=0;i<numFlows;i++) {       flowFile >> flowData[i];        }
  49             updateEndFlow();
  50             translateFlow();
  51             m->gobble(flowFile);
  52                 }
  53
  54                 if(flowFile){   return 1;       }
  55                 else            {       return 0;       }
  56         }
  57         catch(exception& e) {
  58                 m->errorOut(e, "FlowData", "getNext");
  59                 exit(1);
  60         }
  61
  62 }
  63 //********************************************************************************************************************
  64 string FlowData::getSequenceName(ifstream& flowFile) {
  65         try {
  66                 string name = "";
  67
  68         flowFile >> name;
  69
  70                 if (name.length() != 0) {
  71             m->checkName(name);
  72         }else{ m->mothurOut("Error in reading your flowfile, at position " + toString(flowFile.tellg()) + ". Blank name."); m->mothurOutEndLine(); m->control_pressed = true;  }
  73
  74                 return name;
  75         }
  76         catch(exception& e) {
  77                 m->errorOut(e, "FlowData", "getSequenceName");
  78                 exit(1);
  79         }
  80 }
  81
  82 //**********************************************************************************************************************
  83
  84 void FlowData::updateEndFlow(){
  85         try{
  86
  87         if (baseFlow.length() > 4) { return; }
  88
  89                 //int currLength = 0;
  90                 float maxIntensity = (float) maxHomoP + 0.49;
  91
  92                 int deadSpot = 0;
  93
  94                 while(deadSpot < endFlow){
  95                         int signal = 0;
  96                         int noise = 0;
  97
  98                         for(int i=0;i<baseFlow.length();i++){
  99                                 float intensity = flowData[i + deadSpot];
 100                                 if(intensity > signalIntensity){
 101                                         signal++;
 102
 103                                         if(intensity  < noiseIntensity || intensity > maxIntensity){
 104                                                 noise++;
 105                                         }
 106                                 }
 107                         }
 108
 109                         if(noise > 0 || signal == 0){
 110                                 break;
 111                         }
 112
 113                         deadSpot += baseFlow.length();
 114                 }
 115                 endFlow = deadSpot;
 116
 117         }
 118         catch(exception& e) {
 119                 m->errorOut(e, "FlowData", "findDeadSpot");
 120                 exit(1);
 121         }
 122 }
 123
 124 //**********************************************************************************************************************
 125 //TATGCT
 126 //1 0 0 0 0 1
 127 //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
 128 //Becomes TNT
 129 void FlowData::translateFlow(){
 130         try{
 131         sequence = "";
 132         set<char> charInMiddle;
 133         int oldspot = -1;
 134         bool updateOld = false;
 135
 136         for(int i=0;i<endFlow;i++){
 137                         int intensity = (int)(flowData[i] + 0.5);
 138                         char base = baseFlow[i % baseFlow.length()];
 139
 140             if (intensity == 0) { //are we in the middle
 141                 if (oldspot != -1) { charInMiddle.insert(base); }
 142             }else if (intensity >= 1) {
 143                 if (oldspot == -1) { updateOld = true;  }
 144                 else {  //check for bases inbetween two 1's
 145                     if (charInMiddle.count(base) != 0) { //we want to covert to an N
 146                         sequence = sequence.substr(0, oldspot+1);
 147                         sequence += 'N';
 148                     }
 149                     updateOld = true;
 150                     charInMiddle.clear();
 151                 }
 152             }
 153
 154                         for(int j=0;j<intensity;j++){
 155                                 sequence += base;
 156                         }
 157
 158             if (updateOld) { oldspot = sequence.length()-1;  updateOld = false; }
 159                 }
 160
 161                 if(sequence.size() > 4){
 162                         sequence = sequence.substr(4);
 163                 }
 164                 else{
 165                         sequence = "NNNN";
 166                 }
 167     }
 168         catch(exception& e) {
 169                 m->errorOut(e, "FlowData", "translateFlow");
 170                 exit(1);
 171         }
 172 }
 173
 174 //**********************************************************************************************************************
 175
 176 void FlowData::capFlows(int mF){
 177
 178         try{
 179
 180                 maxFlows = mF;
 181                 if(endFlow > maxFlows){ endFlow = maxFlows;     }
 182         translateFlow();
 183
 184         }
 185         catch(exception& e) {
 186                 m->errorOut(e, "FlowData", "capFlows");
 187                 exit(1);
 188         }
 189 }
 190
 191 //**********************************************************************************************************************
 192
 193 bool FlowData::hasGoodHomoP(){
 194
 195         try{
 196
 197         float maxIntensity = (float) maxHomoP + 0.49;
 198
 199         for(int i=0;i<endFlow;i++){
 200             if(flowData[i] > maxIntensity){
 201                 return 0;
 202             }
 203         }
 204                 return 1;
 205         }
 206         catch(exception& e) {
 207                 m->errorOut(e, "FlowData", "hasMinFlows");
 208                 exit(1);
 209         }
 210 }
 211
 212 //**********************************************************************************************************************
 213
 214 bool FlowData::hasMinFlows(int minFlows){
 215
 216         try{
 217                 bool pastMin = 0;
 218                 if(endFlow >= minFlows){        pastMin = 1;    }
 219
 220                 return pastMin;
 221         }
 222         catch(exception& e) {
 223                 m->errorOut(e, "FlowData", "hasMinFlows");
 224                 exit(1);
 225         }
 226 }
 227
 228 //**********************************************************************************************************************
 229
 230 Sequence FlowData::getSequence(){
 231
 232         try{
 233                 return Sequence(seqName, sequence);
 234         }
 235         catch(exception& e) {
 236                 m->errorOut(e, "FlowData", "getSequence");
 237                 exit(1);
 238         }
 239 }
 240
 241 //**********************************************************************************************************************
 242
 243 void FlowData::printFlows(ofstream& outFlowFile){
 244         try{
 245 //      outFlowFile << '>' << seqName << locationString << " length=" << seqLength << " numflows=" << maxFlows << endl;
 246                 outFlowFile << seqName << ' ' << endFlow << ' ' << setprecision(2);
 247
 248                 for(int i=0;i<maxFlows;i++){
 249                         outFlowFile << flowData[i] << ' ';
 250                 }
 251                 outFlowFile << endl;
 252         }
 253         catch(exception& e) {
 254                 m->errorOut(e, "FlowData", "printFlows");
 255                 exit(1);
 256         }
 257 }
 258
 259 //**********************************************************************************************************************
 260
 261 void FlowData::printFlows(ofstream& outFlowFile, string scrapCode){
 262         try{
 263                 outFlowFile << seqName << '|' << scrapCode << ' ' << endFlow << ' ' << setprecision(2);
 264
 265                 for(int i=0;i<numFlows;i++){
 266                         outFlowFile << flowData[i] << ' ';
 267                 }
 268                 outFlowFile << endl;
 269         }
 270         catch(exception& e) {
 271                 m->errorOut(e, "FlowData", "printFlows");
 272                 exit(1);
 273         }
 274 }
 275
 276 //**********************************************************************************************************************
 277
 278 string FlowData::getName(){
 279
 280         try{
 281                 return seqName;
 282         }
 283         catch(exception& e) {
 284                 m->errorOut(e, "FlowData", "getName");
 285                 exit(1);
 286         }
 287 }
 288
 289 //**********************************************************************************************************************