GC 함량 대량 계산기 — for 루프를 버려야 하는 이유

이 토픽을 마치면

교과서에서 배운 벡터화와 Big-O를 엮어서, 수만~수백만 개 DNA 서열의 GC 함량을 for 루프 없이 순식간에 계산하는 도구를 만들 수 있습니다. 그리고 "왜 df * 2는 0.1초인데 내 for 루프는 10분인가?"라는 오래된 의문에 하드웨어 수준에서 답할 수 있게 됩니다.

이 글은 교육용 일반 예제입니다. GC 함량은 프라이머 설계, 시퀀싱 품질 관리, 종 분류 등 분자생물학 어디에나 나오는 기본 지표라서 소재로 삼았습니다.

GC 함량이 뭐길래 매번 계산할까

DNA는 A·T·G·C 네 글자로 되어 있습니다. 이 중 G와 C의 비율이 GC 함량입니다. 별거 아닌 것 같지만 실험에서 엄청나게 중요합니다.

• G와 C는 삼중 수소결합(A·T는 이중)이라, GC가 높으면 이중 가닥이 더 단단히 붙습니다. → 프라이머의 녹는점(Tm) 을 좌우합니다.
• GC 함량이 극단적이면 PCR이 잘 안 됩니다. 그래서 프라이머는 보통 40~60% GC를 노립니다.
• 시퀀싱 데이터에서 리드들의 GC 분포가 이상하면 오염이나 편향의 신호입니다.

그래서 "서열 하나의 GC 함량"은 이렇게 간단히 구합니다.

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

문제는 서열이 하나가 아니라는 데 있습니다. 시퀀싱 한 번이면 리드가 수백만 개입니다. 이걸 어떻게 처리하느냐에서 초보와 고수가 갈립니다.

먼저 완성품을 봅시다 (블랙박스 먼저 돌리기)

우리가 만들 도구는 수백만 개 서열의 GC 함량을 한 번에 받아, 통계와 분포까지 뽑아 줍니다.

result = gc_report(reads)      # reads = 100만 개 서열

print(result.summary())

=== 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 루프 버전을 먼저 만들어 왜 느린지 보고, 그걸 벡터화로 갈아엎어 이 속도에 도달하는 과정을 밟습니다.

이 도구는 어떤 부품으로 조립되는가 (부품 분해도)

GC 함량 대량 계산기
   ┌────────────────────────────────────────┐
   │  [기준선] for 루프 버전 ── 부품: 반복문   │  ← 완성 제공 (일부러 느린 대조군)
   │              │                           │
   │              ▼                           │
   │  [핵심] 서열 → 숫자 행렬 ── 부품: 벡터화  │  ← 직접 만든다 ★
   │              │                           │
   │              ▼                           │
   │  [이해] 왜 빠른가 ──────── 부품: Big-O   │  ← 직접 만든다 ★
   │              │            + 메모리 구조   │
   │              ▼                           │
   │  [출력] 통계·분포 리포트                  │
   └────────────────────────────────────────┘

📌 이 개념들 처음 본다면 (상단 진입 링크)

• 벡터화 — for 루프를 버려야 하는 이유
• Big-O 표기법
• 연속 메모리와 배열 인덱싱 — 벡터화가 왜 빠른지의 뿌리

이번 편은 새로 직접 만드는 개념이 딱 2개(벡터화, Big-O)라 마이크로급입니다. 대신 "왜 빠른가"를 아주 깊이 팝니다.

만들기 1단계 — 느린 기준선부터 (반복문, 완성 제공)

먼저 순진한 for 루프 버전을 만듭니다. 이건 학습 목표가 아니라 대조군이라 완성해서 드립니다. loops-for-while-break에서 배운 그대로예요.

