scikit-bio

scikit-bio

概要

scikit-bio は、生物学データ分析のための包括的な Python ライブラリです。配列操作、アラインメント、系統発生学、微生物生態学、多変量統計を扱います。より広範な Python 科学技術スタックと連携する、専用のデータ構造（DNA、RNA、Protein、DistanceMatrix、TreeNode、TabularMSA）を提供します。

使用する場面

• マイクロバイオームデータのアルファ／ベータ多様性を計算し、PERMANOVA を実行する
• 距離行列から系統樹（NJ、UPGMA）を構築する
• 群集組成データに PCoA 順序付けを実行する
• 生物学的な形式（FASTA、FASTQ、Newick、BIOM）を読み書きする
• ペアワイズ配列アラインメント（Smith-Waterman、Needleman-Wunsch）を実行する
• 系統学的ベータ多様性の UniFrac 距離を計算する
• 生態学的な距離行列で統計検定（ANOSIM、Mantel）を実行する
• QIIME 2 のアーティファクトやマイクロバイオームのパイプラインを扱う
• ハイスループット NGS のアラインメント／バリアントコーリングには、代わりに STAR/BWA を使用してください
• タンパク質構造予測には、代わりに AlphaFold/ESMFold を使用してください

前提条件

pip install scikit-bio
# Optional: pip install biom-format  — HDF5 BIOM table support
# Optional: pip install matplotlib seaborn  — visualization

クイックスタート

import skbio
from skbio.diversity import alpha_diversity, beta_diversity
from skbio.stats.distance import permanova
from skbio.stats.ordination import pcoa
import numpy as np

# Sample OTU counts (samples × features)
counts = np.array([[10, 20, 30], [15, 25, 5], [5, 10, 40], [20, 5, 15]])
sample_ids = ['S1', 'S2', 'S3', 'S4']
grouping = ['control', 'control', 'treatment', 'treatment']

# Alpha diversity
shannon = alpha_diversity('shannon', counts, ids=sample_ids)
print(f"Shannon diversity: {shannon.values}")  # [1.09, 1.04, 0.94, 1.03]

# Beta diversity → PCoA → PERMANOVA
bc_dm = beta_diversity('braycurtis', counts, ids=sample_ids)
pcoa_results = pcoa(bc_dm)
results = permanova(bc_dm, grouping, permutations=999)
print(f"PERMANOVA p-value: {results['p-value']}")

コア API

1. 配列操作

from skbio import DNA, RNA, Protein

# Create and manipulate sequences
dna = DNA('ATCGATCGATCG', metadata={'id': 'gene1', 'description': 'test'})
rc = dna.reverse_complement()
rna = dna.transcribe()
protein = rna.translate()
print(f"DNA: {dna}, RC: {rc}, Protein: {protein}")

# Motif finding and k-mer analysis
motif_positions = dna.find_with_regex('ATG.{3}')
kmer_freqs = dna.kmer_frequencies(k=3)
print(f"3-mer frequencies: {dict(list(kmer_freqs.items())[:3])}")

# Sequence properties
print(f"Has degenerates: {dna.has_degenerates()}")
print(f"GC content: {dna.gc_content():.2f}")
degapped = dna.degap()  # Remove gap characters

# Metadata: sequence-level, positional, interval
dna = DNA('ATCGATCG', metadata={'id': 'seq1'},
          positional_metadata={'quality': [30, 35, 40, 38, 32, 36, 34, 33]})
dna.interval_metadata.add([(0, 4)], metadata={'type': 'promoter'})
print(f"Quality scores: {list(dna.positional_metadata['quality'])}")

2. 配列アラインメント

from skbio import DNA
from skbio.alignment import local_pairwise_align_ssw, TabularMSA

# Pairwise local alignment (Smith-Waterman via SSW)
seq1 = DNA('ACTCGATCGATCGATCGATCG')
seq2 = DNA('ATCGATCGATCGATCGATCGA')
alignment, score, start_end = local_pairwise_align_ssw(seq1, seq2)
print(f"Score: {score}, Positions: {start_end}")

