856. 生信软件开发最佳实践¶

# 规范的Python生信工具结构

# 项目结构
# my_bio_tool/
# ├── pyproject.toml         # 项目配置（替代setup.py）
# ├── src/
# │   └── my_bio_tool/
# │       ├── __init__.py    # 包初始化
# │       ├── cli.py         # 命令行接口
# │       ├── io.py          # 输入输出
# │       ├── analysis.py    # 核心分析
# │       └── utils.py       # 工具函数
# ├── tests/
# │   ├── test_io.py         # IO测试
# │   └── test_analysis.py   # 分析测试
# └── README.md              # 说明文档

# === analysis.py 示例 ===
"""核心分析模块：差异丰度分析"""

import numpy as np                         # 数值计算
from scipy import stats                    # 统计检验
from dataclasses import dataclass          # 数据类
from typing import List, Optional          # 类型提示


@dataclass                                 # 用数据类存储结果
class DiffResult:
    """差异分析结果"""
    feature: str                           # 特征名称
    log2fc: float                          # log2 fold change
    pvalue: float                          # P值
    padj: float                            # 校正后P值


def calculate_log2fc(
    group1: np.ndarray,                    # 组1数据
    group2: np.ndarray,                    # 组2数据
    pseudocount: float = 1.0               # 伪计数（防止log(0)）
) -> float:
    """
    计算两组之间的log2 fold change。

    Parameters
    ----------
    group1 : np.ndarray
        对照组的丰度数据
    group2 : np.ndarray
        实验组的丰度数据
    pseudocount : float
        伪计数，防止除零错误，默认1.0

    Returns
    -------
    float
        log2 fold change值

    Examples
    --------
    >>> calculate_log2fc(np.array([10, 12, 8]), np.array([40, 45, 38]))
    1.95  # 约等于
    """
    mean1 = np.mean(group1) + pseudocount  # 组1均值+伪计数
    mean2 = np.mean(group2) + pseudocount  # 组2均值+伪计数
    return np.log2(mean2 / mean1)          # 计算log2FC


def differential_abundance(
    abundance_matrix: np.ndarray,          # 丰度矩阵(特征×样本)
    group_labels: List[str],               # 分组标签
    feature_names: List[str],              # 特征名称
    method: str = "wilcoxon"               # 统计方法
) -> List[DiffResult]:
    """
    差异丰度分析。

    Parameters
    ----------
    abundance_matrix : np.ndarray
        丰度矩阵，行=特征，列=样本
    group_labels : List[str]
        每个样本的分组标签
    feature_names : List[str]
        特征名称列表
    method : str
        统计方法：'wilcoxon'(默认) 或 'ttest'

    Returns
    -------
    List[DiffResult]
        差异分析结果列表
    """
    # 输入验证
    if abundance_matrix.shape[1] != len(group_labels):  # 检查维度
        raise ValueError(
            f"样本数不匹配: 矩阵{abundance_matrix.shape[1]}列 "
            f"vs 标签{len(group_labels)}个"
        )

    groups = list(set(group_labels))       # 获取唯一分组
    if len(groups) != 2:                   # 必须恰好2组
        raise ValueError(f"需要2个分组，得到{len(groups)}个")

    # 分组索引
    idx1 = [i for i, g in enumerate(group_labels) if g == groups[0]]
    idx2 = [i for i, g in enumerate(group_labels) if g == groups[1]]

    results = []                           # 结果列表
    pvalues = []                           # P值列表（用于多重检验校正）

    for i, feature in enumerate(feature_names):
        g1 = abundance_matrix[i, idx1]     # 组1数据
        g2 = abundance_matrix[i, idx2]     # 组2数据

        # 计算log2FC
        lfc = calculate_log2fc(g1, g2)     # fold change

        # 统计检验
        if method == "wilcoxon":           # Wilcoxon秩和检验
            _, pval = stats.mannwhitneyu(g1, g2, alternative="two-sided")
        elif method == "ttest":            # t检验
            _, pval = stats.ttest_ind(g1, g2)
        else:
            raise ValueError(f"未知方法: {method}")

        pvalues.append(pval)               # 收集P值
        results.append(DiffResult(         # 存储结果
            feature=feature,
            log2fc=lfc,
            pvalue=pval,
            padj=0.0                       # 后面填充
        ))

    # BH多重检验校正
    from statsmodels.stats.multitest import multipletests
    _, padj, _, _ = multipletests(pvalues, method="fdr_bh")

    for r, p in zip(results, padj):        # 填充校正P值
        r.padj = p

    return results

# === cli.py —— 命令行接口 ===
"""命令行接口：让工具像samtools/bwa一样好用"""

import argparse                            # 命令行参数解析
import sys                                 # 系统模块
import logging                             # 日志模块

