Scientific Infectious Disease

2. AMR 遺伝子検出

def detect_amr_genes(assembly_fasta, database="resfinder"):
    """
    薬剤耐性（AMR）遺伝子の検出。

    データベース:
      - ResFinder: 後天性耐性遺伝子
      - CARD (RGI): 包括的 AMR データベース
      - AMRFinderPlus: NCBI 統合 AMR 検出

    結果カテゴリ:
      - 耐性遺伝子（acquired resistance genes）
      - 点変異（point mutations）
      - 耐性表現型予測
    """
    import subprocess
    import json

    if database == "resfinder":
        cmd = (f"python -m resfinder -ifa {assembly_fasta} "
               f"--acquired --point -o resfinder_results/")
        subprocess.run(cmd, shell=True, check=True)

        with open("resfinder_results/ResFinder_results_tab.txt") as f:
            lines = f.readlines()
        results = parse_resfinder_output(lines)

    elif database == "card":
        cmd = f"rgi main -i {assembly_fasta} -o rgi_results -t contig -a BLAST"
        subprocess.run(cmd, shell=True, check=True)
        results = pd.read_csv("rgi_results.txt", sep="\t")

    n_genes = len(results) if isinstance(results, list) else len(results)
    print(f"  AMR: {n_genes} resistance genes detected ({database})")
    return results


def parse_resfinder_output(lines):
    """ResFinder 出力をパースする。"""
    results = []
    for line in lines[1:]:
        fields = line.strip().split("\t")
        if len(fields) >= 6:
            results.append({
                "gene": fields[0],
                "identity": float(fields[1]),
                "coverage": float(fields[2]),
                "phenotype": fields[5] if len(fields) > 5 else "Unknown",
            })
    return results

3. 分子疫学タイピング

def molecular_typing(assembly_fasta, organism, scheme="mlst"):
    """
    分子疫学タイピング。

    scheme:
      - "mlst": Multi-Locus Sequence Typing（7 遺伝子座）
      - "cgmlst": core genome MLST（数百〜数千遺伝子座）
      - "wgmlst": whole genome MLST

    MLST:
      各ハウスキーピング遺伝子座のアリル番号の組み合わせで
      Sequence Type（ST）を決定する。
    """
    import subprocess

    if scheme == "mlst":
        cmd = f"mlst {assembly_fasta} --scheme {organism}"
        result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
        fields = result.stdout.strip().split("\t")
        typing = {
            "file": fields[0],
            "scheme": fields[1],
            "ST": fields[2],
            "alleles": fields[3:],
        }
    elif scheme == "cgmlst":
        cmd = f"chewbbaca AlleleCall -i {assembly_fasta} -g schema/ -o cgmlst_results/"
        subprocess.run(cmd, shell=True, check=True)
        typing = {"scheme": "cgMLST", "results_dir": "cgmlst_results/"}

    print(f"  Typing: ST={typing.get('ST', 'N/A')} ({scheme})")
    return typing

4. 系統解析・伝播推定

def phylogenetic_analysis(alignment_fasta, method="iqtree", model="GTR+G"):
    """
    病原体系統解析パイプライン。

    method:
      - "iqtree": IQ-TREE 2 — 最尤法（ModelFinder 自動モデル選択）
      - "raxml": RAxML-NG — 最尤法
      - "beast": BEAST 2 — ベイズ系統年代学

    アウトブレイク推定:
      - SNP 距離行列 → 最小スパニングツリー
      - tMRCA (最近共通祖先時間) 推定
    """
    import subprocess
    from Bio import Phylo

    if method == "iqtree":
        cmd = (f"iqtree2 -s {alignment_fasta} -m {model} "
               f"-bb 1000 -alrt 1000 -nt AUTO")
        subprocess.run(cmd, shell=True, check=True)
        tree = Phylo.read(f"{alignment_fasta}.treefile", "newick")

    return tree


def transmission_network(snp_matrix, max_snp_distance=10):
    """
    SNP 距離ベースの伝播ネットワーク推定。

    基準:
      - 直接伝播: SNP 距離 ≤ max_snp_distance
      - 近縁クラスタ: SNP 距離 ≤ 2 × max_snp_distance

    アルゴリズム:
      1. ペアワイズ SNP 距離計算
      2. 閾値以下のペアをエッジとして接続
      3. 最小スパニングツリーで伝播方向推定
    """
    import networkx as nx

    G = nx.Graph()
    samples = snp_matrix.index.tolist()
    G.add_nodes_from(samples)

    for i, s1 in enumerate(samples):
        for j, s2 in enumerate(samples):
            if i < j:
                dist = snp_matrix.iloc[i, j]
                if dist <= max_snp_distance:
                    G.add_edge(s1, s2, weight=dist, snp_distance=dist)

    mst = nx.minimum_spanning_tree(G)
    clusters = list(nx.connected_components(G))

    print(f"  Transmission: {G.number_of_edges()} links, "
          f"{len(clusters)} clusters")
    return G, mst, clusters

5. SIR / SEIR コンパートメントモデル

from scipy.integrate import odeint

def sir_model(y, t, beta, gamma, N):
    """
    SIR コンパートメントモデル。

    dS/dt = -β · S · I / N
    dI/dt =  β · S · I / N - γ · I
    dR/dt =  γ · I

    R₀ = β / γ (基本再生産数)
    """
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]


def seir_model(y, t, beta, sigma, gamma, N):
    """
    SEIR コンパートメントモデル（潜伏期あり）。

    dS/dt = -β · S · I / N
    dE/dt =  β · S · I / N - σ · E
    dI/dt =  σ · E - γ · I
    dR/dt =  γ · I

    σ: 潜伏期の逆数 (1/incubation_period)
    """
    S, E, I, R = y
    dSdt = -beta * S * I / N
    dEdt = beta * S * I / N - sigma * E
    dIdt = sigma * E - gamma * I
    dRdt = gamma * I
    return [dSdt, dEdt, dIdt, dRdt]


def run_epidemic_simulation(model="SIR", N=1e6, I0=10, R0=2.5,
                              gamma=1/10, sigma=1/5, days=180):
    """
    感染症拡大シミュレーション。

    Parameters:
        R0: 基本再生産数
        gamma: 回復率 (1/感染期間)
        sigma: 発症率 (1/潜伏期間、SEIR のみ)
        days: シミュレーション日数
    """
    beta = R0 * gamma
    t = np.linspace(0, days, days * 10)

    if model == "SIR":
        y0 = [N - I0, I0, 0]
        sol = odeint(sir_model, y0, t, args=(beta, gamma, N))
        df = pd.DataFrame(sol, columns=["S", "I", "R"])
    elif model == "SEIR":
        y0 = [N - I0, 0, I0, 0]
        sol = odeint(seir_model, y0, t, args=(beta, sigma, gamma, N))
        df = pd.DataFrame(sol, columns=["S", "E", "I", "R"])

    df["t"] = t
    peak_I = df["I"].max()
    peak_day = df.loc[df["I"].idxmax(), "t"]

    print(f"  {model}: R₀={R0:.1f}, peak infection={peak_I:.0f} at day {peak_day:.0f}")
    return df

References

Output Files

利用可能ツール

ToolUniverse SMCP 経由で利用可能な外部ツール。

参照スキル

依存パッケージ

• biopython, ete3, scipy, networkx, subprocess (fastp, bwa-mem2, freebayes, iqtree2)