# Multiple sequence alignment from file
msa = TabularMSA.read('alignment.fasta', constructor=DNA)
consensus = msa.consensus()
conservation = msa.conservation()
print(f"Consensus: {consensus[:20]}, Conservation: {conservation[:5]}")

3. 系統樹

from skbio import TreeNode, DistanceMatrix
from skbio.tree import nj, upgma

# Build tree from distance matrix
data = [[0, 5, 9, 9], [5, 0, 10, 10], [9, 10, 0, 8], [9, 10, 8, 0]]
dm = DistanceMatrix(data, ids=['A', 'B', 'C', 'D'])
tree = nj(dm)
print(tree.ascii_art())

# Tree operations
subtree = tree.shear(['A', 'B', 'C'])  # Prune to subset
tips = [node.name for node in tree.tips()]
lca = tree.lowest_common_ancestor(['A', 'B'])
print(f"Tips: {tips}, LCA children: {len(list(lca.children))}")

# Tree comparison
tree2 = upgma(dm)
rf_dist = tree.compare_rfd(tree2)
cophenetic_dm = tree.cophenetic_matrix()
print(f"Robinson-Foulds distance: {rf_dist}")

4. 多様性分析

from skbio.diversity import alpha_diversity, beta_diversity
import numpy as np

counts = np.array([[10, 20, 30, 0], [15, 25, 5, 10], [5, 10, 40, 2]])
sample_ids = ['S1', 'S2', 'S3']

# Alpha diversity (multiple metrics)
for metric in ['shannon', 'simpson', 'observed_otus', 'pielou_e']:
    alpha = alpha_diversity(metric, counts, ids=sample_ids)
    print(f"{metric}: {alpha.values.round(3)}")

# Beta diversity
bc_dm = beta_diversity('braycurtis', counts, ids=sample_ids)
jaccard_dm = beta_diversity('jaccard', counts, ids=sample_ids)
print(f"Bray-Curtis S1-S2: {bc_dm['S1', 'S2']:.3f}")

# Phylogenetic diversity (requires tree + OTU IDs)
from skbio.diversity import alpha_diversity, beta_diversity

faith_pd = alpha_diversity('faith_pd', counts, ids=sample_ids,
                           tree=tree, otu_ids=feature_ids)
unifrac_dm = beta_diversity('unweighted_unifrac', counts,
                            ids=sample_ids, tree=tree, otu_ids=feature_ids)
w_unifrac_dm = beta_diversity('weighted_unifrac', counts,
                              ids=sample_ids, tree=tree, otu_ids=feature_ids)
print(f"Faith PD: {faith_pd.values}")

5. 順序付け

from skbio.stats.ordination import pcoa, cca

# PCoA from distance matrix
pcoa_results = pcoa(bc_dm)
pc1 = pcoa_results.samples['PC1']
pc2 = pcoa_results.samples['PC2']
prop = pcoa_results.proportion_explained
print(f"PC1 explains {prop.iloc[0]:.1%}, PC2 explains {prop.iloc[1]:.1%}")

# CCA with environmental variables (constrained ordination)
# species_matrix: samples × species counts
# env_matrix: samples × environmental variables
cca_results = cca(species_matrix, env_matrix)
biplot_scores = cca_results.biplot_scores
print(f"Biplot scores shape: {biplot_scores.shape}")

# Save/load ordination results
pcoa_results.write('pcoa_results.txt')
loaded = skbio.OrdinationResults.read('pcoa_results.txt')

6. 統計検定

from skbio.stats.distance import permanova, anosim, permdisp, mantel

grouping = ['control', 'control', 'treatment']

# PERMANOVA — test group differences
perm_results = permanova(bc_dm, grouping, permutations=999)
print(f"PERMANOVA: F={perm_results['test statistic']:.3f}, "
      f"p={perm_results['p-value']:.4f}")

# ANOSIM — alternative group test
anosim_results = anosim(bc_dm, grouping, permutations=999)
print(f"ANOSIM: R={anosim_results['test statistic']:.3f}, "
      f"p={anosim_results['p-value']:.4f}")

