Aritalab:Lecture/Programming/Cpp/Genome
From Metabolomics.JP
< Aritalab:Lecture | Programming | Cpp(Difference between revisions)
(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) | + | まず [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によるプログラム作成=== | ||
− | + | 長さ100文字のバッファを使って行を読み、50MBのヒープ領域に21番染色体のデータを格納しましょう。 | |
+ | Visual Studioインストールしたら、File→新規作成→プロジェクトを選択して、新しいプロジェクトを開始します。 | ||
プロジェクト名.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++) | ||
+ | chr[size++] = line[j]; | ||
+ | } | ||
+ | chr[size] = 0; | ||
+ | cout << "Read " << size << " characters." << endl; | ||
+ | } | ||
+ | </pre> | ||
+ | プログラムを書いたら、ビルド→"プロジェクト名"のビルド でコンパイルします。 | ||
+ | うまく動くことが確認できたら、書いた部分をクラスとして抽象化しましょう。 | ||
+ | 他の染色体を読み込めるよう、 | ||
+ | * 染色体のサイズ (上のプログラムでは50MBの配列chr) | ||
+ | * 実際の文字数 (size) | ||
+ | をクラス化します。 | ||
+ | またゲノムの文字は大文字と小文字が混在しているので、全て小文字に変換する関数 toLowercase も用意します。 | ||
+ | |||
+ | ;21番染色体を読むサンプル | ||
<pre> | <pre> | ||
#include "stdafx.h" | #include "stdafx.h" | ||
#include <iostream> | #include <iostream> | ||
#include <fstream> | #include <fstream> | ||
− | |||
− | |||
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) { | |
− | + | 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: | 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; | << "%" << endl; | ||
− | + | } | |
}; | }; | ||
Line 108: | Line 141: | ||
} | } | ||
</pre> | </pre> | ||
− |
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(); }