Python数据分析（Python for Bioinformatics）

---

一句话说明

Python是生信数据处理和统计分析的主要编程语言，掌握pandas/scipy/matplotlib是分析工程师的基本要求。

---

核心知识点

2.1 pandas基础操作

常用操作速查表

操作	方法	示例	生信应用
读取TSV	pd.read_csv()	pd.read_csv("otu.tsv", sep="\t", index_col=0)	读取OTU丰度表
读取Excel	pd.read_excel()	pd.read_excel("metadata.xlsx")	读取样本信息表
写入文件	df.to_csv()	df.to_csv("result.tsv", sep="\t")	保存分析结果
查看前N行	df.head()	df.head(10)	快速预览数据
查看数据概况	df.describe()	df.describe()	了解数值分布
查看形状	df.shape	df.shape	确认样本数和特征数
查看列名	df.columns	df.columns.tolist()	获取所有列名
按条件筛选	df[条件]	df[df["abundance"] > 0.01]	筛选丰度>1%的物种
按列排序	df.sort_values()	df.sort_values("abundance", ascending=False)	按丰度降序排列
分组聚合	df.groupby()	df.groupby("phylum").sum()	按门水平汇总丰度
合并表格	pd.merge()	pd.merge(otu, meta, on="sample_id")	合并丰度表和元数据
纵向拼接	pd.concat()	pd.concat([df1, df2], axis=0)	合并多个样本的结果
横向拼接	pd.concat()	pd.concat([df1, df2], axis=1)	按列合并（如多个指标）
转置	df.T	otu_t = otu.T	将样本×物种转为物种×样本
透视表	df.pivot_table()	df.pivot_table(values="count", index="sample", columns="taxa")	重塑数据格式
缺失值处理	df.fillna()	df.fillna(0)	将缺失丰度填充为0
应用函数	df.apply()	df.apply(lambda x: x/x.sum(), axis=0)	按列归一化为相对丰度

代码示例：读取和处理OTU表

import pandas as pd
import numpy as np

# ============================================================
# 读取OTU丰度表
# 假设格式：行=OTU/物种，列=样本，值=reads计数
# ============================================================

# 读取TSV文件，第一列作为索引（OTU ID）
otu = pd.read_csv("otu_table.tsv", sep="\t", index_col=0)
print(f"OTU表形状: {otu.shape}")      # 输出: (500, 30) 表示500个OTU × 30个样本
print(f"样本数: {otu.shape[1]}")       # 列数 = 样本数
print(f"OTU数: {otu.shape[0]}")        # 行数 = OTU数

# 查看前5行
print(otu.head())

# ============================================================
# 计算相对丰度（relative abundance）
# 每个样本的所有OTU丰度除以该样本的总reads数
# ============================================================

# 方法一：apply函数（推荐，语义清晰）
rel_abundance = otu.apply(lambda col: col / col.sum(), axis=0)
# axis=0 表示对每一列（每个样本）操作
# 结果：每列之和 = 1.0

# 方法二：等价写法
rel_abundance = otu.div(otu.sum(axis=0), axis=1)

# 验证：每个样本的相对丰度之和应该为1
print(rel_abundance.sum(axis=0))      # 应该全是1.0

# ============================================================
# 筛选低丰度OTU
# 去除在所有样本中平均相对丰度<0.001(0.1%)的OTU
# ============================================================

mean_abundance = rel_abundance.mean(axis=1)         # 每个OTU在所有样本中的平均丰度
filtered_otu = rel_abundance[mean_abundance >= 0.001]  # 保留平均丰度>=0.1%的OTU
print(f"筛选前OTU数: {otu.shape[0]}")
print(f"筛选后OTU数: {filtered_otu.shape[0]}")

# ============================================================
# 按分类层级汇总
# 将OTU level的丰度汇总到门(Phylum)水平
# ============================================================

# 假设有分类信息表：OTU_ID → 门、纲、目、科、属、种
taxonomy = pd.read_csv("taxonomy.tsv", sep="\t", index_col=0)

# 合并OTU表和分类信息
otu_with_tax = otu.join(taxonomy)      # 按索引(OTU ID)合并

# 按门水平汇总
phylum_abundance = otu_with_tax.groupby("Phylum").sum()
# 结果：行=门，列=样本，值=该门下所有OTU的reads总和

# 保存结果
phylum_abundance.to_csv("phylum_abundance.tsv", sep="\t")

---

2.2 scipy统计检验

