Aritalab:Lecture/Programming/Cpp/Genome
From Metabolomics.JP
< Aritalab:Lecture | Programming | Cpp
ゲノムデータを解析する
ゲノムデータを解析するプログラムの作成法を、ステップ毎に紹介します。
データのダウンロード
まず 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(); }
プログラムを書いたら、ビルド→"プロジェクト名"のビルド でコンパイルします。