def gc_content_loop(sequences):
    """서열마다 for 루프로 GC 함량을 계산. O(N × L)이지만 파이썬 루프라 상수가 크다."""
    result = []
    for seq in sequences:
        gc = 0
        for base in seq:            # 글자 하나하나 파이썬이 돈다
            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이면 파이썬 인터프리터가 N × L번 루프를 돕니다. 파이썬 루프는 한 바퀴마다 "이 변수 타입이 뭐지? 이 연산 어떻게 하지?"를 매번 확인하느라 느립니다. 리드 100만 개 × 길이 150이면 1억 5천만 번의 파이썬 루프. 여기서 커피가 필요해집니다.

만들기 2단계 — 서열을 숫자 행렬로 ★ (벡터화)

✍️ 직접 채우는 구간. 부품 = 벡터화. 목표: 파이썬 루프를 없애고, numpy에게 "전부 한 번에 계산해"라고 시킨다.

벡터화의 핵심 발상은 이겁니다. 글자를 하나씩 세지 말고, 서열 전체를 숫자 배열로 만든 뒤 배열 연산 한 방으로 처리한다.

길이가 같은 서열들(시퀀싱 리드는 보통 길이가 같습니다)을 세로로 쌓으면 N × L 크기의 2차원 문자 행렬이 됩니다. 이걸 바이트(숫자)로 바꾸면 numpy가 통째로 다룰 수 있습니다.

🔎 왜 numpy가 빠른가? (드로어 — contiguous-memory-array-indexing) 파이썬 리스트는 값들이 메모리 여기저기 흩어진 "주소록"입니다. numpy 배열은 값들이 메모리에 한 줄로 빽빽이 붙어 있습니다(연속 메모리). 그래서 CPU가 SIMD(한 명령으로 여러 데이터 동시 처리) 가속을 쓸 수 있습니다. for 루프가 "한 명씩 줄 세워 계산"이라면, 벡터화는 "한 반 전체를 동시에 채점"하는 겁니다.

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인 칸"을 한 번에 찾습니다. 여기가 벡터화의 심장입니다.

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

def gc_content_vectorized(sequences):
    mat = to_matrix(sequences)
    # (mat == G) | (mat == C) → 같은 크기의 True/False 행렬. 파이썬 루프 없음!
    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) 이 한 줄에는 for가 없습니다. numpy가 내부적으로 C 언어 수준에서 전체 행렬을 훑죠. 우리 눈엔 한 줄, CPU에겐 SIMD 한 방입니다.

🤔 자기설명 프롬프트 is_gc.sum(axis=1)에서 axis=1은 "행 방향으로 더하라"는 뜻입니다. 만약 axis=0으로 바꾸면 무엇을 세게 될까요? (힌트: 행=서열, 열=위치. 열 방향으로 더하면 "각 위치에서 GC인 서열이 몇 개인가"가 됩니다 — 이건 GC 함량이 아니라 위치별 GC 편향이죠.)

만들기 3단계 — 얼마나 빨라졌는지 재기 ★ (Big-O를 손으로)

✍️ 직접 채우는 구간. 부품 = Big-O. 목표: 두 방식의 속도 격차를 숫자로 확인하고 왜 그런지 설명한다.

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}x")

# 결과는 같고, 벡터화가 훨씬 빨라야 한다
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)이지만, 한 번의 연산에 드는 상수 비용이 하늘과 땅입니다.

• for 루프: 한 글자마다 파이썬 인터프리터가 타입 확인·객체 생성·분기. 상수가 큽니다.
• 벡터화: 연속 메모리 + SIMD로 한 번에 수십 개씩. 상수가 작습니다.

그래서 실무에서는 "Big-O가 같아도 벡터화하라"가 격언입니다. 알고리즘 복잡도(증가율)를 줄이는 것과, 상수를 줄이는 것(벡터화)은 둘 다 성능의 축입니다.

🤔 자기설명 프롬프트 서열 개수를 5만에서 50만으로 10배 늘리면, for 루프 시간과 벡터화 시간은 각각 대략 몇 배가 될까요? (힌트: 둘 다 O(N)이라 이론상 10배. 하지만 상수 차이 때문에 절대 시간은 여전히 벡터화가 압도합니다.)

