Rosalindを解く - DNA逆相補鎖変換

Rosalindの紹介はこちらから。

www.kimoton.com

本日は基礎中の基礎第3弾、DNA配列を逆相補鎖に変換するコードを書いていきたいと思います。

rosalind.info

生物学的知識のおさらい


参考:相補性 (分子生物学) - Wikipedia

核酸の一次構造は、核酸ポリマーの結合を構成する糖リン酸骨格に沿った核酸塩基(A, T, G, C)の順序によって決定されることを学びました。

しかし、一次構造のみからは分子に関する三次元構造について把握することはできません。

ジェームズ・ワトソンフランシス・クリックは、ロザリンド・フランクリンレイモンド・ゴスリングによって作成された高解像度X線画像を参考に、DNAの下記の構造を提案しました。

  1. DNA分子は、逆方向に走る2本の鎖で構成されている。
  2. 各塩基は反対側の鎖の塩基に結合する。アデニンは常にチミンと結合し、シトシンは常にグアニンと結合する。
  3. 2本の鎖は、二重らせんと呼ばれる長いらせん階段のような構造に撚り合わされている。

1, 2はDNAの二次構造を規定します。3は、DNA分子の三次構造を規定します。

ワトソンとクリックのモデルに倣って、2つの相補的な塩基の結合は塩基対(bp: base pair)と呼ばれます。
相補性により、片方の鎖の塩基の順序がわかれば、相補鎖の塩基配列を推測することができます。

問題

DNA配列は、「A」、「C」、「G」、および「U」を含むアルファベットから形成された配列です。DNA配列中では、「A」と「T」は互いに相補的であり、「C」と「G」も同様です。

DNAストリング $s$ の逆相補的配列は、$s$ の前後を逆にして、各塩基の相補的な塩基に変換することによって形成される配列 $s^c$ です(たとえば、「GTCA」の逆相補的配列は「TGAC」です)。

「長さが最大1000塩基長のDNA配列 $s$ が与えられる。 $s$ の逆相補鎖配列 $s^c$ を出力しなさい」

Given: A DNA string s of length at most 1000 bp.
Return: The reverse complement $s^c$ of $s$.

0. ディレクトリ構成

今回のディレクトリ構成は下記とします。

$ tree 
.
├── main.py # 実行ファイル
├── data 
│   └── rosalind_revc.txt # 入力データ
└── utils
    └── file_handler.py # 実行ファイルで使用するモジュール

1. 塩基配列ファイルを読むためのutil関数を使いまわす

以前のDNA塩基数カウントの際に使用したfile_handlerモジュールのload_seq関数を使い回します。

再記するとこんな関数。塩基配列のファイルを読んでSeqオブジェクトを返しています。

#!/usr/bin/env python

import os
from Bio.Seq import Seq

def load_seq(file_name: str, dir_name: str="data") -> Seq:
    """
    Arguments:
        file_name: Input file name
        dir_name: Input directory name

    Returns:
        Bio.Seq.Seq instance
    """
    # os.pathモジュールでディレクトリ名・ファイル名からファイルパスを作成
    file_path = os.path.join(dir_name, file_name)
    # with句で安全にファイルを開く
    with open(file_path, "r") as rf:
        # .strip()メソッドで余計な空白を除去
        seq_str = rf.read().strip()
    # Bio.Seq.Seqインスタンスを生成
    seq = Seq(seq_str)
    return seq

2. メイン関数を作成する

main.py ではload_seq()関数を使用して、data/rosalind_revc.txtファイルに対してDNA塩基配列の逆相補鎖変換を実行しています。
Seqオブジェクトでは、.reverse_complement()メソッドが用意されており、こちらを使用します。

今回は触れませんが、スクラッチで書く場合は相補鎖変換の辞書を用意してあげると良いでしょう。

from utils.file_handler import load_seq
from Bio.Seq import Seq

def transcribe_DNA_to_RNA(seq:Seq) -> str:
    """ Transcribe DNA to RNA

    Arguments: Bio.Seq.Seq instance
    Returns: str 
    """
    return seq.reverse_complement()

print(transcribe_DNA_to_RNA(load_seq("rosalind_revc.txt")))

実行

$ cat data/rosalind_revc.txt
AAAACCCGGT

$ python main.py
## ACCGGGTTTT

コード例

下記で公開しているので参考にしてみて下さい。

repl.it