scikit-bio

# scikit-bio - Bioinformatics and Ecology scikit-bio provides the data structures and statistical methods needed for rigorous biological analysis. It excels in calculating alpha/beta diversity, performing ordination (PCoA), and handling complex phylogenetic trees. ## When to Use - Analyzing microbiome data (taxonomic composition, community structure). - Calculating ecological diversity metrics (Shannon, Simpson, UniFrac). - Performing ordination for visualization (PCoA, DCA). - High-level sequence manipulation (DNA, RNA, Protein with metadata). - Reading and writing phylogenetic trees (Newick format). - Pairwise and multiple sequence alignment analysis. - Statistical testing of community differences (PERMANOVA, ANOSIM). ## Reference Documentation **Official docs**: http://scikit-bio.org/ **Tutorials**: http://scikit-bio.org/docs/latest/ **Search patterns**: `skbio.sequence`, `skbio.stats.distance`, `skbio.diversity`, `skbio.stats.ordination` ## Core Principles ### Grammar of Biological Sequences Instead of using raw strings, scikit-bio uses typed objects (DNA, RNA, Protein). These objects know their alphabet, can handle quality scores (Phred), and support biological operations (transcription, translation). ### Distance Matrices Microbiome research often boils down to comparing samples. The DistanceMatrix object is central, allowing for easy indexing, sub-setting, and statistical testing (e.g., PERMANOVA). ### Diversity Metrics Provides a standardized implementation of hundreds of diversity metrics used in ecology, ensuring reproducibility across studies. ## Quick Reference ### Installation ```bash pip install scikit-bio ``` ### Standard Imports ```python 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 ``` ### Basic Pattern - Sequence Manipulation ```python 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 ``` ## Critical Rules ### ✅ DO - **Use Typed Sequences** - Always use DNA, RNA, or Protein instead of generic Sequence to enable alphabet-specific methods. - **Set Metadata** - Use the metadata and positional_metadata (for quality scores) attributes to keep data self-contained. - **Check Alphabet** - Use `.is_valid()` when importing data from untrusted sources to catch non-IUPAC characters. - **Project Distance Matrices** - Use `pcoa()` to visualize high-dimensional community data. - **Specify Reference for UniFrac** - When calculating UniFrac diversity, ensure your phylogenetic tree contains all taxa present in your abundance table. - **Validate Trees** - Use `skbio.tree.TreeNode` for tree manipulations as it provides robust traversal methods. ### ❌ DON'T - **Mix Sequence Types** - Don't try to align DNA with RNA objects. - **Ignore Quality Scores** - If you have FASTQ data, store Phred scores in positional_metadata. - **Manually Calculate Distance** - Don't use raw NumPy for biological distances; use `skbio.diversity.beta` to access UniFrac or Bray-Curtis. - **Forget ID Matching** - When performing ordination or PERMANOVA, ensure IDs in your distance matrix and metadata mapping match exactly. ## Anti-Patterns (NEVER) ```python 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') ``` ## Sequence Analysis (skbio.sequence) ### Advanced Sequence Operations ```python 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() ``` ## Diversity Analysis (skbio.diversity) ### Alpha Diversity (Within a sample) ```python 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}") ``` ### Beta Diversity (Between samples) ```python 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) ``` ## Phylogenetics (skbio.tree) ### Handling Trees ```python 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') ``` ## Statistics and Ordination ### PCoA (Principal Coordinates Analysis) ```python 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 ``` ### Community Testing (PERMANOVA) ```python 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']) ``` ## Practical Workflows ### 1. Microbiome Distance Visualization ```python 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. Sequence QC and Filtering ```python 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. Phylogenetic Distance Calculation ```python def get_tip_distances(tree): """Returns a distance matrix of all tips in a tree.""" dm = tree.tip_tip_distances() return dm ``` ## Performance Optimization ### Vectorized Diversity When calculating diversity for many samples, pass the entire 2D array to `beta.pw_distances` instead of looping through pairs. ### Tree Traversal Use `tree.preorder()` or `tree.postorder()` instead of manual recursion for much better performance on large (10,000+ tip) trees. ## Common Pitfalls and Solutions ### The "Missing ID" in Distance Matrix PERMANOVA and PCoA will fail if your metadata index and DistanceMatrix IDs don't match. ```python # ✅ 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 Memory Usage UniFrac calculation can be memory-intensive. For massive datasets, consider using the `skbio.diversity.beta.unweighted_unifrac` optimized versions or external tools like Unifrac-Binaries. ### IUPAC Ambiguity Standard DNA objects don't allow arbitrary characters. ```python # ❌ 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 is the mathematical heart of modern microbiome and evolutionary research. By enforcing strict biological typing and providing validated ecological metrics, it ensures that your biological insights are grounded in statistical rigor.