GC含有量を大量に計算する — `for`ループを捨てるべき理由

このトピックを終えると

教科書で学んだベクトル化とBig-Oを組み合わせて、数千～数百万個のDNA配列のGC含量をforループを使わずに瞬時に計算するツールを作成できます。また、「なぜdf * 2は0.1秒なのに、自分のforループは10分もかかるのか？」という長年の疑問に、ハードウェアレベルで答えることができるようになります。

この記事は教育目的の一般的な例です。GC含量は、プライマー設計、シーケンシング品質管理、種分類など、分子生物学のあらゆる場面で登場する基本的な指標であるため、題材として選びました。

GC含量とは何か？なぜ毎回計算するのか

DNAはA、T、G、Cの4つの文字で構成されています。このうち、GとCの割合がGC含量です。大したことないように見えますが、実験では非常に重要です。

GとCは三重水素結合（AとTは二重）であるため、GC含量が高いほど二重鎖がより強く結合します。→ プライマーの**融解温度（Tm）**に影響を与えます。
GC含量が極端であると、PCRがうまくいきません。そのため、プライマーは通常、**40～60％**のGC含量を目指します。
シーケンシングデータにおいて、リードのGC分布が異常な場合、汚染やバイアスの兆候です。

そのため、「配列1つのGC含量」は次のように簡単に計算します。

python

def gc_content_one(seq):
    gc = seq.count("G") + seq.count("C")
    return gc / len(seq)

assert abs(gc_content_one("GGCC") - 1.0) < 1e-9
assert abs(gc_content_one("ATAT") - 0.0) < 1e-9
assert abs(gc_content_one("ATGC") - 0.5) < 1e-9

問題は、配列が1つではないことです。シーケンシング1回で、何百万ものリードが得られます。これをどのように処理するかが、初心者と熟練者の違いです。

まずは完成品を見てみましょう（ブラックボックスを先に実行する）

私たちが作るツールは、数百万個の配列のGC含量を一度に受け取り、統計と分布を抽出します。

python

result = gc_report(reads)      # reads = 100万個の配列

print(result.summary())

text

=== GC含量レポート (1,000,000配列) ===
平均GC含量 : 0.502
標準偏差 : 0.071
40〜60％の範囲 : 92.4％
40％未満 : 3.8％（AT-rich）
60％超 : 3.8％（GC-rich）
処理時間 : 0.09秒 ← forループだった場合、数十秒〜分

重要なのは、この0.09秒です。これから、まずforループ版を作成し、なぜ遅いのかを確認し、それをベクトル化してこの速度に到達するまでの過程を説明します。

このツールはどのような部品で構成されているか（部品分解図）

text

GC含有量大量計算ツール
   ┌────────────────────────────────────────────────────┐
   │  [基準線] for ループバージョン ── 部品：反復処理   │  ← 完成品として提供（意図的に遅い対照群）
   │              │                                     │
   │              ▼                                     │
   │  [核心] シーケンス → 数値行列 ── 部品：ベクトル化  │  ← 自分で作成 ★
   │              │                                     │
   │              ▼                                     │
   │  [理解] なぜ高速なのか ──────── 部品：Big-O        │  ← 自分で作成 ★
   │              │            + メモリ構造             │
   │              ▼                                     │
   │  [出力] 統計・分布レポート                         │
   └────────────────────────────────────────────────────┘

部品	どこで学んだか	このツールで行うこと
反復処理	`loops-for-while-break`	遅い基準線（対照群）を作成
ベクトル化	`vectorization-why-no-for-loop`	シーケンスを行列に変換して一括計算
Big-O	`big-o-notation`	2つの方式の速度差を説明

📌 これらの概念を初めて見る場合（上部のリンク）

ベクトル化 — forループを捨てるべき理由

Big-O記法

連続メモリと配列インデックス — ベクトル化がなぜ高速なのかの根本

今回の内容は、新たに自分で作成する概念がちょうど2つ(ベクトル化、Big-O)なので、マイクロレベルです。代わりに、「なぜ高速なのか」を深く掘り下げます。

ステップ1：まずは単純なforループから（反復処理、完成版を提供）

