Aritalab:Lecture/Programming/Cpp/Genome

From Metabolomics.JP
< Aritalab:Lecture | Programming | Cpp(Difference between revisions)
Jump to: navigation, search
(New page: ==ゲノムデータを解析する== ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。 ===データのダウンロード=== まず [ht...)
 
m
Line 2: Line 2:
 
ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。
 
ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。
 
===データのダウンロード===
 
===データのダウンロード===
まず [http://hgdownload.cse.ucsc.edu/downloads.html UCSD Download]サイトに行って、ヒトゲノムのデータをダウンロードしてみます。ここではData set by chromosomeのセクションに進んで日本がシーケンスした21版染色体 (chr21.fa.gz) を取得します。faというのはFastA形式の意味です。96万行もあるファイルで、headコマンドにより各行50文字程度あることがわかります。
+
まず [http://hgdownload.cse.ucsc.edu/downloads.html UCSD Download]サイトに行って、ヒトゲノムのデータをダウンロードしてみます。ここではData set by chromosomeのセクションに進んで日本がシーケンスした21版染色体 (chr21.fa.gz) を取得します。faというのはFastA形式の意味です。wcコマンドにより96万行、49MBのファイルであることがわかります。headコマンドにより各行50文字程度あることがわかります。
 
<pre>
 
<pre>
 
$ gunzip chr21.fa.gz
 
$ gunzip chr21.fa.gz
Line 21: Line 21:
  
 
===Visual Studioによるプログラム作成===
 
===Visual Studioによるプログラム作成===
インストールしたら、File&rarr;新規作成&rarr;プロジェクトを選択して、新しいプロジェクトを開始します。
+
長さ100文字のバッファを使って行を読み、50MBのヒープ領域に21番染色体のデータを格納しましょう。
 +
Visual Studioインストールしたら、File&rarr;新規作成&rarr;プロジェクトを選択して、新しいプロジェクトを開始します。
 
プロジェクト名.cppというソースを書くウィンドウが作成されるので、ここにプログラムを書いていきます。
 
プロジェクト名.cppというソースを書くウィンドウが作成されるので、ここにプログラムを書いていきます。
 +
まず、ファイルから1行ずつ読み込んで結合するプログラムを作成します。
  
;ゲノムを読むサンプル
+
;ファイルから1行ずつ読み込むサンプル
 +
<pre>
 +
#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&#43;&#43;)
 +
      chr[size&#43;&#43;] = line[j];
 +
  }
 +
  chr[size] = 0;
 +
  cout << "Read " << size << " characters." << endl;
 +
}
 +
</pre>
 +
プログラムを書いたら、ビルド&rarr;"プロジェクト名"のビルド でコンパイルします。
 +
うまく動くことが確認できたら、書いた部分をクラスとして抽象化しましょう。
 +
他の染色体を読み込めるよう、
 +
* 染色体のサイズ (上のプログラムでは50MBの配列chr)
 +
* 実際の文字数 (size)
 +
をクラス化します。
 +
またゲノムの文字は大文字と小文字が混在しているので、全て小文字に変換する関数 toLowercase も用意します。
 +
 
 +
;21番染色体を読むサンプル
 
<pre>
 
<pre>
 
#include "stdafx.h"
 
#include "stdafx.h"
 
#include <iostream>
 
#include <iostream>
 
#include <fstream>
 
#include <fstream>
#include <time.h>
 
#include <cstdlib>
 
  
 
using namespace std;
 
using namespace std;
Line 40: Line 74:
 
int bufsize;
 
int bufsize;
 
int size;
 
int size;
 +
 
char toLowercase(char c)
 
char toLowercase(char c)
{
+
  {
switch (c) {
+
    switch (c) {
case 'A': return 'a';
+
    case 'A': return 'a';
case 'T': return 't';
+
    case 'T': return 't';
case 'G': return 'g';
+
    case 'G': return 'g';
case 'C': return 'c';
+
    case 'C': return 'c';
case 'N': return 'n';
+
    case 'N': return 'n';
case 'a': case 't': case 'g': case 'c': case 'n':
+
    case 'a': case 't': case 'g': case 'c': case 'n':
return c;
+
      return c;
default:
+
    default:
cout << "Unknown char " << c << endl;
+
      cout << "Unknown char " << c << endl;
exit(0);
+
      exit(0);
}
+
    }
}
+
  }
  
 
public:
 
public:
Chromosome(int siz) { chr = new char[siz]; bufsize = siz; }
+
  Chromosome(int siz) { chr = new char[siz]; bufsize = siz; }
~Chromosome() { delete[] chr; }
+
  ~Chromosome() { delete[] chr; }
Chromosome(const Chromosome&)
+
  Chromosome(const Chromosome&)
{ cout << "Error: copy constructor is called." << endl; exit(0); }
+
    { cout << "Error: copy constructor is called." << endl; exit(0); }
Chromosome& operator=(const Chromosome& x)
+
  Chromosome& operator=(const Chromosome& x)
{ cout << "Error: assignment constructor is called." << endl; exit(0); }
+
    { cout << "Error: assignment constructor is called." << endl; exit(0); }
  
void open(char* filename)
+
  void open(char* filename)
{
+
  {
char line[101];
+
    char line[100];
line[100] = 0; // 最後にヌル文字を入れておく
+
    size=0; // クラスのメンバー変数
size=0;
+
    ifstream fin(filename);
ifstream fin(filename);
+
    if (!fin)
if (!fin)
+
    { cout << "File not found" << endl; exit(0); }
{ cout << "File not found" << endl; exit(0); }
+
    fin.getline(line, 100);
// skip the first line
+
    while (!fin.getline(line, 100).eof()) {
fin.getline(line, 100);
+
      for(int j=0; line[j] != 0; j&#43;&#43;) {
while (!fin.getline(line, 100).eof())
+
        chr[size&#43;&#43;] = toLowercase(line[j]);
{
+
        if (size >= bufsize)
for(int j=0; line[j] != 0; j&#43;&#43;)
+
        { cout << "Error: bufsize overflow" << endl; exit(0); }
chr[size&#43;&#43;] = toLowercase(line[j]);
+
      }
if (size >= bufsize)
+
    }
{ cout << "Error: bufsize overflow" << endl; exit(0); }
+
    chr[size] = 0; // 最後にヌル文字
}
+
  }
chr[size] = 0;
+
}
+
  
void printGC()
+
  void printGC()
{
+
  {
int baseA =0, baseT =0, baseG = 0, baseC =0, baseN=0;
+
    int baseA =0, baseT =0, baseG = 0, baseC =0, baseN=0;
for (int i=0; i < size; i&#43;&#43;)
+
    for (int i=0; i < size; i&#43;&#43;)
switch (chr[i])
+
      switch (chr[i])
{  
+
      {  
case 'a': baseA&#43;&#43;; break;
+
        case 'a': baseA&#43;&#43;; break;
case 't': baseT&#43;&#43;; break;
+
        case 't': baseT&#43;&#43;; break;
case 'g': baseG&#43;&#43;; break;
+
        case 'g': baseG&#43;&#43;; break;
case 'c': baseC&#43;&#43;; break;
+
        case 'c': baseC&#43;&#43;; break;
case 'n': baseN&#43;&#43;; break;
+
        case 'n': baseN&#43;&#43;; break;
}
+
      }
cout << "GC content: " << (baseG+baseC)*100 / (baseA+baseT)  
+
      cout << "GC content: " << (baseG+baseC)*100 / (baseA+baseT)  
 
                     << "%" << endl;
 
                     << "%" << endl;
}
+
  }
 
};
 
};
  
Line 108: Line 141:
 
}
 
}
 
</pre>
 
</pre>
プログラムを書いたら、ビルド&rarr;"プロジェクト名"のビルド でコンパイルします。
 

Revision as of 16:18, 19 October 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というソースを書くウィンドウが作成されるので、ここにプログラムを書いていきます。 まず、ファイルから1行ずつ読み込んで結合するプログラムを作成します。

ファイルから1行ずつ読み込むサンプル
#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)

をクラス化します。 またゲノムの文字は大文字と小文字が混在しているので、全て小文字に変換する関数 toLowercase も用意します。

21番染色体を読むサンプル
#include "stdafx.h"
#include <iostream>
#include <fstream>

using namespace std;

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

char 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);
    }
  }

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); }
    fin.getline(line, 100);
    while (!fin.getline(line, 100).eof()) {
      for(int j=0; line[j] != 0; j++) {
        chr[size++] = toLowercase(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': baseA++; break;
        case 't': baseT++; break;
        case 'g': baseG++; break;
        case 'c': baseC++; break;
        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();
}
Personal tools
Namespaces

Variants
Actions
Navigation
metabolites
Toolbox