source: plugins/base2/net.sf.basedb.normalizers/trunk/src/c++/bin/qQN.cc @ 1220

// $Id: qQN.cc 1220 2010-04-07 14:38:43Z jari $

/*
  Copyright (C) 2009 Jari Häkkinen

  This file is part of the Normalizers plug-in package for BASE
  (net.sf.based.normalizers). The package is available at
  http://baseplugins.thep.lu.se/ BASE main site is
  http://base.thep.lu.se/

  This is free software; you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  The software is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h> // this header file is created by configure

#include <yat/normalizer/ColumnNormalizer.h>
#include <yat/normalizer/qQuantileNormalizer.h>

#include <yat/utility/CommandLine.h>
#include <yat/utility/MatrixWeighted.h>
#include <yat/utility/OptionHelp.h>
#include <yat/utility/OptionInFile.h>
#include <yat/utility/OptionOutFile.h>
#include <yat/utility/OptionSwitch.h>
#include <yat/utility/stl_utility.h>
#include <yat/statistics/utility.h>

#include <cstdlib>
#include <fstream>
#include <functional>
#include <iostream>
#include <limits>
#include <stdexcept>

using namespace theplu::yat::normalizer;
using namespace theplu::yat::utility;


class CleanUpMatrix : std::unary_function<DataWeight, DataWeight>
{
public:
  CleanUpMatrix(void) {}

  inline DataWeight operator()(DataWeight x) const
    { return ( x.data()>0 ?
               x : DataWeight(std::numeric_limits<double>::quiet_NaN(),0.0) ); }
};


void create_matrix_until_183_fixed(MatrixWeighted**, const MatrixWeighted&,
                                   const std::string&);
void create_target(std::vector<double>&, const MatrixWeighted&, bool use_median);
void create_target(std::vector<double>&, const MatrixWeighted&,
                   const std::string&, bool use_median);

/**
   writes the data values in the matrix ignoring the weights, i.e.,
   produces the same output as the Matrix output operator does.
 */
std::ostream& operator<< (std::ostream&, const MatrixWeighted&);


int main(int argc, char* argv[])
{
  CommandLine cmd;
  OptionInFile assay(cmd, "assay-data", "assay annotations");
  OptionInFile indata(cmd, "in-data", "data to be normalized");
  OptionOutFile outdata(cmd, "out-data", "normalized data");
  OptionSwitch target_median(cmd, "median",
                             "use median instead of mean for targets", false);
  OptionHelp help(cmd);
  help.synopsis()=(std::string("See ") +
                   "http://baseplugins.thep.lu.se/net.sf.basedb.normalizers " +
                   "for\ndetails on this program.\n\n" +
                   "All assays are used to calculate the target if no assay " +
                   "annotion file is set.\n" +
                   "The assay file format is pair information: " +
                   "assay_no\\tyes/no\n" +
                   "where assay_no is the assay number starting at 1, yes " +
                   "indicates that it should\n" +
                   "be used in target calculation, and no indicates that it " +
                   "should not be a part\n" +
                   "of target.\n");
  OptionSwitch version(cmd, "version", "output version and exit");
  std::stringstream copyright;
  copyright << PACKAGE_STRING << '\n'
            << "Copyright (C) 2009 Jari Häkkinen\n\n"
            << "This is free software see the source for copying "
            << "conditions. There is NO\nwarranty; not even for "
            << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\n";
  try {
    cmd.parse(argc, argv);
  }
  catch (cmd_error& e) {
    if (version.present()) {
      std::cout << copyright.str();
      return EXIT_SUCCESS;
    }
    std::cerr << e.what() << std::endl;
    return EXIT_FAILURE;
  }
  if (version.present()) {
    std::cout << copyright.str();
    return EXIT_SUCCESS;
  }
  std::ifstream* infile=NULL;
  std::streambuf* cin_buffer=NULL;
  if (indata.present()) {
    infile=new std::ifstream(indata.value().c_str());
    cin_buffer = std::cin.rdbuf(); // save cin's input buffer
    std::cin.rdbuf(infile->rdbuf());
  }

  // The BASE core API deprecated exportPlainMatrix method of the
  // deprecated BioAssaySetExporter does not behave well when
  // exporting bioassaysets where assays have been filtered out in
  // some filter step. Variables and functions with 183 in their name
  // should be removed when http://baseplugins.thep.lu.se/ticket/183
  // is fixed.

  // m is a pointer below, should be changed to a normal variable when
  // http://baseplugins.thep.lu.se/ticket/183 is fixed.

  // Change to m when http://baseplugins.thep.lu.se/ticket/183 fixed.
  MatrixWeighted m_tmp_until_183_fixed(std::cin,'\t');
  if (indata.present()) {
    std::cin.rdbuf(cin_buffer); // restore old input buffer
    infile->close();
    delete infile;
  }
  // Remove until HERE below when
  // http://baseplugins.thep.lu.se/ticket/183 fixed.
  MatrixWeighted* m=NULL;
  if (assay.present())
    create_matrix_until_183_fixed(&m, m_tmp_until_183_fixed, assay.value());
  else
    m = &m_tmp_until_183_fixed;
  // HERE remove until here.

  std::transform(m->begin(), m->end(), m->begin(), CleanUpMatrix());
  std::vector<double> target;
  ( assay.present() ?
    create_target(target, *m, assay.value(), target_median.value()) :
    create_target(target, *m, target_median.value()) );
  std::transform(target.begin(), target.end(),
                 target.begin(), theplu::yat::utility::Log<double>());
  std::transform(data_iterator(m->begin()), data_iterator(m->end()),
                 data_iterator(m->begin()), theplu::yat::utility::Log<double>());
  // q = min(100,target_size/10) but no smaller than 10
  unsigned int q=target.size()/10;
  q = q>100 ? 100 : ( q<10 ? 10 : q );
  qQuantileNormalizer qqn(target.begin(), target.end(), q);
  ColumnNormalizer<qQuantileNormalizer> cn(qqn);
  MatrixWeighted result(m->rows(),m->columns());
  cn(*m,result);
  std::transform(data_iterator(result.begin()), data_iterator(result.end()),
                 data_iterator(result.begin()), theplu::yat::utility::Exp<double>());

  std::ofstream* outfile=NULL;
  std::streambuf* cout_buffer = std::cout.rdbuf();
  if (outdata.present()) {
    outfile=new std::ofstream(outdata.value().c_str());
    cout_buffer = std::cout.rdbuf(); // save cout's output buffer
    std::cout.rdbuf(outfile->rdbuf());
  }
  std::cout << result << std::endl;
  if (outdata.present()) {
    std::cout.rdbuf(cout_buffer); // restore old output buffer
    outfile->close();
    delete outfile;
  }

  // clean up until http://baseplugins.thep.lu.se/ticket/183 fixed.
  if (assay.present())
    delete m;

  return EXIT_SUCCESS;
}


