Python数据分析（Python for Bioinformatics）¶

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

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)}个")

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))

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")

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")

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()

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()

# ============================================================
# 这些写法在新版本中已移除或废弃，面试时要知道！
# ============================================================

# ---------- pd.append() 已在 pandas 2.0 移除 ----------
# 错误（pandas 2.0+ 会报 AttributeError）：
# df1 = df1.append(df2)

# 正确：改用 pd.concat()
df1 = pd.concat([df1, df2], ignore_index=True)

# ---------- np.int / np.float 已在 NumPy 1.24+ 移除 ----------
# 错误（NumPy 1.24+ 会报 AttributeError）：
# x = np.int(5)
# y = np.float(3.14)

# 正确：使用 Python 内置类型或明确的 NumPy 类型
x = int(5)          # 或 np.int64(5)
y = float(3.14)     # 或 np.float64(3.14)

# ---------- df.iteritems() 已在 pandas 2.0 移除 ----------
# 错误：
# for col_name, col_data in df.iteritems(): ...

# 正确：改用 df.items()
for col_name, col_data in df.items():
    pass

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 = 简单拼接（堆叠）

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 按行除

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

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分析环境
# ============================================================

# 创建新的conda环境（scikit-bio 0.6+ 已支持 Python 3.12）
conda create -n bioinfo_py python=3.12 -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__}')
"