常用统计方法一览

检验方法	适用场景	scipy函数	前提假设
t检验	两组连续变量比较（正态分布）	scipy.stats.ttest_ind()	正态分布、方差齐性
Wilcoxon秩和检验	两组连续变量比较（非正态）	scipy.stats.mannwhitneyu()	无分布假设
配对t检验	配对样本比较（正态）	scipy.stats.ttest_rel()	配对、正态分布
Wilcoxon符号秩检验	配对样本比较（非正态）	scipy.stats.wilcoxon()	配对
卡方检验	分类变量关联性	scipy.stats.chi2_contingency()	期望频数>=5
Kruskal-Wallis检验	多组比较（非正态）	scipy.stats.kruskal()	无分布假设
ANOVA	多组比较（正态）	scipy.stats.f_oneway()	正态、方差齐性
Pearson相关	线性相关性	scipy.stats.pearsonr()	正态分布
Spearman相关	单调相关性	scipy.stats.spearmanr()	无分布假设

检验方法选择流程

数据类型？
├── 连续变量（如丰度值）
│   ├── 正态分布？
│   │   ├── 是 → 两组: t检验 / 多组: ANOVA
│   │   └── 否 → 两组: Wilcoxon(Mann-Whitney U) / 多组: Kruskal-Wallis
│   └── 配对数据？
│       ├── 正态 → 配对t检验
│       └── 非正态 → Wilcoxon符号秩
└── 分类变量（如有/无）
    └── 卡方检验 或 Fisher精确检验

生信数据的关键提醒：微生物丰度数据通常不满足正态分布，优先使用非参数检验（Wilcoxon）。

多重检验校正

from statsmodels.stats.multitest import multipletests

# ============================================================
# 多重检验校正（Multiple Testing Correction）
# 为什么需要：同时检验多个物种时，P值需要校正，否则假阳性率会很高
# 常用方法：Benjamini-Hochberg (FDR校正)
# ============================================================

# 假设我们对100个物种分别做了Wilcoxon检验
p_values = [0.001, 0.03, 0.04, 0.05, 0.15]  # 示意，实际为100个p值

# BH校正（控制FDR，False Discovery Rate）
reject, pvals_corrected, _, _ = multipletests(
    p_values,
    alpha=0.05,                        # 显著性阈值
    method='fdr_bh'                    # Benjamini-Hochberg方法
)
# reject: 布尔数组，True表示该检验显著
# pvals_corrected: 校正后的p值（q值）

# 筛选显著差异的物种
significant = [i for i, r in enumerate(reject) if r]
print(f"校正前显著(p<0.05): {sum(p < 0.05 for p in p_values)}个")
print(f"校正后显著(q<0.05): {len(significant)}个")

---

2.3 scikit-learn基础

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import pandas as pd
import numpy as np

# ============================================================
# 数据标准化（Standardization）
# 目的：消除量纲差异，让不同特征在同一尺度上比较
# ============================================================

# 读取丰度矩阵：行=样本，列=物种
abundance = pd.read_csv("abundance_samples_by_species.tsv", sep="\t", index_col=0)

scaler = StandardScaler()              # 创建标准化器
abundance_scaled = scaler.fit_transform(abundance)
# fit_transform: 计算均值和标准差 → 对每列做 (x - mean) / std
# 结果：每列均值=0，标准差=1

# ============================================================
# PCA降维（Principal Component Analysis，主成分分析）
# 目的：将高维数据投影到低维空间，用于可视化和探索
# ============================================================

pca = PCA(n_components=2)              # 降到2维
pca_result = pca.fit_transform(abundance_scaled)
# pca_result 形状: (样本数, 2)

# 获取解释方差比例（每个主成分解释了多少数据变异）
print(f"PC1解释方差: {pca.explained_variance_ratio_[0]:.2%}")
print(f"PC2解释方差: {pca.explained_variance_ratio_[1]:.2%}")

# 转为DataFrame方便后续绘图
pca_df = pd.DataFrame(
    pca_result,
    columns=["PC1", "PC2"],
    index=abundance.index               # 保留样本名
)

# ============================================================
# 随机森林分类（Random Forest）
# 目的：根据微生物组成预测样本分组（如健康vs疾病）
# ============================================================

# 准备数据
X = abundance_scaled                    # 特征矩阵（样本×物种）
y = metadata["group"].values            # 标签（如 "healthy", "disease"）