# PERMDISP — test homogeneity of dispersions
permdisp_results = permdisp(bc_dm, grouping, permutations=999)

# Mantel test — correlation between distance matrices
r, p_value, n = mantel(bc_dm, jaccard_dm, method='pearson', permutations=999)
print(f"Mantel: r={r:.3f}, p={p_value:.4f}")

7. ファイル I/O

import skbio

# Read various formats
seq = skbio.DNA.read('input.fasta', format='fasta')
tree = skbio.TreeNode.read('tree.nwk')
dm = skbio.DistanceMatrix.read('distances.txt')

# Memory-efficient reading (generator for large files)
for seq in skbio.io.read('large.fasta', format='fasta', constructor=skbio.DNA):
    print(f"{seq.metadata['id']}: {len(seq)} bp")

# Write and convert formats
seq.write('output.fasta', format='fasta')
seqs = list(skbio.io.read('input.fastq', format='fastq', constructor=skbio.DNA))
skbio.io.write(seqs, format='fasta', into='output.fasta')

8. 距離行列

from skbio import DistanceMatrix
import numpy as np

# Create from array
data = np.array([[0, 1, 2], [1, 0, 3], [2, 3, 0]])
dm = DistanceMatrix(data, ids=['A', 'B', 'C'])

# Access and slice
print(f"A-B distance: {dm['A', 'B']}")
subset = dm.filter(['A', 'C'])
condensed = dm.condensed_form()  # scipy-compatible
df = dm.to_data_frame()
print(f"Shape: {df.shape}")

主要な概念

データ構造

• DNA/RNA/Protein：メタデータ（配列レベル、塩基ごとの位置情報、区間領域）を持つ、検証済みの生物学的配列オブジェクト
• TabularMSA：表形式の多重配列アラインメント。コンセンサス、保存性、位置エントロピーに対応
• TreeNode：トラバース、枝刈り、再ルート化、距離計算に対応する系統樹
• DistanceMatrix：ID に基づくインデックス付けが可能な対称距離行列。順序付けや統計機能と連携
• Table：マイクロバイオームのカウントデータ用の BIOM 特徴量テーブル（サンプル × 特徴量）。pandas、polars、AnnData と相互運用可能

系統学的多様性と非系統学的多様性

メトリックの種類	非系統学的	系統学的（系統樹が必要）
アルファ多様性	Shannon、Simpson、observed_otus、Pielou の均等度	Faith の PD
ベータ多様性	Bray-Curtis、Jaccard	非加重 UniFrac、加重 UniFrac

カウントは整数である必要があります

多様性関数には相対頻度ではなく、整数の存在量カウントが必要です。割合がある場合は、逆算して乗算します：

# Wrong: relative abundance
# counts = np.array([0.1, 0.5, 0.4])
# Correct: integer counts
counts = np.array([10, 50, 40])

よくあるワークフロー

ワークフロー 1：マイクロバイオームの多様性分析

import skbio
from skbio import TreeNode
from skbio.diversity import alpha_diversity, beta_diversity
from skbio.stats.ordination import pcoa
from skbio.stats.distance import permanova
import numpy as np

# 1. Load data
counts = np.loadtxt('otu_table.tsv', delimiter='\t', dtype=int, skiprows=1)
sample_ids = ['S1', 'S2', 'S3', 'S4', 'S5', 'S6']
feature_ids = ['OTU1', 'OTU2', 'OTU3', 'OTU4']
tree = TreeNode.read('phylogeny.nwk')
grouping = ['control', 'control', 'control', 'treatment', 'treatment', 'treatment']

# 2. Alpha diversity
shannon = alpha_diversity('shannon', counts, ids=sample_ids)
faith = alpha_diversity('faith_pd', counts, ids=sample_ids,
                        tree=tree, otu_ids=feature_ids)
print(f"Mean Shannon - Control: {shannon[:3].mean():.2f}, "
      f"Treatment: {shannon[3:].mean():.2f}")

