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(ifstream& flowFile, float signal, float noise, int maxHomoP){
  19
  20         try {
  21                 m = MothurOut::getInstance();
  22
  23                 baseFlow = "TACG";
  24                 seqName = "";
  25                 numFlows = 0;
  26                 locationString = "";
  27                 seqLength = 0;
  28
  29                 string lengthString;
  30                 string flowString;
  31
  32                 flowFile >> seqName >> locationString >> lengthString >> flowString;
  33
  34                 convert(lengthString.substr(7), seqLength);
  35                 convert(flowString.substr(9), numFlows);
  36
  37                 flowData.resize(numFlows);
  38
  39                 if (seqName == "") {
  40                         m->mothurOut("Error reading quality file, name blank at position, " + toString(flowFile.tellg()));
  41                         m->mothurOutEndLine();
  42                 }
  43                 else{
  44                         seqName = seqName.substr(1);
  45                         for(int i=0;i<numFlows;i++)     {       flowFile >> flowData[i];        }
  46                 }
  47
  48                 findDeadSpot(signal, noise, maxHomoP);
  49                 translateFlow();
  50
  51                 m->gobble(flowFile);
  52         }
  53         catch(exception& e) {
  54                 m->errorOut(e, "FlowData", "FlowData");
  55                 exit(1);
  56         }
  57
  58 }
  59
  60 //**********************************************************************************************************************
  61
  62 void FlowData::findDeadSpot(float signalIntensity, float noiseIntensity, int maxHomoP){
  63         try{
  64
  65                 int currLength = 0;
  66                 float maxIntensity = (float) maxHomoP + 0.49;
  67
  68                 deadSpot = 0;
  69                 while(currLength < seqLength + 4){
  70                         int signal = 0;
  71                         int noise = 0;
  72
  73                         for(int i=0;i<4;i++){
  74                                 float intensity = flowData[i + 4 * deadSpot];
  75                                 if(intensity > signalIntensity){
  76                                         signal++;
  77
  78                                         if(intensity  < noiseIntensity || intensity > maxIntensity){
  79                                                 noise++;
  80                                         }
  81                                 }
  82                                 currLength += (int)(intensity+0.5);
  83                         }
  84
  85                         if(noise > 0 || signal == 0){
  86                                 break;
  87                         }
  88
  89                         deadSpot++;
  90                 }
  91                 deadSpot *= 4;
  92                 seqLength = currLength;
  93
  94         }
  95         catch(exception& e) {
  96                 m->errorOut(e, "FlowData", "findDeadSpot");
  97                 exit(1);
  98         }
  99 }
 100
 101 //**********************************************************************************************************************
 102
 103 void FlowData::translateFlow(){
 104
 105         try{
 106                 sequence = "";
 107                 for(int i=0;i<deadSpot;i++){
 108                         int intensity = (int)(flowData[i] + 0.5);
 109                         char base = baseFlow[i % 4];
 110
 111                         for(int j=0;j<intensity;j++){
 112                                 sequence += base;
 113                         }
 114                 }
 115
 116                 if(sequence.size() > 4){
 117                         sequence = sequence.substr(4);
 118                 }
 119                 else{
 120                         sequence = "NNNN";
 121                 }
 122         }
 123         catch(exception& e) {
 124                 m->errorOut(e, "FlowData", "translateFlow");
 125                 exit(1);
 126         }
 127 }
 128
 129 //**********************************************************************************************************************
 130
 131 void FlowData::capFlows(int maxFlows){
 132
 133         try{
 134
 135                 numFlows = maxFlows;
 136                 if(deadSpot > maxFlows){        deadSpot = maxFlows;    }
 137
 138         }
 139         catch(exception& e) {
 140                 m->errorOut(e, "FlowData", "capFlows");
 141                 exit(1);
 142         }
 143 }
 144
 145 //**********************************************************************************************************************
 146
 147 bool FlowData::hasMinFlows(int minFlows){
 148
 149         try{
 150                 bool pastMin = 0;
 151
 152                 if(deadSpot >= minFlows){       pastMin = 1;    }
 153                 return pastMin;
 154         }
 155         catch(exception& e) {
 156                 m->errorOut(e, "FlowData", "hasMinFlows");
 157                 exit(1);
 158         }
 159 }
 160
 161 //**********************************************************************************************************************
 162
 163 Sequence FlowData::getSequence(){
 164
 165         try{
 166                 return Sequence(seqName, sequence);
 167         }
 168         catch(exception& e) {
 169                 m->errorOut(e, "FlowData", "getSequence");
 170                 exit(1);
 171         }
 172 }
 173
 174 //**********************************************************************************************************************
 175
 176 int FlowData::getSeqLength(){
 177
 178         try{
 179                 return seqLength;
 180         }
 181         catch(exception& e) {
 182                 m->errorOut(e, "FlowData", "getSeqLength");
 183                 exit(1);
 184         }
 185 }
 186
 187 //**********************************************************************************************************************
 188
 189 void FlowData::printFlows(ofstream& outFlowFile){
 190         try{
 191         //      outFlowFile << '>' << seqName << locationString << " length=" << seqLength << " numflows=" << maxFlows << endl;
 192                 outFlowFile << seqName << ' ' << deadSpot << ' ' << setprecision(2);
 193
 194                 for(int i=0;i<numFlows;i++){
 195                         outFlowFile << flowData[i] << ' ';
 196                 }
 197                 outFlowFile << endl;
 198         }
 199         catch(exception& e) {
 200                 m->errorOut(e, "FlowData", "printFlows");
 201                 exit(1);
 202         }
 203 }
 204
 205 //**********************************************************************************************************************
 206
 207 void FlowData::printFlows(ofstream& outFlowFile, string scrapCode){
 208         try{
 209
 210         //      outFlowFile << '>' << seqName << locationString << " length=" << seqLength << " numflows=" << maxFlows << endl;
 211
 212                 outFlowFile << seqName << '|' << scrapCode << ' ' << deadSpot << ' ' << setprecision(2);
 213
 214                 for(int i=0;i<numFlows;i++){
 215                         outFlowFile << flowData[i] << ' ';
 216                 }
 217                 outFlowFile << endl;
 218         }
 219         catch(exception& e) {
 220                 m->errorOut(e, "FlowData", "printFlows");
 221                 exit(1);
 222         }
 223 }
 224
 225 //**********************************************************************************************************************
 226
 227 void FlowData::printFASTA(ofstream& outFASTA){
 228         try{
 229
 230                 outFASTA << '>' << seqName << endl;
 231                 outFASTA << sequence << endl;
 232
 233         }
 234         catch(exception& e) {
 235                 m->errorOut(e, "FlowData", "printFlows");
 236                 exit(1);
 237         }
 238 }
 239
 240
 241 //**********************************************************************************************************************