# 创建随机森林分类器
rf = RandomForestClassifier(
    n_estimators=100,                   # 100棵决策树
    random_state=42,                    # 随机种子（保证结果可重复）
    n_jobs=-1                           # 使用所有CPU核
)

# 交叉验证评估模型效果
scores = cross_val_score(
    rf, X, y,
    cv=5,                               # 5折交叉验证
    scoring='accuracy'                  # 准确率
)
print(f"5折交叉验证准确率: {scores.mean():.2%} (+/- {scores.std():.2%})")

# 训练模型并查看特征重要性
rf.fit(X, y)
importances = pd.Series(
    rf.feature_importances_,            # 每个特征的重要性得分
    index=abundance.columns             # 对应的物种名
).sort_values(ascending=False)

print("Top 10重要物种:")
print(importances.head(10))

---

2.4 可视化（matplotlib/seaborn）

基本图表类型

图表类型	函数	适用场景
柱状图	plt.bar() / sns.barplot()	组间均值比较
箱线图	sns.boxplot()	数据分布比较
散点图	plt.scatter() / sns.scatterplot()	PCoA/PCA结果展示
热图	sns.heatmap()	相关性矩阵、丰度热图
堆叠柱状图	df.plot(kind="bar", stacked=True)	物种组成展示
小提琴图	sns.violinplot()	数据分布形状比较

---

实战代码

计算Alpha多样性

import pandas as pd
import numpy as np

# ============================================================
# Alpha多样性指数计算
# Alpha多样性衡量单个样本内的物种多样性
# ============================================================

def shannon_diversity(counts):
    """
    计算Shannon多样性指数
    Shannon Index = -sum(pi * ln(pi))
    其中 pi = 第i个物种的相对丰度

    参数:
        counts: 一维数组，每个物种的reads计数
    返回:
        float: Shannon指数值（越高表示多样性越高）
    """
    # 去除0值（0不参与计算，因为 0*log(0) 在数学上定义为0）
    counts = np.array(counts)
    counts = counts[counts > 0]         # 只保留>0的值

    # 计算相对丰度 pi
    total = counts.sum()
    proportions = counts / total        # pi = ni / N

    # 计算Shannon指数
    shannon = -np.sum(proportions * np.log(proportions))
    return shannon


def simpson_diversity(counts):
    """
    计算Simpson多样性指数
    Simpson Index = 1 - sum(pi^2)
    取值范围 [0, 1)，越接近1表示多样性越高

    参数:
        counts: 一维数组，每个物种的reads计数
    返回:
        float: Simpson指数值
    """
    counts = np.array(counts)
    counts = counts[counts > 0]

    total = counts.sum()
    proportions = counts / total

    simpson = 1 - np.sum(proportions ** 2)
    return simpson


def observed_species(counts):
    """
    计算观测到的物种数（Observed Species / Richness）
    即计数>0的物种数量
    """
    counts = np.array(counts)
    return np.sum(counts > 0)


# ============================================================
# 应用到OTU表
# ============================================================

# 读取OTU表（行=OTU, 列=样本）
otu = pd.read_csv("otu_table.tsv", sep="\t", index_col=0)

# 对每个样本计算多样性指数
alpha_diversity = pd.DataFrame({
    "Shannon": otu.apply(shannon_diversity, axis=0),     # 对每列（样本）计算
    "Simpson": otu.apply(simpson_diversity, axis=0),
    "Observed_Species": otu.apply(observed_species, axis=0)
})

print(alpha_diversity.head())
# 输出示例:
#           Shannon  Simpson  Observed_Species
# sample_01    3.42     0.93              187
# sample_02    2.89     0.88              156

# 保存结果
alpha_diversity.to_csv("alpha_diversity.tsv", sep="\t")

Bray-Curtis距离矩阵计算

import pandas as pd
import numpy as np
from scipy.spatial.distance import pdist, squareform

# ============================================================
# Bray-Curtis距离（Beta多样性）
# Bray-Curtis衡量两个样本之间物种组成的差异
# 取值范围 [0, 1]：0=完全相同，1=完全不同
# ============================================================

# 读取相对丰度表（行=样本, 列=物种）
# 注意：这里要转置，因为pdist要求行是样本
rel_abundance = pd.read_csv("relative_abundance.tsv", sep="\t", index_col=0)
rel_abundance_T = rel_abundance.T      # 转置：行变成样本

# 计算Bray-Curtis距离
# pdist: pairwise distances（成对距离）
bc_condensed = pdist(rel_abundance_T, metric="braycurtis")
# 返回压缩格式的距离向量

