text: /*
   * window.cpp for MSIsensor1.1
   * Copyright (c) 2013 Beifang Niu && Kai Ye CNIC && Xi’an Jiaotong University All Rights Reserved.
   *
   * Permission is hereby granted, free of charge, to any person
   * obtaining a copy of this software and associated documentation
   * files (the "Software"), to deal in the Software without
   * restriction, including without limitation the rights to use,
   * copy, modify, merge, publish, distribute, sublicense, and/or sell
   * copies of the Software, and to permit persons to whom the
   * Software is furnished to do so, subject to the following
   * conditions:
   *
   * The above copyright notice and this permission notice shall be
   * included in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
   * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
   * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
   * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
   * OTHER DEALINGS IN THE SOFTWARE.
   */
  
  #include <iostream>
  #include <sstream>
  #include <bitset>
  #include <omp.h>
  
  #include "window.h"
  #include "param.h"
  #include "polyscan.h"
  
  extern Param paramd;
  extern PolyScan polyscan;
  extern char homo_code[];
  extern char uhomo_code[];
  extern bit8_t alphabet[];
  
  Window::Window() 
      : _start(0)
      , _end(0)
      , _chr("")
      , _siteCount(0)
      , _startSite(NULL)
      , _endSite(NULL)
  {
      //xxxx
      //
  };
  
  Window::~Window() {
      //xxx
  };
  
  void Window::InitialDisW() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->InitialDis();
      }
  };
  
  
  void Window::InitialTumorDisW() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->InitialTumorDis();
      }
  };
  
  void Window::ClearDis() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->ReleaseMemory();
      }
  };
  
  
  void Window::ClearTumorDis() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->ReleaseTumorMemory();
      }
  };
  
  void Window::OutputDisW() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->OutputDis();
      }
  };
  
  void Window::OutputTumorDisW() {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->OutputTumorDis();
      }
  };
  
  void Window::PouroutDisW(Sample &oneSample) {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->PouroutDis(oneSample);
      }
  };
  
  void Window::PourTumoroutDisW(Sample &oneSample) {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->PourTumoroutDis(oneSample);
      }
  };
  
  void Window::DisGenotypingW(Sample &oneSample) {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++) {
          p = _startSite + i;
          p->DisGenotyping(oneSample);
      }
  };
  
  void Window:: PouroutTumorSomatic(Sample &oneSample) {
      HomoSite *p = NULL;
      for (unsigned short i=0; i<_siteCount; i++){
          p = _startSite + i;
          p->DisTumorSomatic(oneSample);
      }
  };
  
  // change start
  void Window::ChangeStart() {
      if ((_start - MAX_SPAN_SIZE) < 0) {
          _start = 0;
      } else { _start -= MAX_SPAN_SIZE; }
  }
  
  void Window::GetDistribution(std::vector <SPLIT_READ> &readsInWindow) {
      for (unsigned short j=0; j<polyscan.totalBamPairsNum; j++) {
          // normal
          if (!polyscan.totalBamPairs[j].normal_bam.empty()) {
              // extract reads
              LoadReads(readsInWindow, polyscan.totalBamPairs[j].normal_bam.c_str());
              ScanReads(readsInWindow, j, false);
              readsInWindow.clear();
          }
          // tumor
          if (!polyscan.totalBamPairs[j].tumor_bam.empty()) {
              // extract reads
              LoadReads(readsInWindow, polyscan.totalBamPairs[j].tumor_bam.c_str());
              ScanReads(readsInWindow, j, true);
              readsInWindow.clear();
          }
      }
  }
  
  void Window::GetTumorDistribution(std::vector <SPLIT_READ> &readsInWindow) {
      for (unsigned short j=0; j<polyscan.totalBamTumorsNum; j++) {
          // tumor
          if (!polyscan.totalBamTumors[j].tumor_bam.empty()) {
              // extract reads
              LoadReads(readsInWindow, polyscan.totalBamTumors[j].tumor_bam.c_str());
              ScanReads(readsInWindow, j, true);
              readsInWindow.clear();
          }
      }
  }
  void Window::LoadReads( std::vector <SPLIT_READ> &readsInWindow, 
                          const std::string bam ) {
      std::string tag = "";
      if (!bam.empty()) {
          // extract reads
          ReadInBamReads(bam.c_str(), _chr, _start, _end, readsInWindow, tag);
      }
  }
  
  void Window::ScanReads( const std::vector <SPLIT_READ> &readsInWindow, 
                          unsigned short bamIndex, 
                          bool isTumor) {
  
      // openmp parallel 
      omp_set_num_threads( paramd.numberThreads );
  #pragma omp parallel for
      for (unsigned short i=0; i<_siteCount; i++) {
          HomoSite *p = _startSite + i;
          unsigned long tsize = readsInWindow.size();
          for (unsigned long j=0; j<tsize; j++) {
              if ( readsInWindow[j].Mapped ) {
                  if ( (readsInWindow[j].MatchedRelPos < p->lowcut) || (readsInWindow[j].MatchedRelPos > p->highcut) ) continue; 
              }
              unsigned short tCount = DoOneRead(readsInWindow[j].ReadSeq, p);
              if ( (tCount > 0) && (tCount < paramd.s_dispots) ) {
                  if (isTumor) {
                      p->tumorDis[bamIndex][tCount-1]++;
                  } else {
                      p->normalDis[bamIndex][tCount-1]++;
                  }
              } else {
                  // don't scan reverse if mapped
                  if ( readsInWindow[j].Mapped ) continue; 
                  // reverse
                  std::string tStr = readsInWindow[j].ReadSeq;
                  ReverseComplement(tStr);
                  tCount = DoOneRead(tStr, p);
                  if ( (tCount > 0) && (tCount < paramd.s_dispots) ) {
                      if (isTumor) {
                          p->tumorDis[bamIndex][tCount-1]++;
                      } else {
                          p->normalDis[bamIndex][tCount-1]++;
                      }
                  }
              }
          }
      }
  }
  
  unsigned short Window::DoOneRead( const std::string &oneRead, const HomoSite *p ) {
      std::string::size_type startPos = 0;
      unsigned short count = 0;
      while (std::string::npos != (startPos = oneRead.find(p->fbases, startPos))) {
          //std::cout<<startPos<<std::endl;
          count = 0;
          std::string::size_type tstart0 = startPos+p->fbases.length();
          std::string::size_type tstart  = tstart0;
          while ( tstart0 == (tstart = oneRead.find(p->bases, tstart)) ) {
              count++;
              tstart += p->bases.length();
              tstart0 = tstart;
          }
          // if get one
          if ( (p->typeLen == 1) && (count >= paramd.MininalHomoSize)
               ||
               (p->typeLen >  1) && (count >= paramd.MinMicrosate)
             ) {
              tstart = tstart0;
              if (tstart == (tstart0 = oneRead.find(p->ebases, tstart0))) {
                  return count;
              }
          }
          startPos++;
      }
      return 0;
  }
  
  // reverse complement string
  void Window::ReverseComplement(std::string &theWord) {       
      char tempChar;
      unsigned int t_uint0;
      unsigned int t_uint1;    
      for (int i=0; i<theWord.length()/2; i++) {
          tempChar = theWord[i];
          t_uint0 = alphabet[tempChar];
          t_uint1 = alphabet[theWord[theWord.length()-i-1]];
          theWord[i] = uhomo_code[t_uint1];
          theWord[theWord.length()-i-1] = uhomo_code[t_uint0];
      }
      if (theWord.length()%2) {
          tempChar = theWord[theWord.length()/2];
          t_uint0 = alphabet[tempChar];
          theWord[theWord.length()/2] = uhomo_code[t_uint0];
      }
  }  
  
  