scikit-bio

scikit-bio - バイオインフォマティクスと生態学

scikit-bio は、厳密な生物学的解析に必要なデータ構造と統計手法を提供します。アルファ／ベータ多様性の計算、序列化（PCoA）、複雑な系統樹の処理に優れています。

使用する場面

• マイクロバイオームデータの解析（分類群の構成、群集構造）。
• 生態学的多様性指標（Shannon、Simpson、UniFrac）の計算。
• 可視化のための序列化（PCoA、DCA）の実行。
• 高度な配列操作（メタデータ付きの DNA、RNA、Protein）。
• 系統樹の読み込みと書き出し（Newick 形式）。
• ペアワイズおよび多重配列アラインメント解析。
• 群集の差異に関する統計的検定（PERMANOVA、ANOSIM）。

参考ドキュメント

公式ドキュメント: http://scikit-bio.org/
チュートリアル: http://scikit-bio.org/docs/latest/
検索パターン: skbio.sequence、skbio.stats.distance、skbio.diversity、skbio.stats.ordination

基本原則

生物学的配列の文法

生の文字列を使う代わりに、scikit-bio は型付きオブジェクト（DNA、RNA、Protein）を使用します。これらのオブジェクトはアルファベットを把握しており、クオリティスコア（Phred）を扱えるほか、生物学的な操作（転写、翻訳）をサポートします。

距離行列

マイクロバイオーム研究では、多くの場合、サンプル同士の比較が中心になります。DistanceMatrix オブジェクトは中心的な役割を果たし、インデックス指定、サブセット化、統計的検定（例：PERMANOVA）を簡単に行えます。

多様性指標

生態学で使用される数百の多様性指標を標準化された方法で実装し、研究間の再現性を確保します。

クイックリファレンス

インストール

pip install scikit-bio

標準インポート

import skbio
import numpy as np
import pandas as pd
from skbio import DNA, RNA, Protein, Sequence
from skbio.stats.distance import DistanceMatrix
from skbio.diversity import alpha, beta
from skbio.stats.ordination import pcoa

基本パターン - 配列操作

from skbio import DNA

# 1. Create a sequence with metadata
seq = DNA("ACC--GTT", metadata={'id': 'sample1', 'desc': 'gene A'})

# 2. Biological operations
rc = seq.reverse_complement()
degapped = seq.degap()

# 3. Validation
print(f"Is valid? {seq.is_valid()}") # Checks against DNA alphabet

重要なルール

✅ 推奨事項

• 型付き配列を使用する - アルファベット固有のメソッドを利用できるように、汎用の Sequence ではなく、必ず DNA、RNA、または Protein を使用します。
• メタデータを設定する - データを自己完結させるために、metadata 属性と positional_metadata 属性（クオリティスコア用）を使用します。
• アルファベットを確認する - 信頼できないソースからデータをインポートする場合は、.is_valid() を使用して IUPAC 以外の文字を検出します。
• 距離行列を射影する - 高次元の群集データを可視化するには、pcoa() を使用します。
• UniFrac の参照を指定する - UniFrac 多様性を計算する際は、系統樹に存在量テーブル内のすべての分類群が含まれていることを確認します。
• 樹を検証する - 系統樹の操作には skbio.tree.TreeNode を使用します。堅牢な走査メソッドを利用できます。

❌ 禁止事項

• 配列の型を混在させない - DNA オブジェクトと RNA オブジェクトをアラインメントしようとしないでください。
• クオリティスコアを無視しない - FASTQ データがある場合は、Phred スコアを positional_metadata に保存します。
• 距離を手動で計算しない - 生物学的な距離の計算に生の NumPy を使わず、skbio.diversity.beta を使用して UniFrac または Bray-Curtis にアクセスします。
• ID の照合を忘れない - 序列化または PERMANOVA を実行する際は、距離行列の ID とメタデータの対応が完全に一致していることを確認します。

アンチパターン（絶対に避ける）

from skbio import DNA

# ❌ BAD: Manual reverse complement with string logic
rc_str = seq_str[::-1].replace('A', 't').replace('T', 'a')... # Fragile

# ✅ GOOD: Built-in validated method
rc_seq = DNA(seq_str).reverse_complement()

# ❌ BAD: Using raw lists for community analysis
# data = [[1, 0, 5], [2, 1, 0]]
# dist = manual_bray_curtis(data)

# ✅ GOOD: Using DistanceMatrix
from skbio.stats.distance import DistanceMatrix
dm = DistanceMatrix(matrix_data, ids=['S1', 'S2', 'S3'])

# ❌ BAD: Stripping metadata to run scikit-learn
# vals = seq.values # Metadata lost!

# ✅ GOOD: Process within skbio or use standard IO
seq.write('output.fasta')

配列解析（skbio.sequence）

高度な配列操作

from skbio import DNA, Protein

# Sequence with quality scores
seq = DNA("ACGT", positional_metadata={'quality': [30, 35, 40, 20]})

# Slicing preserves metadata
sub = seq[1:3] 
print(sub.positional_metadata) # {'quality': [35, 40]}

# K-mer frequencies
freqs = seq.kmer_frequencies(k=2)

# Translation (DNA -> Protein)
protein = DNA("ATGCGA").translate()

多様性解析（skbio.diversity）

アルファ多様性（サンプル内）

from skbio.diversity import alpha