# 转换为方阵（距离矩阵）
bc_matrix = squareform(bc_condensed)

# 转为DataFrame，添加样本名
bc_df = pd.DataFrame(
    bc_matrix,
    index=rel_abundance_T.index,       # 行名=样本名
    columns=rel_abundance_T.index      # 列名=样本名
)

print("Bray-Curtis距离矩阵:")
print(bc_df.head())

# 保存距离矩阵
bc_df.to_csv("bray_curtis_distance.tsv", sep="\t")

差异检验（Wilcoxon）

import pandas as pd
import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

# ============================================================
# Wilcoxon秩和检验（Mann-Whitney U检验）
# 比较两组样本之间每个物种的丰度差异
# 适用于：非正态分布的微生物丰度数据
# ============================================================

# 读取数据
rel_abundance = pd.read_csv("relative_abundance.tsv", sep="\t", index_col=0)
metadata = pd.read_csv("metadata.tsv", sep="\t", index_col=0)

# 按分组拆分样本
# 假设metadata有一列"group"，值为"healthy"或"disease"
group1_samples = metadata[metadata["group"] == "healthy"].index
group2_samples = metadata[metadata["group"] == "disease"].index

# 对每个物种进行Wilcoxon检验
results = []

for species in rel_abundance.index:    # 遍历每个物种
    # 提取两组的丰度值
    group1_values = rel_abundance.loc[species, group1_samples].values
    group2_values = rel_abundance.loc[species, group2_samples].values

    # Mann-Whitney U检验（即Wilcoxon秩和检验）
    try:
        stat, pvalue = stats.mannwhitneyu(
            group1_values,
            group2_values,
            alternative='two-sided'     # 双侧检验
        )
    except ValueError:
        # 如果所有值都相同，检验会报错
        pvalue = 1.0
        stat = 0

    # 计算效应大小（两组均值的差异）
    mean_g1 = np.mean(group1_values)
    mean_g2 = np.mean(group2_values)
    log2fc = np.log2((mean_g2 + 1e-10) / (mean_g1 + 1e-10))
    # 加1e-10避免除以0
    # log2FC > 0: disease组丰度更高
    # log2FC < 0: healthy组丰度更高

    results.append({
        "species": species,
        "mean_healthy": mean_g1,
        "mean_disease": mean_g2,
        "log2FC": log2fc,
        "p_value": pvalue
    })

# 转为DataFrame
results_df = pd.DataFrame(results)

# 多重检验校正（BH方法）
reject, q_values, _, _ = multipletests(
    results_df["p_value"],
    alpha=0.05,
    method="fdr_bh"                     # Benjamini-Hochberg
)
results_df["q_value"] = q_values       # 添加校正后的q值
results_df["significant"] = reject      # 是否显著

# 按q值排序
results_df = results_df.sort_values("q_value")

# 筛选显著差异物种
sig_species = results_df[results_df["significant"]]
print(f"显著差异物种数: {len(sig_species)} / {len(results_df)}")
print(sig_species.head(10))

# 保存结果
results_df.to_csv("differential_abundance.tsv", sep="\t", index=False)

堆叠柱状图绘制

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['font.sans-serif'] = ['SimHei']    # 显示中文
matplotlib.rcParams['axes.unicode_minus'] = False       # 显示负号

# ============================================================
# 堆叠柱状图（Stacked Bar Chart）
# 展示各样本的门水平物种组成
# ============================================================

# 读取门水平相对丰度表（行=门, 列=样本）
phylum_rel = pd.read_csv("phylum_relative_abundance.tsv", sep="\t", index_col=0)

# 只保留Top 10的门，其余归为"Others"
top_phyla = phylum_rel.mean(axis=1).nlargest(10).index  # 按平均丰度取前10
others = phylum_rel.loc[~phylum_rel.index.isin(top_phyla)].sum()  # 其余汇总
plot_data = pd.concat([
    phylum_rel.loc[top_phyla],
    pd.DataFrame(others, columns=["Others"]).T           # 添加Others行
])

# 转置：行=样本，列=门（pandas绘制柱状图的格式要求）
plot_data_T = plot_data.T

# 定义配色（选择区分度高的颜色）
colors = [
    '#E64B35', '#4DBBD5', '#00A087', '#3C5488', '#F39B7F',
    '#8491B4', '#91D1C2', '#DC0000', '#7E6148', '#B09C85',
    '#CCCCCC'                                              # Others用灰色
]

