二つの遺伝子発現データを比較し，発現領域が近い遺伝子のペアを出力する．#DBCLS

nearest_gene.cc

#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <vector>

using namespace std;

// 遺伝子発現領域データ用構造体
// gene expression domain
struct gene {
  int file_id;
  int left;
  int right;
  int center;
  double fpkm;

  gene() {}
  gene(int id, int l, int r, double f)
    : file_id(id), left(l), right(r), center((l + r) / 2), fpkm(f) {}

  // sort用
  // centerの値を基準にしてsortされるようにする
  bool operator <(const gene &rhs) const {
    if (center  != rhs.center)  { return center  < rhs.center; }
    if (left    != rhs.left)    { return left    < rhs.left; }
    if (right   != rhs.right)   { return right   < rhs.right; }
    if (file_id != rhs.file_id) { return file_id < rhs.file_id; }
    return fpkm < rhs.fpkm;
  }
};

typedef vector<gene> gene_list;

// 遺伝子発現データファイルを読み込む
gene_list read_gene_list_file(const char *filename, int file_id) {
  ifstream fin(filename);
  string s;
  getline(fin, s); // 1行目はヘッダーなので捨てる

  // 2行目以降の各行を処理する
  gene_list genes;
  while (getline(fin, s)) {
    istringstream sin(s);
    string t;
    int l, r;
    double f;
    // タブ区切りの4, 5, 6列目がそれぞれleft, right, fpkm．
    sin >> t >> t >> t >> l >> r >> f;

    // Bのgene listに追加する
    genes.push_back(gene(file_id, l, r, f));
  }
  fin.close();

  return genes;
}

// genesのデータのうちcenterがgのcenterに最も近いものを探す
gene_list::const_iterator find_nearest_gene(const gene &g, const gene_list &genes) {
  // upper: genesのcenterのうち，gのcenter以上のもので最小の要素
  // exampleの場合: 17
  gene_list::const_iterator upper = lower_bound(genes.begin(), genes.end(), g);
  if (upper == genes.end()) { --upper; }

  // lower: genesのcenterのうち，gのcenter以下のもので最大の要素
  // exampleの場合: 10
  gene_list::const_iterator lower = upper;
  if (upper != genes.begin() && g.center < upper->center) { --lower; }

  // lowerとupperのうちgのcenterに近い方を返す
  return ((abs(g.center - lower->center) <
        abs(g.center - upper->center))) ? lower : upper;
}

int main(int argc, char const* argv[]) {
  // 引数チェック
  if (argc < 3) {
    cout << "usage: $ " << argv[0] <<
      " gene_expression_domain_A.tsv gene_expression_domain_B.tsv ..." << endl;
    return 1;
  }

  // for each tsv file
  int input_files = argc - 1;
  const char **filenames = argv + 1;
  gene_list genes[input_files];

  // file id = 0, 1, ..., input_files - 1
  for (int i = 0; i < input_files; ++i) {
    // ファイルの遺伝子発現データをgenesに追加する
    genes[i] = read_gene_list_file(filenames[i], i);
  }

  // gene expression domainのデータをcenterでsortする
  for (int i = 0; i < input_files; ++i) {
    sort(genes[i].begin(), genes[i].end());
  }

  // 各ファイルのgenesごとに最も近い他のファイルのgeneを出力する
  for (int i = 0; i < input_files; ++i) {
    ofstream fout((string(filenames[i]) + ".out").c_str());

    for (gene_list::const_iterator g = genes[i].begin(), end = genes[i].end(); g != end; ++g) {
      int id = g->file_id;
      fout << filenames[i] << "\t" << g->left << "\t" << g->right << "\t" << g->fpkm;

      // 他のファイルのgeneの出力
      for (int j = 0; j < input_files; ++j) {
        if (j != id) {
          gene_list::const_iterator nearest = find_nearest_gene(*g, genes[j]);

          fout << "\t" << filenames[nearest->file_id] << "\t" << nearest->left << "\t" <<
            nearest->right << "\t" << nearest->fpkm << "\t" << abs(g->center - nearest->center);
        }
      }
      fout << endl;
    }
    fout.close();
  }

  return 0;
}