# 3. Beta diversity + PCoA
unifrac_dm = beta_diversity('unweighted_unifrac', counts,
                            ids=sample_ids, tree=tree, otu_ids=feature_ids)
pcoa_results = pcoa(unifrac_dm)
print(f"PC1: {pcoa_results.proportion_explained.iloc[0]:.1%} variance")

# 4. Statistical testing
results = permanova(unifrac_dm, grouping, permutations=999)
print(f"PERMANOVA p={results['p-value']:.4f}")

ワークフロー 2：系統樹の構築

from skbio import DNA, DistanceMatrix
from skbio.tree import nj
from skbio.alignment import local_pairwise_align_ssw
import numpy as np

# 1. Read sequences
seqs = list(skbio.io.read('sequences.fasta', format='fasta', constructor=DNA))
n = len(seqs)
print(f"Loaded {n} sequences")

# 2. Compute pairwise distances (Hamming-like via alignment scores)
dist_data = np.zeros((n, n))
for i in range(n):
    for j in range(i+1, n):
        _, score, _ = local_pairwise_align_ssw(seqs[i], seqs[j])
        max_len = max(len(seqs[i]), len(seqs[j]))
        dist_data[i, j] = dist_data[j, i] = 1 - (score / max_len)

# 3. Build tree
ids = [s.metadata.get('id', f'seq{i}') for i, s in enumerate(seqs)]
dm = DistanceMatrix(dist_data, ids=ids)
tree = nj(dm)
print(tree.ascii_art())

# 4. Analyze tree
tree.write('output.nwk', format='newick')
cophenetic = tree.cophenetic_matrix()
print(f"Cophenetic matrix shape: {cophenetic.shape}")

ワークフロー 3：配列処理パイプライン

1. ジェネレーター（constructor=DNA を指定した skbio.io.read）を使用して FASTQ 配列を読み込む
2. positional_metadata['quality'] の品質スコアでフィルタリングする
3. find_with_regex() と配列スライスを使用してアダプターをトリミングする
4. find_with_regex('ATG.{3}') でモチーフを検索する
5. dna.transcribe().translate() で翻訳する
6. skbio.io.write() で出力 FASTA を書き込む

主要なパラメーター

パラメーター	関数	デフォルト	効果
metric	alpha/beta_diversity	必須	多様性メトリック名（例：'shannon'、'braycurtis'）
permutations	permanova/anosim/mantel	999	置換回数。値が大きいほど p 値の精度が高くなります
tree	alpha/beta_diversity	None	系統学的メトリック（Faith PD、UniFrac）で必須
otu_ids	alpha/beta_diversity	None	カウントを系統樹の末端に対応付ける特徴量 ID
method	mantel	'pearson'	相関手法：'pearson' または 'spearman'
constructor	io.read	None	パースに使用する配列クラス（DNA、RNA、Protein）
format	io.read/write	Auto	ファイル形式（fasta、fastq、newick など）
genetic_code	translate	1	NCBI 遺伝暗号（1=標準、11=細菌）
k	kmer_frequencies	必須	頻度計算に使用する k-mer の長さ

ベストプラクティス

1. 大きなファイルにはジェネレーターを使用する — skbio.io.read() はジェネレーターを返します。数百万件の配列に対して list() を使用するのは避けてください
2. 系統多様性には投影 CRS を使用する — 一致する tree と otu_ids を必ず指定し、tree.shear(feature_ids) で特徴量セットに合わせてツリーを枝刈りしてください
3. PERMANOVA と PERMDISP を組み合わせる — PERMANOVA は分散の違いに敏感です。グループの均質性を確認するため、PERMDISP を実行してください
4. 論文掲載に適した p 値には 999 回以上の置換を使用する。探索的解析では 99 回で十分です
5. 大きな特徴量テーブルには JSON より HDF5 BIOM 形式を使用する（I/O が高速で、ファイルサイズも小さくなります）
6. アンチパターン — 多様性関数での相対存在量：多様性関数には割合ではなく整数カウントが必要です。必要に応じて整数カウントに戻してください
7. アンチパターン — k-mer 解析で小さな k を使用する：k < 3 では識別力がほとんどありません。配列比較には k=5-7 を使用してください