# 绘图
fig, ax = plt.subplots(figsize=(14, 6))

plot_data_T.plot(
    kind='bar',                         # 柱状图
    stacked=True,                       # 堆叠
    ax=ax,                              # 指定Axes
    color=colors[:len(plot_data_T.columns)],  # 颜色
    edgecolor='white',                  # 白色边框
    linewidth=0.5                       # 边框宽度
)

# 美化
ax.set_ylabel("Relative Abundance", fontsize=12)
ax.set_xlabel("")
ax.set_title("Phylum-level Composition", fontsize=14)
ax.legend(
    title="Phylum",
    bbox_to_anchor=(1.02, 1),           # 图例放右边
    loc='upper left',
    fontsize=9
)
ax.set_ylim(0, 1)                       # y轴范围 [0, 1]
ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha='right')
plt.tight_layout()

# 保存图片
plt.savefig("phylum_barplot.pdf", dpi=300, bbox_inches='tight')
plt.savefig("phylum_barplot.png", dpi=300, bbox_inches='tight')
plt.show()

PCoA绘图

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform
from skbio.stats.ordination import pcoa    # scikit-bio包
# 安装: conda install -c conda-forge scikit-bio

# ============================================================
# PCoA分析与可视化（Principal Coordinates Analysis）
# PCoA = 主坐标分析，基于距离矩阵进行降维
# 常用于展示样本间的beta多样性
# ============================================================

# 读取数据
rel_abundance = pd.read_csv("relative_abundance.tsv", sep="\t", index_col=0)
metadata = pd.read_csv("metadata.tsv", sep="\t", index_col=0)

# 计算Bray-Curtis距离矩阵
from skbio import DistanceMatrix

rel_T = rel_abundance.T                # 转置：行=样本
bc_condensed = pdist(rel_T, metric="braycurtis")
dm = DistanceMatrix(squareform(bc_condensed), ids=rel_T.index.tolist())

# PCoA分析
pcoa_result = pcoa(dm)

# 提取前两个主坐标轴的坐标
pcoa_coords = pcoa_result.samples[["PC1", "PC2"]]
# 解释方差比例
pc1_var = pcoa_result.proportion_explained["PC1"]
pc2_var = pcoa_result.proportion_explained["PC2"]

# 合并分组信息
pcoa_df = pcoa_coords.join(metadata[["group"]])

# ============================================================
# 绘制PCoA散点图
# ============================================================

fig, ax = plt.subplots(figsize=(8, 6))

# 定义分组颜色
group_colors = {"healthy": "#00A087", "disease": "#E64B35"}
groups = pcoa_df["group"].unique()

for group in groups:
    subset = pcoa_df[pcoa_df["group"] == group]
    ax.scatter(
        subset["PC1"],
        subset["PC2"],
        c=group_colors[group],
        label=group,
        s=80,                           # 点大小
        alpha=0.7,                      # 透明度
        edgecolors='white',             # 白色边框
        linewidth=0.5
    )

# 添加坐标轴标签（含解释方差）
ax.set_xlabel(f"PC1 ({pc1_var:.1%})", fontsize=12)
ax.set_ylabel(f"PC2 ({pc2_var:.1%})", fontsize=12)
ax.set_title("PCoA (Bray-Curtis)", fontsize=14)
ax.legend(title="Group", fontsize=10)

# 添加参考线
ax.axhline(y=0, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
ax.axvline(x=0, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)

plt.tight_layout()
plt.savefig("pcoa_plot.pdf", dpi=300, bbox_inches='tight')
plt.savefig("pcoa_plot.png", dpi=300, bbox_inches='tight')
plt.show()

# ============================================================
# PERMANOVA检验（检验组间距离是否显著不同）
# ============================================================
from skbio.stats.distance import permanova

permanova_result = permanova(dm, metadata["group"], permutations=999)
print(f"PERMANOVA结果:")
print(f"  test statistic (pseudo-F): {permanova_result['test statistic']:.3f}")
print(f"  p-value: {permanova_result['p-value']:.4f}")
# p < 0.05 表示两组的微生物组成有显著差异

相关性热图

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats

# ============================================================
# 物种间相关性热图
# 用Spearman相关分析物种之间的共现关系
# ============================================================

# 读取丰度数据（行=样本, 列=物种）
abundance = pd.read_csv("top_species_abundance.tsv", sep="\t", index_col=0)

# 计算Spearman相关系数矩阵
corr_matrix = abundance.corr(method="spearman")

# 计算p值矩阵（判断相关性是否显著）
n_species = abundance.shape[1]
p_matrix = pd.DataFrame(
    np.ones((n_species, n_species)),
    index=abundance.columns,
    columns=abundance.columns
)

for i in range(n_species):
    for j in range(i+1, n_species):
        rho, pval = stats.spearmanr(
            abundance.iloc[:, i],
            abundance.iloc[:, j]
        )
        p_matrix.iloc[i, j] = pval
        p_matrix.iloc[j, i] = pval

# 创建掩码：只显示下三角
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))

