Aritalab:Lecture/Programming/Cpp/Genome

Revision as of 17:13, 19 October 2010

ここではVisual Studio 2010の利用を前提として、ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。

データの読み込み

まず UCSD Downloadサイトに行って、ヒトゲノムのデータをダウンロードしてみます。ここではData set by chromosomeのセクションに進んで日本がシーケンスした21版染色体 (chr21.fa.gz) を取得します。faというのはFastA形式の意味です。wcコマンドにより96万行、49MBのファイルであることがわかります。headコマンドにより各行50文字程度あることがわかります。

$ gunzip chr21.fa.gz
$ wc chr21.fa
  962599   962599 49092500 chr21.fa
$ head chr21.fa
>chr21
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

Visual Studioによるプログラム作成

長さ100文字のバッファを使って行を読み、50MBのヒープ領域に21番染色体のデータを格納しましょう。 Visual Studioインストールしたら、「File→新規作成→プロジェクト」を選択して、新しいプロジェクトを開始します。 プロジェクト名.cppというソースを書くウィンドウが作成されるので、ここにプログラムを書いていきます。 ここでは例として、lectureというプロジェクト名で解説します。 まず、ファイルから1行ずつ読み込んで結合するプログラムを、main関数の中で作成します。

ファイルから1行ずつ読み込むサンプル lecture.cpp

#include "stdafx.h"
#include <iostream>
#include <fstream>

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
  char line[100];
  char* chr = new char[50000000];
  int size=0;
  ifstream fin("chr21.fa");
  if (!fin)
    { cout << "File not found" << endl; exit(0); }
  fin.getline(line, 100);　// 1行目の ">chr21" を捨てる
  while (!fin.getline(line, 100).eof())
  {
    for(int j=0; line[j] != 0; j++)
      chr[size++] = line[j];
  }
　　chr[size] = 0;
  cout << "Read " << size << " characters." << endl;
}

プログラムを書いたら、「ビルド→"プロジェクト名"のビルド」 でコンパイルします。 コンパイルエラーがでたら、エラーメッセージ毎にソース位置を特定できるので、修正してください。 コンパイルできたら、「デバッグ→デバッグ無しで開始」で実行してみましょう。 動作確認ができたら、書いた部分をクラスとして抽象化します。

処理のクラス化

他の染色体を読み込めるよう、

• 染色体のサイズ　(上のプログラムでは50MBの配列chr)
• 実際の文字数　(size)

をまずクラス化します。この段階では、まだ lecture.cpp のファイルに全てのコードが書かれています。

21番染色体を読むサンプル lecture.cpp

#include "stdafx.h"
#include <iostream>
#include <fstream>

using namespace std;

class Chromosome
{
private:
  char* chr;
  int bufsize;
  int size;

public:
  Chromosome(int siz) { chr = new char[siz]; bufsize = siz; }
  ~Chromosome() { delete[] chr; }
  Chromosome(const Chromosome&)
  { cout << "Error: copy constructor is called." << endl; exit(0); }
  Chromosome& operator=(const Chromosome& x)
  { cout << "Error: assignment constructor is called." << endl; exit(0); }

  void open(char* filename)
  {
    char line[100];
    size=0;
    ifstream fin(filename);
    if (!fin)
    { cout << "File not found" << endl; exit(0); }
    // skip the first line
    fin.getline(line, 100);
    while (!fin.getline(line, 100).eof())
    {
      for(int j=0; line[j] != 0; j++) {
        chr[size++] = line[j];
        if (size >= bufsize)
        { cout << "Error: bufsize overflow" << endl; exit(0); }
      }
    }
    chr[size] = 0;
  }

  void printGC()
  {
    int baseA =0, baseT =0, baseG = 0, baseC =0, baseN=0;
    for (int i=0; i < size; i++)
      switch (chr[i])
    { 
      case 'a': case 'A': baseA++; break;
      case 't': case 'T': baseT++; break;
      case 'g': case 'G': baseG++; break;
      case 'c': case 'C': baseC++; break;
      case 'n': case 'N': baseN++; break;
    }
    cout << "GC content: " << (baseG+baseC)*100 / (baseA+baseT) << "%" << endl;
  }
};

