Aritalab:Lecture/Programming/Cpp/Genome

ゲノムデータを解析する

ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。

データのダウンロード

まず UCSD Downloadサイトに行って、ヒトゲノムのデータをダウンロードしてみます。ここではData set by chromosomeのセクションに進んで日本がシーケンスした21版染色体 (chr21.fa.gz) を取得します。faというのはFastA形式の意味です。96万行もあるファイルで、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によるプログラム作成

インストールしたら、File→新規作成→プロジェクトを選択して、新しいプロジェクトを開始します。 プロジェクト名.cppというソースを書くウィンドウが作成されるので、ここにプログラムを書いていきます。

；ゲノムを読むサンプル

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

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[101];
		line[100] = 0; // 最後にヌル文字を入れておく
		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++] = 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();
}

プログラムを書いたら、ビルド→"プロジェクト名"のビルド でコンパイルします。