# 绘制热图
fig, ax = plt.subplots(figsize=(10, 8))

sns.heatmap(
    corr_matrix,
    mask=mask,                          # 只显示下三角
    cmap="RdBu_r",                      # 红蓝配色（正相关=红，负相关=蓝）
    center=0,                           # 以0为中心
    vmin=-1, vmax=1,                    # 范围[-1, 1]
    square=True,                        # 方格
    linewidths=0.5,                     # 格子边框
    cbar_kws={"shrink": 0.8, "label": "Spearman rho"},
    annot=True,                         # 显示数值
    fmt=".2f",                          # 保留2位小数
    ax=ax
)

# 在显著相关的格子上标记星号
for i in range(n_species):
    for j in range(i):
        if p_matrix.iloc[i, j] < 0.05:
            ax.text(
                j + 0.5, i + 0.8, "*",
                ha='center', va='center',
                color='black', fontsize=10
            )

ax.set_title("Species Correlation (Spearman)", fontsize=14)
plt.tight_layout()
plt.savefig("correlation_heatmap.pdf", dpi=300, bbox_inches='tight')
plt.show()

---

面试常问点

★ pandas的DataFrame和Series有什么区别？

参考答案：

Series 是一维数据结构，类似于一列数据，有索引（index）和值（values）。可以理解为 Excel 表格中的一列。

DataFrame 是二维数据结构，由多个 Series 组成，有行索引（index）和列名（columns）。可以理解为一张完整的 Excel 表格。

# Series：一维
s = pd.Series([3.42, 2.89, 4.01], index=["S1", "S2", "S3"])

# DataFrame：二维
df = pd.DataFrame({
    "Shannon": [3.42, 2.89, 4.01],
    "Simpson": [0.93, 0.88, 0.95]
}, index=["S1", "S2", "S3"])

从 DataFrame 中取一列就得到一个 Series：df["Shannon"] 返回 Series。

---

★ 如何用Python计算Shannon多样性指数？

参考答案：

Shannon 多样性指数的公式是：H' = -sum(pi * ln(pi))，其中 pi 是第 i 个物种的相对丰度。

import numpy as np

def shannon(counts):
    counts = np.array(counts)
    counts = counts[counts > 0]        # 去除0值
    proportions = counts / counts.sum() # 计算相对丰度
    return -np.sum(proportions * np.log(proportions))

Shannon 指数越高，表示群落多样性越高。一般微生物群落的 Shannon 指数在 2-5 之间。在我实习做 16S 分析时，会用 QIIME2 自动计算，但也会用 Python 自己验证结果。

---

★ Wilcoxon和t检验的区别？什么时候用哪个？

参考答案：

核心区别在于对数据分布的假设：

对比项	t检验	Wilcoxon秩和检验
分布假设	要求正态分布	不要求，非参数方法
比较对象	比较均值差异	比较秩次（排名）差异
敏感性	对异常值敏感	对异常值稳健
检验效能	正态数据下效能更高	非正态数据更可靠

在生信中，微生物丰度数据几乎都不是正态分布（通常是零膨胀的、右偏的），所以优先使用 Wilcoxon检验（也叫 Mann-Whitney U 检验）。

在我的实习项目中，比较健康组和疾病组的菌群差异时，就是用的 Wilcoxon 秩和检验，然后用 BH 方法做多重检验校正。

---

★ PCA和PCoA的区别？

参考答案：

对比项	PCA	PCoA
全称	Principal Component Analysis	Principal Coordinates Analysis
输入	原始数据矩阵（样本x特征）	距离矩阵（样本x样本）
距离度量	固定使用欧氏距离	可以使用任意距离（如Bray-Curtis）
适用场景	基因表达、环境因子	微生物群落组成

为什么微生物分析常用PCoA而不是PCA？