int _tmain(int argc, _TCHAR* argv[])
{

  Chromosome C21(50000000);
  C21.open("chr21.fa");
  C21.printGC();
}

これを実行すると、21番染色体のGC含量は69%とわかりました。

ファイルの分割

クラスChromosomeを他のプログラムでも利用できるように別ファイルにしてみましょう。 Javaと同じように、1ファイル1クラスにしておくと便利です。

ファイルを分割するには、ヘッダーファイルとソースファイルを追加します。 Visual Studioの画面左側にあるソリューションエクスプローラーをクリックします。

1. ヘッダーファイルフォルダを右クリックして、「追加→新しい項目→ヘッダーファイル」を選びます。ファイル名は chromosome.h とします。
2. Chromosomeクラスの定義を、そっくりそのまま chromosome.h の中に移します。先頭のインクルード文と、using namespace std;という行も移します。
3. プロジェクト名.cpp のヘッダー部分には、代わりに #include "chromosome.h" と書きます。

この時点で、一度再コンパイルしてみてください。問題なくコンパイルできるはずです。 ファイルを#includeすることは、その部分にファイルの中身をそのまま展開することにあたるので、当たり前といえば当たり前です。

しかし、このままではChromosomeクラスの中身が全てヘッダーに入ったままなので、一部をcppファイルとして置き換えましょう。再びソリューションエクスプローラーをクリックします。

1. ソースファイルフォルダを右クリックして、「追加→新しい項目→C++ソース」を選びます。ファイル名は chromosome.cpp とします。
2. chromosome.cpp ファイルにヘッダー　#include "chromosome.h" を書きます。

これからヘッダーファイル内のプログラムを順次、ソースファイル内に移します。 原則は以下の通りです。

• 処理コードが短いコンストラクタ等はヘッダーに残す。
• 処理コードが長いメンバー関数は、ヘッダー部分に以下の型情報だけを書く。

戻り値　関数名( 引数 );

• メンバー関数の中身は、ソースファイルに以下の形式で書く。

 戻り値 クラス名::関数名( 引数 )
  { 関数定義 }

具体的には、以下のようになります。

chromosome.h サンプル

#include <iostream>
using namespace std;

class Chromosome
{
private:
  char* chr;
  int bufsize;
  int size;

public:
  Chromosome(int siz) { chr = new char[siz]; bufsize = siz; }
  ~Chromosome() { delete[] chr; }
  Chromosome(const Chromosome&)
  { cout << "Error: copy constructor is called." << endl; exit(0); }
  Chromosome& operator=(const Chromosome& x)
  { cout << "Error: assignment constructor is called." << endl; exit(0); }

  void open(char *);
  void printGC();
};

chromosome.cpp サンプル

#include "stdafx.h"
#include <iostream>
#include <fstream>
#include "chromosome.h"
using namespace std;

void Chromosome::open(char* filename)
  { ..関数定義.. }

void Chromosome::printGC()
  { ..関数定義.. }

lecture.cpp サンプル

#include "stdafx.h"
#include "chromosome.h"
using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
  Chromosome C21(50000000);
  C21.open("chr21.fa");
  C21.printGC();
}

これでファイルの分割は完了です。先頭に書く#include文は、そのファイルの中で利用する情報だけを書けば十分です。

クラスメンバー関数の充実

ここまでできれば、後はクラスのメンバー関数を順次充実させられます。 例えば、ゲノムの文字は大文字と小文字が混在しているので、全て小文字に変換する関数 toLowercase を用意してみましょう。

chromosome.cpp への追加コード

char Chromosome::toLowercase(char c)
  {
    switch (c) {
    case 'A': return 'a';
    case 'T': return 't';
    case 'G': return 'g';
    case 'C': return 'c';
    case 'N': return 'n';
    case 'a': case 't': case 'g': case 'c': case 'n':
      return c;
    default:
      cout << "Unknown char " << c << endl;
      exit(0);
    }
  }

ヘッダーファイル chromosome.h には、 char toLowerase(char); の１行を追加し、メンバー関数 open(char* filename) の中で

としておきましょう。