counts = [10, 0, 5, 2, 20] # Species abundances
otus = ['OTU1', 'OTU2', 'OTU3', 'OTU4', 'OTU5']

shannon = alpha.shannon(counts)
simpson = alpha.simpson(counts)
observed_otus = alpha.observed_otus(counts)

print(f"Shannon Index: {shannon:.3f}")

ベータ多様性（サンプル間）

from skbio.diversity import beta
import numpy as np

# Abundance table (samples x taxa)
data = np.array([[10, 20, 0], 
                 [5, 15, 2], 
                 [0, 1, 30]])
ids = ['Sample1', 'Sample2', 'Sample3']

# Bray-Curtis distance
bc_dm = beta.pw_distances(data, ids=ids, metric='braycurtis')

# Weighted UniFrac (Requires a tree)
# unifrac_dm = beta.weighted_unifrac(data, otu_ids, tree)

系統学（skbio.tree）

系統樹の処理

from skbio import TreeNode
from io import StringIO

# Load Newick tree
tree_str = "((A:0.1, B:0.2)C:0.3, D:0.4)E;"
tree = TreeNode.read(StringIO(tree_str))

# Traverse
for node in tree.tips():
    print(f"Leaf: {node.name}, dist: {node.length}")

# Find Lowest Common Ancestor
lca = tree.find_lca(['A', 'B'])
print(f"LCA of A and B is {lca.name}")

# Rooting
tree.root_at('D')

統計と序列化

PCoA（主座標分析）

from skbio.stats.ordination import pcoa
from skbio.stats.distance import DistanceMatrix

# Create DM
dm = DistanceMatrix([[0, 0.5, 0.8], [0.5, 0, 0.2], [0.8, 0.2, 0]], 
                    ids=['A', 'B', 'C'])

# Perform PCoA
results = pcoa(dm)

# View proportion of variance explained
print(results.proportion_explained)

# Access coordinates for plotting
coords = results.samples

群集検定（PERMANOVA）

from skbio.stats.distance import permanova

# metadata linking IDs to groups
metadata = pd.DataFrame({
    'BodySite': ['Gut', 'Gut', 'Skin', 'Skin']}, 
    index=['S1', 'S2', 'S3', 'S4'])

# Assuming dm is a DistanceMatrix of S1..S4
# results = permanova(dm, metadata, column='BodySite', permutations=999)
# print(results['p-value'])

実践的なワークフロー

1. マイクロバイオーム距離の可視化

def plot_microbiome_pcoa(abundance_df, metadata_df, group_col):
    """Full workflow: Abundance -> Distance -> PCoA -> Table."""
    from skbio.diversity import beta
    from skbio.stats.ordination import pcoa
    
    # 1. Calculate Bray-Curtis distance
    dm = beta.pw_distances(abundance_df.values, ids=abundance_df.index, metric='braycurtis')
    
    # 2. PCoA
    pc = pcoa(dm)
    
    # 3. Merge with metadata for plotting
    plot_data = pc.samples[['PC1', 'PC2']].join(metadata_df)
    
    return plot_data # Ready for Seaborn/Plotly

2. 配列の QC とフィルタリング

def filter_low_quality_sequences(sequences, min_avg_qual=30):
    """Filters DNA sequences based on Phred scores."""
    valid_seqs = []
    for s in sequences:
        # Assuming quality is in positional_metadata
        avg_q = np.mean(s.positional_metadata['quality'])
        if avg_q >= min_avg_qual:
            valid_seqs.append(s)
    return valid_seqs

3. 系統距離の計算

def get_tip_distances(tree):
    """Returns a distance matrix of all tips in a tree."""
    dm = tree.tip_tip_distances()
    return dm

パフォーマンスの最適化

ベクトル化された多様性計算

多数のサンプルの多様性を計算する場合は、ペアごとにループする代わりに、2D 配列全体を beta.pw_distances に渡します。

系統樹の走査

大規模な（先端が 10,000 個以上ある）系統樹では、手動で再帰処理を行うよりも、tree.preorder() または tree.postorder() を使用すると、パフォーマンスが大幅に向上します。

よくある落とし穴と解決策

距離行列の「ID 欠落」

メタデータのインデックスと DistanceMatrix の ID が一致しない場合、PERMANOVA と PCoA は失敗します。

# ✅ Solution: Ensure alignment
common_ids = set(dm.ids).intersection(metadata.index)
dm_sub = dm.filter(common_ids)
metadata_sub = metadata.loc[list(common_ids)]

UniFrac のメモリ使用量

UniFrac の計算はメモリを大量に消費する場合があります。大規模なデータセットでは、最適化された skbio.diversity.beta.unweighted_unifrac のバージョンや、Unifrac-Binaries のような外部ツールの利用を検討してください。

IUPAC の曖昧性

標準の DNA オブジェクトでは、任意の文字を使用できません。

# ❌ Problem: DNA("ACGN") raises error if 'N' isn't handled
# ✅ Solution: scikit-bio DNA supports IUPAC (N, R, Y, etc.) by default.
# But if you have non-standard characters:
from skbio import Sequence
custom = Sequence("ACG-X") # Generic sequence permits anything

scikit-bio は、現代のマイクロバイオーム研究と進化研究の数学的な中核です。厳密な生物学的型付けを適用し、検証済みの生態学的指標を提供することで、生物学的な知見が統計的な厳密さに基づくことを保証します。