因为微生物丰度数据是"成分数据"（compositional data），使用欧氏距离不合适。而 Bray-Curtis 距离专门设计用于物种丰度数据的比较，它考虑了物种的有无和丰度差异。所以我们先计算 Bray-Curtis 距离矩阵，再用 PCoA 降维可视化。

在 QIIME2 中，beta 多样性分析默认就是用 Bray-Curtis + PCoA 的组合。

---

★ 你在项目中用Python做了什么分析？

参考答案：

在百迈客实习期间，我主要用 Python 做以下分析：

1. 数据预处理：用 pandas 读取和整理 OTU 丰度表、合并样本元数据、计算相对丰度、按分类层级汇总（如从 OTU 汇总到门水平）。

2. Alpha 多样性计算：用 numpy 实现 Shannon、Simpson 指数的计算（虽然 QIIME2 也能算，但有时需要自定义分析或验证）。

3. 差异分析：用 scipy.stats 的 Mann-Whitney U 检验比较两组间的物种丰度差异，用 statsmodels 做 BH 多重检验校正。

4. 可视化：用 matplotlib 和 seaborn 绘制堆叠柱状图（物种组成）、PCoA 散点图（beta 多样性）、箱线图（alpha 多样性组间比较）、热图（物种相关性）。

5. 批量处理：写 Python 脚本批量整理客户交付的分析报告数据，自动生成标准化的结果表格。

---

易错/易混淆点

1. merge vs join vs concat

import pandas as pd

# ============================================================
# 三种合并方式的区别
# ============================================================

# ---------- merge：按指定列合并（类似SQL的JOIN） ----------
# 最灵活，可以指定合并键
result = pd.merge(
    otu_table,                          # 左表
    metadata,                           # 右表
    on="sample_id",                     # 合并键（两表共有的列名）
    how="inner"                         # inner/left/right/outer
)
# 也可以用不同列名合并：
# pd.merge(df1, df2, left_on="id", right_on="sample")

# ---------- join：按索引合并 ----------
# 基于DataFrame的index进行合并，更简洁
result = otu_table.join(metadata)       # 默认按index合并
# 等价于 pd.merge(otu_table, metadata, left_index=True, right_index=True)

# ---------- concat：简单拼接 ----------
# 纵向拼接（追加行）
all_samples = pd.concat([batch1, batch2], axis=0)
# axis=0: 上下拼接（追加行）

# 横向拼接（追加列）
combined = pd.concat([otu_table, diversity_index], axis=1)
# axis=1: 左右拼接（追加列）

# 记忆口诀：
# merge = 按列合并（灵活），join = 按索引合并（简便），concat = 简单拼接（堆叠）

2. axis=0 vs axis=1

import pandas as pd
import numpy as np

# ============================================================
# axis参数的含义（新手最容易搞混的地方！）
# ============================================================

# 假设有一个 3x4 的DataFrame
df = pd.DataFrame(
    np.random.rand(3, 4),               # 3行4列
    index=["sample_A", "sample_B", "sample_C"],    # 行 = 样本
    columns=["sp1", "sp2", "sp3", "sp4"]           # 列 = 物种
)

# axis=0：沿着行的方向操作（结果是对每列操作）
df.sum(axis=0)        # 对每一列求和 → 得到4个值（每个物种的总丰度）
df.mean(axis=0)       # 每列的均值 → 每个物种在所有样本中的平均丰度
df.drop("sample_A", axis=0)  # 删除一行

# axis=1：沿着列的方向操作（结果是对每行操作）
df.sum(axis=1)        # 对每一行求和 → 得到3个值（每个样本的总reads数）
df.mean(axis=1)       # 每行的均值 → 每个样本的平均物种丰度
df.drop("sp1", axis=1)       # 删除一列

# 记忆技巧：
# axis=0 → 上下方向（垂直，像0的形状） → 操作行 → 结果按列
# axis=1 → 左右方向（水平，像1的形状） → 操作列 → 结果按行

# 生信实战：计算相对丰度
# 每个样本的丰度除以该样本的总丰度（行=样本，列=物种）
relative = df.div(df.sum(axis=1), axis=0)  # axis=1 对行求和，axis=0 按行除

3. 深拷贝 vs 浅拷贝

import pandas as pd

# ============================================================
# 深拷贝 vs 浅拷贝
# 这是Python中非常重要的概念，搞混了会修改原始数据！
# ============================================================

otu_original = pd.read_csv("otu_table.tsv", sep="\t", index_col=0)