まず、単純なforループを使用したバージョンを作成します。これは学習の目的というよりも、比較対象として完成版を提供します。loops-for-while-breakで学んだことをそのまま適用します。

python

def gc_content_loop(sequences):
    """各配列に対してforループでGC含量を計算。O(N × L)だが、Pythonのループであるためオーバーヘッドが大きい。"""
    result = []
    for seq in sequences:
        gc = 0
        for base in seq:            # 各文字に対してPythonがループ処理を行う
            if base == "G" or base == "C":
                gc += 1
        result.append(gc / len(seq))
    return result

reads = ["ATGCGC", "AAATTT", "GCGCGC", "ATATGC"]
loop_out = gc_content_loop(reads)

assert abs(loop_out[0] - 4/6) < 1e-9
assert abs(loop_out[1] - 0.0) < 1e-9
assert abs(loop_out[2] - 1.0) < 1e-9

このコードは正しく動作します。結果も正確です。問題は速度だけです。配列の数がN個、各配列の長さがLの場合、PythonインタープリターはN × L回ループを処理します。Pythonのループは、各ループで「この変数の型は何か？」「この演算をどのように行うか？」を毎回確認するため、処理に時間がかかります。配列が100万個、各配列の長さが150の場合、1億5千万回のPythonループ処理が発生します。ここでコーヒーが必要になります。

ステップ2：配列を数値行列に変換 ★（ベクトル化）

✍️ 手動で実装する部分。 部品 = ベクトル化。目的：Pythonのループをなくし、NumPyに「すべてを一度に計算して」と指示する。

ベクトル化の基本的な考え方はこれです。文字を1つずつ数えるのではなく、配列全体を数値配列に変換し、配列演算を1回で処理する。

同じ長さの配列（シーケンスリードは通常同じ長さです）を縦に積み重ねると、N × L サイズの2次元文字行列になります。これをバイト（数値）に変換すると、NumPyでまとめて処理できるようになります。

🔎 NumPyが高速である理由（補足 — contiguous-memory-array-indexing） Pythonリストは、値がメモリのあちこちに分散した「アドレス帳」です。NumPy配列は、値がメモリに1列にぎっしりと並んでいます（連続メモリ）。そのため、CPUはSIMD（1つの命令で複数のデータを同時に処理）の高速化を使用できます。forループが「1人ずつ列に並べて計算」であるのに対し、ベクトル化は「1つのクラス全体を同時に採点」するようなものです。

python

import numpy as np

def to_matrix(sequences):
    """同じ長さの配列を(N, L) uint8行列に変換。"""
    # 各配列をバイトに変換して積み重ねる。"ATGC" → [65, 84, 71, 67]
    return np.array([np.frombuffer(s.encode("ascii"), dtype=np.uint8) for s in sequences])

mat = to_matrix(["ATGC", "GGCC"])
assert mat.shape == (2, 4)
assert mat[0, 0] == ord("A")     # 65
assert mat[1, 0] == ord("G")     # 71

次に、この行列から「GまたはCであるセル」を一度に見つけます。ここがベクトル化の核心です。

python

G, C = ord("G"), ord("C")

def gc_content_vectorized(sequences):
    mat = to_matrix(sequences)
    # (mat == G) | (mat == C) → 同じサイズのTrue/False行列。Pythonのループなし！
    is_gc = (mat == G) | (mat == C)
    # 行（=配列）ごとにTrueの数を数え、長さで割る
    return is_gc.sum(axis=1) / mat.shape[1]

vec_out = gc_content_vectorized(["ATGCGC", "AAATTT", "GCGCGC", "ATATGC"])

assert abs(vec_out[0] - 4/6) < 1e-9
assert abs(vec_out[1] - 0.0) < 1e-9
assert abs(vec_out[2] - 1.0) < 1e-9
# 何よりも：forループバージョンと結果が完全に同じでなければならない
loop_out = gc_content_loop(["ATGCGC", "AAATTT", "GCGCGC", "ATATGC"])
assert np.allclose(vec_out, loop_out)

最後の assert np.allclose(...) が重要です。結果は同じで、速度だけが異なる — これがベクトル化の約束です。(mat == G) | (mat == C) この1行にはforループはありません。NumPyは内部的にC言語レベルで、行列全体を処理します。私たちが目にするのは1行ですが、CPUにとってはSIMDによる一回の処理です。