void create_matrix_until_183_fixed(MatrixWeighted** m,
                                   const MatrixWeighted& m_tmp_183,
                                   const std::string& assay)
{
  std::ifstream is(assay.c_str());
  std::vector<size_t> idx;
  while (is.good()) {
    size_t i;
    is >> i;
    if (is.good()) idx.push_back(i);
    std::string line;
    getline(is, line);
    if (idx.size()>m_tmp_183.columns())
      throw std::runtime_error(std::string("183_tmp: Too many annotation ") +
                               "columns wrt data matrix");
  }

  if (*m) delete *m;
  size_t rows=m_tmp_183.rows();
  size_t cols=idx.size();
  *m=new MatrixWeighted(rows,cols);
  for (size_t i=0; i<rows; ++i)
    for (size_t j=0; j<cols; ++j)
      (**m)(i,j)=m_tmp_183(i,idx[j]-1);
}


void create_target(std::vector<double>& t, const MatrixWeighted& m,
                   const std::string& assay, bool use_median)
{
  std::vector< std::vector<double> > calc_support(m.rows());
  for (std::vector< std::vector<double> >::iterator i=calc_support.begin();
       i!=calc_support.end(); ++i)
    i->reserve(m.columns());

  std::ifstream is(assay.c_str());
  std::string line;
  size_t column=0;
  while (getline(is, line)) {
    size_t found=line.find("yes");
    if (found!=std::string::npos) // use this assay as a part of reference
      for (size_t row=0; row<m.rows(); ++row)
        if (m(row,column).weight())       // weight either 0 or 1
          calc_support[row].push_back(m(row,column).data());
    ++column;
    if (column>m.columns())
      throw std::runtime_error("Too many annotation columns wrt data matrix");
  }

  t.reserve(m.rows());
  for (size_t row=0; row<m.rows(); ++row) {
    const std::vector<double>& cs=calc_support[row];
    if (!cs.empty())
      t.push_back( use_median ?
                   theplu::yat::statistics::median(cs.begin(), cs.end()) :
                   std::accumulate(cs.begin(), cs.end(), 0.0) / cs.size() );
  }
  if (!t.size())
    throw std::runtime_error("Not a well defined reference, aborting");
}


void create_target(std::vector<double>& t, const MatrixWeighted& m,
                   bool use_median)
{
  t.reserve(m.rows());
  std::vector<double> cs;
  cs.reserve(m.columns());
  for (size_t row=0; row<m.rows(); ++row) {
    cs.clear();
    for (size_t column=0; column<m.columns(); ++column)
      if (m(row,column).weight())      // weight either 0 or 1
        cs.push_back(m(row,column).data());
    if (!cs.empty())
      t.push_back( use_median ?
                   theplu::yat::statistics::median(cs.begin(), cs.end()) :
                   std::accumulate(cs.begin(), cs.end(), 0.0) / cs.size() );
  }
  if (!t.size())
    throw std::runtime_error("Not a well defined reference, aborting");
}


std::ostream& operator<< (std::ostream& s, const MatrixWeighted& m)
{
  s.setf(std::ios::dec);
  s.precision(12);
  for(size_t i=0, j=0; i<m.rows(); i++)
    for (j=0; j<m.columns(); j++) {
      s << m(i,j).data();
      if (j<m.columns()-1)
        s << s.fill();
      else if (i<m.rows()-1)
        s << "\n";
    }
  return s;
}