def create_parser():
    """创建命令行参数解析器"""
    parser = argparse.ArgumentParser(
        prog="mybiotool",                  # 程序名
        description="微生物组差异丰度分析工具",  # 描述
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
示例用法:
  # 基本用法
  mybiotool diff -i abundance.tsv -m metadata.tsv -o results/

  # 指定分组列和方法
  mybiotool diff -i abundance.tsv -m metadata.tsv -g Group -method wilcoxon
        """
    )

    # 子命令
    subparsers = parser.add_subparsers(dest="command")

    # diff子命令：差异分析
    diff_parser = subparsers.add_parser("diff", help="差异丰度分析")
    diff_parser.add_argument(
        "-i", "--input",                   # 输入文件
        required=True,
        help="丰度表文件(TSV格式，行=特征，列=样本)"
    )
    diff_parser.add_argument(
        "-m", "--metadata",                # 元数据文件
        required=True,
        help="样本元数据文件(TSV格式)"
    )
    diff_parser.add_argument(
        "-g", "--group",                   # 分组列名
        default="Group",
        help="元数据中的分组列名(默认: Group)"
    )
    diff_parser.add_argument(
        "-o", "--output",                  # 输出目录
        default="./results",
        help="输出目录(默认: ./results)"
    )
    diff_parser.add_argument(
        "--method",                        # 统计方法
        choices=["wilcoxon", "ttest"],
        default="wilcoxon",
        help="统计检验方法(默认: wilcoxon)"
    )
    diff_parser.add_argument(
        "--padj-cutoff",                   # P值阈值
        type=float,
        default=0.05,
        help="校正P值阈值(默认: 0.05)"
    )
    diff_parser.add_argument(
        "-v", "--verbose",                 # 详细输出
        action="store_true",
        help="输出详细日志"
    )

    return parser


def main():
    """主函数入口"""
    parser = create_parser()               # 创建解析器
    args = parser.parse_args()             # 解析参数

    if args.command is None:               # 没有子命令
        parser.print_help()                # 打印帮助
        sys.exit(1)

    # 设置日志级别
    level = logging.DEBUG if args.verbose else logging.INFO
    logging.basicConfig(
        level=level,
        format="%(asctime)s [%(levelname)s] %(message)s"
    )

    if args.command == "diff":             # 差异分析
        logging.info(f"读取输入文件: {args.input}")
        # ... 调用分析函数 ...
        logging.info("分析完成！")

if __name__ == "__main__":
    main()

# === tests/test_analysis.py ===
"""单元测试：确保分析逻辑正确"""

import pytest                              # 测试框架
import numpy as np                         # 数值计算
from my_bio_tool.analysis import (
    calculate_log2fc,
    differential_abundance
)


def test_log2fc_basic():
    """测试基本的log2FC计算"""
    g1 = np.array([10, 10, 10])            # 组1: 均值10
    g2 = np.array([20, 20, 20])            # 组2: 均值20
    result = calculate_log2fc(g1, g2)      # 计算
    assert abs(result - 0.93) < 0.1        # log2(21/11) ≈ 0.93


def test_log2fc_equal():
    """测试两组相等时log2FC应约为0"""
    g1 = np.array([50, 50, 50])
    g2 = np.array([50, 50, 50])
    result = calculate_log2fc(g1, g2)
    assert abs(result) < 0.01             # 应该接近0


def test_diff_abundance_dimension_mismatch():
    """测试维度不匹配时应报错"""
    matrix = np.random.rand(5, 3)          # 5特征×3样本
    labels = ["A", "A", "B", "B"]          # 4个标签(不匹配)
    names = [f"feature_{i}" for i in range(5)]

    with pytest.raises(ValueError):        # 应该抛出ValueError
        differential_abundance(matrix, labels, names)


def test_diff_abundance_result_count():
    """测试结果数量应等于特征数"""
    np.random.seed(42)                     # 固定种子
    matrix = np.random.rand(10, 6)         # 10特征×6样本
    labels = ["A", "A", "A", "B", "B", "B"]
    names = [f"f{i}" for i in range(10)]

    results = differential_abundance(matrix, labels, names)
    assert len(results) == 10             # 应有10个结果

# .github/workflows/test.yml —— GitHub Actions CI/CD
name: Tests                               # 工作流名称

on:                                        # 触发条件
  push:
    branches: [main]                       # 推送到main分支
  pull_request:
    branches: [main]                       # PR到main分支

jobs:
  test:
    runs-on: ubuntu-latest                 # 运行环境
    strategy:
      matrix:
        python-version: ["3.10", "3.11", "3.12"]  # 多版本测试

    steps:
    - uses: actions/checkout@v4            # 检出代码

    - name: Set up Python                 # 设置Python
      uses: actions/setup-python@v5
      with:
        python-version: ${{ matrix.python-version }}

    - name: Install dependencies          # 安装依赖
      run: |
        pip install -e ".[test]"           # 安装包+测试依赖

    - name: Run tests                     # 运行测试
      run: |
        pytest tests/ -v --cov=src/       # 带覆盖率的测试

    - name: Check coverage               # 检查覆盖率
      run: |
        pytest --cov=src/ --cov-fail-under=80  # 覆盖率≥80%

# pyproject.toml —— 现代Python包配置（替代setup.py）
[build-system]
requires = ["setuptools>=68.0", "setuptools-scm"]
build-backend = "setuptools.backends._legacy:_Backend"

[project]
name = "my-bio-tool"
version = "0.1.0"
description = "微生物组差异丰度分析工具"
readme = "README.md"
license = {text = "MIT"}
requires-python = ">=3.10"
dependencies = [
    "numpy>=1.24",
    "scipy>=1.10",
    "pandas>=2.0",
    "statsmodels>=0.14",
]

[project.optional-dependencies]
test = ["pytest", "pytest-cov"]

[project.scripts]
mybiotool = "my_bio_tool.cli:main"     # 命令行入口点

# 代码规范
Python: PEP 8 + type hints + docstrings
R:      tidyverse style + roxygen2 文档
命令行:  argparse(Python) / optparse(R)

# 测试框架
Python: pytest + pytest-cov
R:      testthat + covr

# CI/CD平台
GitHub Actions: .github/workflows/
GitLab CI:      .gitlab-ci.yml
Bitbucket:      bitbucket-pipelines.yml

# 发布渠道
Python: PyPI (pip install)
R:      CRAN / Bioconductor
Conda:  Bioconda (生信专用频道)

# 开发工具
代码检查: ruff(Python) / lintr(R)
格式化:   black(Python) / styler(R)
类型检查: mypy(Python)
文档生成: Sphinx(Python) / pkgdown(R)