부품을 하나로 — 완성된 리포트

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가 좋은 습관입니다. 세 구간(범위 내 / AT-rich / GC-rich)의 비율을 다 더하면 정확히 1이어야 합니다. 안 그러면 어딘가 경계 조건이 새고 있다는 뜻이죠. 이렇게 "합이 1"처럼 반드시 성립해야 하는 불변식을 assert로 박아 두면, 나중에 코드를 고쳐도 실수를 즉시 잡습니다.

다른 길도 있다 (멀티패스 성찰)

• Biopython gc_fraction: 서열 하나의 GC를 안전하게 구해 줍니다(소문자·N 처리 포함). 다만 리스트를 넘기면 내부적으로 파이썬 루프라 대량엔 느립니다. 우리 방식이 나은 때: 길이가 같은 리드가 수백만 개일 때.
• Pandas str 접근자: pd.Series(reads).str.count("G") 도 됩니다. 편하지만 문자열 연산이라 numpy 바이트 행렬보다 느립니다.
• 길이가 다른 서열: 우리 to_matrix는 길이가 같아야 합니다. 길이가 제각각이면? 짧은 쪽을 패딩하거나, 서열별 str.count를 numpy로 감싸는 하이브리드를 씁니다. 트레이드오프: 벡터화의 순수한 속도 ↔ 유연성.

핵심: 벡터화는 "길이가 같은 대량 데이터"에서 최강입니다. 데이터 모양을 먼저 보고 무기를 고르는 게 진짜 실력입니다.

다음 단계로 (하단 출구 링크)

• 벡터화가 왜 빠른지 더 깊이 → 연속 메모리와 배열 인덱싱
• 여러 서열을 표로 묶어 조건별 통계를 내고 싶다면 → 벡터화 다음에 groupby를 엮은 응용편으로
• 중복부터 걸러야 한다면 → 응용편 프라이머 중복 제거기

직접 해보기 (독립 문제)

1. AT 함량 동시 계산: GC와 함께 AT 함량도 한 번의 행렬 연산으로 구하세요. ((mat==A)|(mat==T)) GC+AT가 1이 되는지 assert로 확인.
2. 위치별 GC 편향: axis=0으로 더해서 "리드의 각 위치에서 GC 비율"을 구하세요. 시퀀싱 앞부분이 편향됐는지 볼 수 있습니다.
3. 품질 필터: GC가 20~80% 밖인 리드의 인덱스를 벡터 연산으로 뽑으세요. (np.where)
4. 도전: 길이가 다른 서열이 섞여 있을 때, for 루프 없이 GC를 구하는 방법을 하나 설계해 보세요. (힌트: 총 GC 개수와 총 길이를 따로 벡터화)

정리

우리는 "GC 함량 계산"이라는 아주 기본적인 작업에서, 같은 정답을 수십 배 빠르게 얻는 법을 배웠습니다.

• 반복문은 정확하지만 파이썬 루프의 상수 비용 때문에 느렸습니다.
• 벡터화는 서열을 연속 메모리 행렬로 만들어, SIMD로 한 번에 처리했습니다.
• Big-O는 "증가율은 같아도 상수가 다르다"를 가르쳐, 왜 빨라졌는지 설명해 줬습니다.

for 루프가 손에 익었다고 그게 항상 답은 아닙니다. 데이터가 커지는 순간, "이걸 배열 연산 하나로 바꿀 수 있나?" 를 먼저 묻는 습관 — 그게 데이터를 다루는 사람의 감각입니다.

이 글은 일반 교육 예제입니다. 실제 파이프라인에서는 품질 점수, 어댑터 제거, 멀티 파일 등 변수가 훨씬 많습니다. 그 상세한 버전은 이 벡터화 뼈대 위에 여러분이 붙이면 됩니다.