🤔 自己説明プロンプト is_gc.sum(axis=1)のaxis=1は、「行方向に加算する」という意味です。もしaxis=0に変更したら、何を数えることになるでしょうか？（ヒント：行=配列、列=位置。列方向に加算すると、「各位置でGCである配列が何個か」となり、これはGC含量ではなく、位置ごとのGC偏りになります。）

ステップ3：どれくらい高速化されたかを計測する ★ (Big-Oを手計算で)

✍️ 自分で埋める部分。 要素 = Big-O。目標：2つの方法の速度差を数値で確認し、その理由を説明する。

python

import random, time

random.seed(0)
# 長さ100の配列を5万個生成
reads = ["".join(random.choice("ATGC") for _ in range(100)) for _ in range(50000)]

t0 = time.time(); out_loop = gc_content_loop(reads);        loop_time = time.time() - t0
t0 = time.time(); out_vec  = gc_content_vectorized(reads);  vec_time  = time.time() - t0

print(f"forループバージョン：{loop_time:.3f}秒")
print(f"ベクトル化バージョン：{vec_time:.3f}秒")
print(f"倍速：{loop_time / vec_time:.0f}倍")

# 結果は同じで、ベクトル化の方がはるかに高速であるべき
assert np.allclose(out_loop, out_vec)
assert vec_time < loop_time

どちらの方法も処理量はO(N × L)で同じです。Big-Oだけを見ると同等です。しかし、なぜ速度が数十倍も違うのでしょうか？

ここがBig-Oの落とし穴であり、魅力でもあります。Big-Oは「増加率」を言い、定数については言いません。 forループとベクトル化はどちらもO(N×L)ですが、1回の演算にかかる定数コストは天と地ほど違います。

forループ：1文字ごとにPythonインタープリタが型チェック、オブジェクト生成、分岐を行います。定数が大きいです。
ベクトル化：連続メモリとSIMDを使用して、一度に数十個を処理します。定数が小さいです。

そのため、実務では「Big-Oが同じでもベクトル化する」が格言です。アルゴリズムの複雑さ（増加率）を減らすことと、定数を減らすこと（ベクトル化）はどちらもパフォーマンスの軸です。

🤔 自己説明プロンプト 配列の数を5万から50万に10倍に増やすと、forループの時間とベクトル化の時間はそれぞれおよそ何倍になるでしょうか？（ヒント：どちらもO(N)なので、理論上は10倍。しかし、定数の違いにより、絶対時間は依然としてベクトル化の方が圧倒的に高速です。）

部品を一つに — 最終レポート

python

def gc_report(sequences):
    gc = gc_content_vectorized(sequences)
    in_range = np.mean((gc >= 0.4) & (gc <= 0.6))
    return {
        "n": len(sequences),
        "mean": float(gc.mean()),
        "std": float(gc.std()),
        "in_range_40_60": float(in_range),
        "at_rich": float(np.mean(gc < 0.4)),
        "gc_rich": float(np.mean(gc > 0.6)),
    }

report = gc_report(reads)
assert report["n"] == 50000
assert 0.45 < report["mean"] < 0.55  # ランダムな ATGC なので、平均は 0.5 に近い
assert abs(report["in_range_40_60"] + report["at_rich"] + report["gc_rich"] - 1.0) < 1e-9

最後の assert は良い習慣です。3つの区間（範囲内 / AT-リッチ / GC-リッチ）の割合をすべて足すと、正確に 1 になるはずです。そうでない場合、どこかで境界条件に問題があることを意味します。このように、「合計が 1」のように必ず成立しなければならない不変条件を assert で記述しておくと、後でコードを修正するときに、すぐに間違いを発見できます。

他の道もある（マルチパス考察）