よくあるレシピ

レシピ：多様性解析の前に希薄化する

from skbio.diversity import subsample_counts, alpha_diversity
import numpy as np

# Rarefy to minimum depth
min_depth = min(counts.sum(axis=1))
rarefied = np.array([subsample_counts(row, n=min_depth) for row in counts])
print(f"Rarefied to {min_depth} counts per sample")

# Calculate diversity on rarefied data
shannon = alpha_diversity('shannon', rarefied, ids=sample_ids)

レシピ：ベータ多様性の一部計算（大規模データセット）

from skbio.diversity import partial_beta_diversity
import itertools

# Only compute specific pairs (saves time on large matrices)
pairs = list(itertools.combinations(sample_ids[:10], 2))
partial_dm = partial_beta_diversity('braycurtis', counts,
                                    ids=sample_ids, id_pairs=pairs)
print(f"Computed {len(pairs)} pairwise distances")

レシピ：BIOM テーブルの操作

from skbio import Table
import numpy as np

# Read BIOM table
table = Table.read('feature_table.biom')
print(f"Samples: {table.shape[1]}, Features: {table.shape[0]}")

# Filter low-abundance features
filtered = table.filter(lambda row, id_, md: row.sum() > 10, axis='observation')

# Convert to pandas
df = table.to_dataframe()
# Normalize to relative abundance
rel_abundance = df.div(df.sum(axis=0), axis=1)

トラブルシューティング

問題	原因	解決策
ValueError: Ids must be unique	DistanceMatrix／配列内に重複 ID がある	ID を重複排除します：ids = list(dict.fromkeys(ids))
ValueError: Counts must be integers	多様性解析に相対存在量を渡している	整数カウントを使用します。元に戻すには (proportions * 1000).astype(int) を実行します
大きな FASTA ファイルでメモリエラーが発生する	すべての配列を一度に読み込んでいる	ジェネレーターを使用します：for seq in skbio.io.read(...)
ツリーの末端／OTU ID の不一致	系統多様性の計算に失敗する	ツリーを枝刈りします：tree = tree.shear(feature_ids)
PERMANOVA の p=0.001 に対し、PCoA ではグループが重なっている	分散に有意差がある	permdisp() を実行して確認します。PERMANOVA は位置と分散の両方を検定します
翻訳が正しくない	細菌にデフォルトの遺伝暗号（標準）が使われている	細菌／古細菌の配列には genetic_code=11 を設定します
NJ ツリーの構築が遅い	大きな n では O(n³) になる	分類群が 1000 を超える場合は GME または BME を使用します：from skbio.tree import gme

同梱リソース

• references/extended_api.md — 高度なアラインメントパラメーター（SSW、ギャップペナルティ、CIGAR）、ツリー構築アルゴリズム（NJ 対 UPGMA 対 GME/BME）、BIOM テーブルの操作、タンパク質エンベディング、QIIME 2、pandas、scikit-learn との統合パターンを網羅した拡張 API リファレンス

元の内容から移行していないもの：元の単一の api_reference.md（749 行）は、インラインで扱っていない機能に焦点を当てて references/extended_api.md にまとめました。省略した内容：Core API と重複する基本的な例、および詳細なトラブルシューティング（上記のトラブルシューティング表で説明済み）。

参考資料

• scikit-bio ドキュメント：https://scikit.bio/docs/latest/
• GitHub：https://github.com/scikit-bio/scikit-bio
• QIIME 2 フォーラム：https://forum.qiime2.org

関連スキル

• scanpy-scrna-seq — シングルセル RNA-seq 解析（データ領域は異なりますが、序列化の概念を共有しています）
• statistical-analysis — 生態統計を補完する一般的な統計検定
• matplotlib-scientific-plotting — PCoA プロットと多様性比較の可視化