DNAseq on C++

実は院からは生物関係を専門にしてます。

各種DNAやらアミノ酸列とかをごにょごにょして、最終的には遺伝子操作による新人類を作り、愛の名の下に旧人類を抹殺する研究です。要するに、分子生物学に対し、情報学の立場からアプローチしようという分野です。バイオインフォマティクス(Bioinfomatics)とも言います。
別に生物に限った話ではなく、化学にも触手を伸ばしてます。医薬用外劇物の新物質の存在可能性をサポートベクターマシンで予想したりとか。こうなると、もはや旧人類に未来はありません。(旧人類の立場から見れば。詳しくはその辺のSFを読むとよいよ!)

さて、冗談は置いておくとして。

話を端的に持っていくと、まずDNAシーケンスをコンピュータ上で表現する必要があります。
DNAは塩基(base)という化学物質が連なったものです。A(アデニン)G(グアニン)C(シトシン)T(チミン)の4種類です。しかし、化学物質はそのままではメモリには入らないので、メモリに入れるような形にせにゃならんです。まず思いつくのはAGCTをそのままchar(or std::string)の配列にする方法です。しかし、これは実は不便。何が不便なのかと言うと、バイオインフォマティクスの分野では「塩基同士の演算」をよく使用するからです。結合したり、ランダム配列を作ったり、欠損したり…。
なので、これを数字にエンコードしたものを配列にします。(A,B,C,D)=(0,1,2,3)みたいに。連続したものにするとさらに便利。





余談だけど、つまりこれは2ビットの列だけあればDNAは表現できることになる。すると、実は人間の全DNAシーケンスはCD-R一枚にすべておさめることができる。その変のエロゲなんかよりは格段に少ないデータ量しかない。これには驚く。
CD1枚とか楽勝やろwww新人類始まってたwwwwとかと思う人もいるだろう。が、実はそうではない。これにも驚く。が、やはり残念ながら新人類の時代はもっとずっと先になりそうである。
というのも、人間の生命情報はDNAがすべてではないし、そもそもDNAすべてが遺伝情報ではないから。
後者から言うと、まずDNA内の遺伝情報は全体の2%程度しかない。他は、その遺伝情報を発現させるかどうかを外因的に決定する情報を持っている(と言われている)。残念ながら、そのほとんどの部分は解明されていない。前者について言えば、ミトコンドリアだって独自のDNAを持っているし、遺伝するRNAの存在も発見されている。こちらは、DNAよりもさらに解明がすすんで"いない"。
新人類の愛は、まだ遥か先に享受することになるだろう。





さて、あとは組むのみ!
オブジェクト思考すると、まず塩基(Base)が必要です。塩基がつながるとシーケンス(Sequence)になります。

その前に、塩基コードを定義しないといけません。

class Basecode {
private:
	std::vector<std::string> basenames;
public:
	Basecode();
	virtual ~Basecode();
	std::string name(const size_t code) const;
	int code(const std::string name) const;
};
Basecode::Basecode() {
	std::string b[] = {"A", "G", "C", "T"};
	for (int i = 0; i < 4; i++) {
		basenames.push_back(b[i]);
	}
}

std::string Basecode::name(const size_t code) const {
	return basenames.at(code);
}

int Basecode::code(const std::string name) const {
	std::vector<std::string>::const_iterator iter = basenames.begin();
	for (int i = 0; iter != basenames.end(); i++) {
		if (*iter == name) {
			return i;
		}
		iter++;
	}
	return -1;
}

これを使って、塩基クラスを作ります。

class Base {
private:
	int c; // 塩基コード
	std::string n; // 塩基名
public:
	Base();
	virtual ~Base();
	void makeBase(const std::string name);
	int code() const;
	std::string name() const;
};

void Base::makeBase(const std::string name) {
	c = Basecode().code(name);
	n = Basecode().name(c);
}
std::string Base::name() const {
	return n;
}
int Base::code() const {
	return c;
}

塩基と塩基コードはis-aでもいいとは思うけど、継承しないことにした。理由は、塩基コードクラスを継承して塩基クラスを作った場合、塩基にアクセスするたびに塩基コード列を舐めることになるから。実は最初継承で組んだんだけど、自然に組んだらその問題が発現することが分かってやめた。シーケンスを舐めるのにはO(n)かかるが、塩基コード列まで毎回舐めるとO(4n)になる。さすがにこれはちょっと…。よって、塩基クラスにコードと塩基名を持たせることにした。

配列のインデックスオーバーフローは例外に任せた。こういうのは言語に備わったものを使うのがシンプルだし安全かと。
塩基クラスに境界判定が無いのはそのせいです。塩基コードクラスが例外を投げてくれる。唯一、塩基コードから例外ではなく-1が帰る場合があるが、これは後でライブラリ内で使う時にはsize_tに入れる以外は特に使用しないので、その時にout of range例外を投げてくれるはず。(以降すべてoperator []とか使わずにatメソッドを使うので、これで安全。)

塩基コードのコンストラクタがカッコ悪いのが短所。改善したいなぁ。

で、シーケンス!

class Sequence {
private:
	std::vector<Base> seq;
public:
	Sequence();
	virtual ~Sequence();
	void readSequence(const std::string seq_str);
	size_t length() const;
	Base at(const size_t i) const;
};

void Sequence::readSequence(const std::string seq_str) {
	std::string::const_iterator iter = seq_str.begin();
	for (int i = 0; iter != seq_str.end(); i++) {
		Base b;
		b.makeBase(seq_str.substr(i, 1));
		seq.push_back(b);
		iter++;
	}
}

size_t Sequence::length() const {
	return this->seq.size();
}

Base Sequence::at(const size_t i) const {
	return seq.at(i);
}

これを使って、シーケンスアラインメントの課題を解いてみたが、ちゃんと動いてくれた。
あとはSequenceクラスを拡張するだけでどんどん便利な俺様ライブラリに!なって!