# ---------- 浅拷贝（引用）—— 危险！ ----------
otu_ref = otu_original                  # 这不是拷贝！只是给同一个对象起了新名字
otu_ref.iloc[0, 0] = 999               # 修改otu_ref...
print(otu_original.iloc[0, 0])         # 输出: 999 ← 原始数据也被修改了！

# ---------- 深拷贝 —— 安全 ----------
otu_copy = otu_original.copy()          # 创建完全独立的副本
otu_copy.iloc[0, 0] = 0                # 修改副本
print(otu_original.iloc[0, 0])         # 原始数据不受影响

# 生信场景：做数据筛选/过滤时，一定要用 .copy()
filtered = otu_original[otu_original.mean(axis=1) > 10].copy()
# 如果不加 .copy()，后续修改filtered可能触发SettingWithCopyWarning

4. loc vs iloc

import pandas as pd

# ============================================================
# loc vs iloc
# loc = 按标签（名称）索引
# iloc = 按位置（整数）索引
# ============================================================

df = pd.DataFrame(
    {"Shannon": [3.42, 2.89, 4.01], "Simpson": [0.93, 0.88, 0.95]},
    index=["sample_A", "sample_B", "sample_C"]
)

# ---------- loc：按标签名称 ----------
df.loc["sample_A", "Shannon"]          # 返回: 3.42
df.loc["sample_A":"sample_B"]          # 返回前两行（包含末尾！与Python切片不同！）
df.loc[df["Shannon"] > 3, "Shannon"]   # 条件筛选

# ---------- iloc：按整数位置 ----------
df.iloc[0, 0]                          # 返回: 3.42（第0行第0列）
df.iloc[0:2]                           # 返回前两行（不包含末尾，标准Python切片）
df.iloc[:, 0]                          # 所有行的第0列

# 记忆：
# loc = label-based（标签）→ loc 里面用"名称"
# iloc = integer-based（整数）→ iloc 里面用"数字"

# 关键区别：loc切片包含末尾，iloc切片不包含末尾！
df.loc["sample_A":"sample_C"]          # 包含sample_C（3行）
df.iloc[0:3]                           # 不包含index=3（也是3行，但规则不同）

---

常用Python包安装

# ============================================================
# 创建生信Python分析环境
# ============================================================

# 创建新的conda环境
conda create -n bioinfo_py python=3.10 -y

# 激活环境
conda activate bioinfo_py

# 安装核心数据分析包
conda install -c conda-forge pandas numpy scipy matplotlib seaborn scikit-learn -y

# 安装统计分析包
conda install -c conda-forge statsmodels -y

# 安装生信特定包
conda install -c conda-forge scikit-bio -y   # 距离矩阵、PCoA、PERMANOVA

# 安装Jupyter（交互式分析）
conda install -c conda-forge jupyterlab -y

# 验证安装
python -c "
import pandas as pd
import numpy as np
import scipy
import matplotlib
import seaborn as sns
import sklearn
print('所有包安装成功！')
print(f'pandas: {pd.__version__}')
print(f'numpy: {np.__version__}')
print(f'scipy: {scipy.__version__}')
"

---

快速复习卡片

问题	一句话答案
pandas读TSV	pd.read_csv("file.tsv", sep="\t", index_col=0)
计算相对丰度	df.apply(lambda col: col / col.sum(), axis=0)
DataFrame vs Series	DataFrame是二维表格，Series是一维列
Shannon指数公式	H' = -sum(pi * ln(pi))
Wilcoxon检验函数	scipy.stats.mannwhitneyu()
多重检验校正	multipletests(p_values, method='fdr_bh')
PCA vs PCoA	PCA用原始矩阵+欧氏距离，PCoA用距离矩阵+任意距离
loc vs iloc	loc按标签名，iloc按整数位置
axis=0 vs axis=1	axis=0沿行方向（结果按列），axis=1沿列方向（结果按行）
深拷贝	df.copy() 创建独立副本，避免修改原数据

---

学习路径建议

1. 先掌握：pandas 读写/筛选/分组 + matplotlib 基本图表（日常工作最常用）
2. 再学习：scipy 统计检验 + seaborn 高级可视化（差异分析和发表级图表）
3. 进阶：scikit-learn PCA/随机森林 + scikit-bio PCoA/PERMANOVA
4. 面试重点：能手写 Shannon 指数计算 + 解释 PCA vs PCoA 的区别