Biopythonのgc_fraction: 1つの配列のGC含量を安全に取得できます（小文字・Nの処理を含む）。ただし、リストを渡すと内部でPythonループを使用するため、大量のデータでは速度が遅くなります。こちらの方法が有利な場合: 長さが同じリードが数百万個ある場合。
Pandasのstrアクセサ: pd.Series(reads).str.count("G")でも可能です。便利ですが、文字列演算であるため、NumPyのバイト行列よりも速度が遅くなります。
長さの異なる配列: to_matrix関数は、配列の長さが同じである必要があります。長さがそれぞれ異なる場合はどうすればよいでしょうか？短い方をパディングするか、配列ごとにstr.countをNumPyでラップするハイブリッドを使用します。トレードオフ: ベクトル化による純粋な速度 ↔ 柔軟性。

核心: ベクトル化は、「長さが同じ大量のデータ」において最も強力です。データの形状を見て、最適な方法を選択することが真のスキルです。

次のステップへ（下部の外部リンク）

ベクトル化がなぜ高速なのかをより深く理解するには → 連続したメモリと配列インデックス
複数の配列をテーブルにまとめ、条件別の統計を計算したい場合は → ベクトル化の後にgroupbyを組み合わせた応用編
まず重複を取り除く必要がある場合は → 応用編プライマー重複除去エンジン

実際に試してみましょう（独立した課題）

AT含量の同時計算: GCとともに、AT含量も一度の行列演算で求めましょう。((mat==A)|(mat==T)) GC+ATが1になることをassertで確認します。
位置ごとのGC偏り: axis=0で合計して「リードの各位置でのGC比率」を求めます。シーケンスの先頭部分が偏っているかどうかを確認できます。
品質フィルター: GCが20〜80%外のリードのインデックスをベクトル演算で抽出します。(np.where)
挑戦: 長さが異なる配列が混ざっている場合、forループなしでGCを計算する方法を1つ設計してください。（ヒント：総GC数と総長さを別々にベクトル化する）

まとめ

「GC含量の計算」という非常に基本的なタスクを通して、同じ正解を数十倍早く得る方法を学びました。

繰り返し処理は正確ですが、Pythonループの定数コストのために遅くなります。
ベクトル化は、一連のデータを連続したメモリ行列にすることで、SIMDによって一度に処理します。
Big-Oは、「増加率は同じでも定数が異なる」ことを教えてくれ、なぜ高速化されたのかを説明してくれます。

forループに慣れているからといって、それが常に最善の選択肢とは限りません。データが大きくなると、**「これを配列演算一つで置き換えることはできないか？」**とまず考える習慣を身につけること——それがデータを取り扱う人のセンスです。

この記事は一般的な教育例です。実際のパイプラインでは、品質スコア、アダプター除去、複数のファイルなど、より多くの変数が存在します。その詳細なバージョンは、このベクトル化のフレームワークの上に構築することで実現できます。

BioPlayground

GC含有量を大量に計算するためのツール — forループの使用を避けるべき理由

GC含有量を大量に計算する — `for`ループを捨てるべき理由

このトピックを終えると

GC含量とは何か？なぜ毎回計算するのか

まずは完成品を見てみましょう（ブラックボックスを先に実行する）

このツールはどのような部品で構成されているか（部品分解図）

ステップ1：まずは単純なforループから（反復処理、完成版を提供）

ステップ2：配列を数値行列に変換 ★（ベクトル化）

ステップ3：どれくらい高速化されたかを計測する ★ (Big-Oを手計算で)

部品を一つに — 最終レポート

他の道もある（マルチパス考察）

次のステップへ（下部の外部リンク）

実際に試してみましょう（独立した課題）

まとめ

💬 質問・コメント

GC含有量を大量に計算するためのツール — forループの使用を避けるべき理由

GC含有量を大量に計算する — forループを捨てるべき理由

このトピックを終えると

GC含量とは何か？なぜ毎回計算するのか

まずは完成品を見てみましょう（ブラックボックスを先に実行する）

このツールはどのような部品で構成されているか（部品分解図）

ステップ1：まずは単純なforループから（反復処理、完成版を提供）

ステップ2：配列を数値行列に変換 ★（ベクトル化）

ステップ3：どれくらい高速化されたかを計測する ★ (Big-Oを手計算で)

部品を一つに — 最終レポート

他の道もある（マルチパス考察）

次のステップへ（下部の外部リンク）

実際に試してみましょう（独立した課題）

まとめ

💬 質問・コメント

GC含有量を大量に計算する — `for`ループを捨てるべき理由