Rosalindを解く - DNA塩基数のカウント

Rosalindの紹介はこちらから。

本日は基礎中の基礎、DNA塩基数を数えるコードを書いていきたいと思います。

rosalind.info

生物学的知識のおさらい

参考：クロマチン - Wikipedia

すべての生物を構成する細胞は、生命の構成要素であると考えられています。
我々真核細胞の構成要素である核は、光学顕微鏡で見るとクロマチンと呼ばれる高分子から構成されていることが分かります。 このクロマチンのほとんどは、有糸分裂（真核細胞の分裂）の間、染色体と呼ばれる長くて細い紐に凝縮します。

クロマチンに含まれる高分子の1つは核酸と呼ばれます。核酸を構成する分子はヌクレオチドと呼ばれ、この分子が連なることでクロマチンが構成されています。核酸はデオキシリボースと呼ばれる糖分子、リン酸と呼ばれる負に帯電したイオン、および核酸塩基（略して「塩基」）と呼ばれる化合物から構成されます。

特定の種類の核酸のヌクレオチドには常に同じ糖分子とリン酸分子が含まれており、塩基のみが異なっています。核酸塩基は、アデニン（A）、シトシン（C）、グアニン（G）、およびチミン（T）から構成されています。

DNAは、バクテリアを含む地球上のすべての生物に見られます。多くのウイルス（生きていないものと見なされることが多い）にも見られます。因みにDNAとゲノムは概念的には ゲノム >> DNAであり、ゲノム は生物の染色体に含まれるDNAの全量を指す用語として使用されています。

問題

「最大1000塩基長のDNA文字列の各塩基数をカウントしなさい」

Given: A DNA string s of length at most 1000 nt.
Return: Four integers (separated by spaces) counting the respective number of times that the symbols 'A', 'C', 'G', and 'T' occur in s.

解

0. ディレクトリ構成

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

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

1. 塩基配列ファイルを読むためのutil関数を作る

util/file_handler.pyとして、塩基配列ファイルを読むためのモジュールを作成します。
ポイントとしては、Biopythonパッケージより、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_dna.txtファイルに対して塩基数のカウントを実行しています。

#!/usr/bin/env python

"""python
A solution to a ROSALIND bioinformatics problem.

Problem Title: Counting DNA Nucleotides: DNA
URL: http://rosalind.info/problems/dna/
"""

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

def calc_ACGT_count(seq:Seq) -> str:
    """ Concatenate nucleotide counts in the order of  A, C, G, T splitted by space (' ')

    Arguments: Bio.Seq.Seq instance
    Returns: str 
    """
    # .count()メソッドで各塩基の数をカウントし、A, C, G, Tの順番のリストを作成
    ACGT_count_list = [seq.count("A"),
                       seq.count("C"),
                       seq.count("G"),
                       seq.count("T")]
    # .join()メソッドにより、スペース（" "）区切りで結合
    return " ".join([str(i) for i in ACGT_count_list])

# 標準出力に出力
print(calc_ACGT_count(load_seq("rosalind_dna.txt")))

実行

$ cat data/rosalind_dna.txt
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC

$ python main.py
## 20 12 17 21

コード例

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

repl.it