生物信息学入门指南：如何用Python分析基因表达数据（附实战代码）

如果你是一名程序员或数据分析师，对生命科学充满好奇，看着实验室里那些神秘的基因数据却不知从何下手，那么这篇文章就是为你准备的。生物信息学早已不是生物学家的专属领域，它正成为连接代码与生命的桥梁。想象一下，用你熟悉的Python，就能解读细胞内部的“对话”，从海量的A、T、C、G字符中，发现疾病的关键或生命的规律。这听起来很酷，对吧？本指南将完全从零开始，抛开复杂的生物学理论，聚焦于最核心的基因表达数据分析流程。我们将使用真实的、可获取的数据集，手把手带你完成从数据下载、清洗、探索到可视化和基础统计推断的全过程。无论你是想转行进入生物科技公司，还是仅仅希望拓展自己的数据科学技能树，这篇实战指南都将为你提供一个坚实、可操作的起点。

1. 环境准备与数据获取：搭建你的生信工作台

工欲善其事，必先利其器。在开始分析之前，我们需要一个稳定、高效的Python数据分析环境。对于生物信息学而言，Anaconda 是一个近乎完美的选择，它集成了Python、R以及众多科学计算和生物信息学专用的包，并能轻松管理环境依赖，避免版本冲突。

1.1 创建专属的Conda环境

我强烈建议为生物信息学项目创建一个独立的环境。打开你的终端（Windows用户请使用Anaconda Prompt），执行以下命令：

# 创建一个名为 bioinfo 的新环境，并指定Python版本
conda create -n bioinfo python=3.9

# 激活该环境
conda activate bioinfo

接下来，安装我们核心的分析库。pandas 和 numpy 是数据处理的基础，scipy 和 statsmodels 用于统计检验，scikit-learn 可能会用于后续的机器学习分析，而 matplotlib 和 seaborn 则是可视化的利器。

# 安装核心数据分析库
conda install pandas numpy scipy statsmodels scikit-learn matplotlib seaborn jupyter

对于生物信息学特有的数据格式（如.bam, .vcf, .fastq等），虽然本文不涉及底层序列处理，但了解 Biopython 这个强大的工具包很有必要，你可以选择安装：

conda install -c conda-forge biopython

1.2 获取实战数据集：基因表达矩阵

理论再好，没有数据也是空谈。为了本次实战，我们将使用一个经典的公共数据集——TCGA（癌症基因组图谱） 中乳腺癌样本的RNA-Seq基因表达数据的一个子集。为了方便初学者，我已经将处理好的数据上传到云端。

注意：在实际科研中，原始测序数据（FASTQ文件）需要经过复杂的比对、定量流程才能得到基因表达矩阵。本指南为了聚焦于分析，直接使用已定量好的“表达矩阵”作为起点。

我们可以用 pandas 直接从网络地址读取数据。这个数据集的行是基因，列是样本，单元格中的数值代表该基因在该样本中的表达水平（通常是经过标准化处理的读数，如FPKM或TPM）。

import pandas as pd
import numpy as np

# 示例数据URL（此处为示意，请替换为真实可用的数据链接）
# 在实际操作中，你可以从 GEO (Gene Expression Omnibus) 或 TCGA 官网下载
# 这里我们假设有一个CSV文件
data_url = "https://example.com/breast_cancer_expression_subset.csv" # 请替换为真实链接

# 为了演示，我们创建一个模拟的小型数据框
np.random.seed(42) # 确保结果可重复
gene_names = [f'Gene_{i}' for i in range(1, 101)] # 100个基因
sample_names = [f'Sample_Normal_{i}' for i in range(1, 6)] + [f'Sample_Tumor_{i}' for i in range(1, 6)] # 5正常+5肿瘤

# 模拟数据：正常样本表达量较低，肿瘤样本某些基因高表达
data = np.random.randn(100, 10) * 5 + 20 # 基础表达水平
# 让前20个基因在肿瘤样本中表达上调
data[:20, 5:] += np.random.randn(20, 5) * 3 + 15

expression_df = pd.DataFrame(data, index=gene_names, columns=sample_names)
print("基因表达数据框维度：", expression_df.shape)
print(expression_df.head()) # 查看前5行

运行上述代码，你会看到一个100行（基因）、10列（样本）的DataFrame。这就是我们后续所有分析的基石。

2. 数据清洗与预处理：从原始矩阵到可靠数据

拿到原始表达矩阵后，直接进行分析往往会得到误导性的结果。数据清洗是确保分析质量的关键一步，主要包括处理缺失值、过滤低表达基因和进行数据标准化。

2.1 审视数据质量与处理缺失值

首先，检查数据中是否存在缺失值（NaN）。

# 检查缺失值
missing_values = expression_df.isnull().sum().sum()
print(f"数据集中缺失值的总数量：{missing_values}")

if missing_values > 0:
    # 策略1：删除包含缺失值的基因（行）
    # expression_df_cleaned = expression_df.dropna(axis=0)
    
    # 策略2：用中位数或均值填充（需谨慎，在生信中不常用）
    # 更常见的做法是，如果缺失过多，直接删除该基因
    # 这里演示填充（仅当缺失很少时）
    from sklearn.impute import SimpleImputer
    imputer = SimpleImputer(strategy='median')
    expression_df_filled = pd.DataFrame(imputer.fit_transform(expression_df),
                                         columns=expression_df.columns,
                                         index=expression_df.index)
    print("已用中位数填充缺失值。")
else:
    expression_df_filled = expression_df.copy()
    print("数据完整，无缺失值。")

在基因表达分析中，一个基因如果在大多数样本中都不表达或表达量极低，它携带的信息量就很少，且会增加统计检验的负担。因此，我们需要过滤掉这些低表达基因。

2.2 过滤低表达基因

一个常用的过滤标准是：保留在至少X%的样本中，表达量高于某个阈值（例如，大于1 TPM或FPKM）的基因。这里我们设定：保留在至少20%的样本中表达量大于1的基因。

# 设定阈值
min_expression = 1
min_samples = int(expression_df_filled.shape[1] * 0.20) # 20%的样本

# 计算每个基因在多少个样本中表达量大于阈值
gene_pass_filter = (expression_df_filled > min_expression).sum(axis=1) >= min_samples

# 应用过滤
expression_df_filtered = expression_df_filled[gene_pass_filter]
print(f"过滤前基因数：{expression_df_filled.shape[0]}")
print(f"过滤后基因数：{expression_df_filtered.shape[0]}")

2.3 数据标准化与转换

不同样本之间的测序深度（总读数）可能不同，不同基因的表达量范围也可能差异巨大。为了进行公平的比较，我们需要进行标准化。对于样本间的标准化，通常在上游定量步骤（如DESeq2, edgeR）已完成。我们这里关注的是基因层面的标准化（Z-score标准化），这在后续的热图聚类等分析中很常见。

from sklearn.preprocessing import StandardScaler

# 对基因进行Z-score标准化（按行标准化）
# 即：对于每个基因，在所有样本中的表达值转换为均值为0，标准差为1。
scaler = StandardScaler()
expression_df_zscore = pd.DataFrame(scaler.fit_transform(expression_df_filtered.T).T, # 注意转置
                                     columns=expression_df_filtered.columns,
                                     index=expression_df_filtered.index)

print("Z-score标准化后的数据（前5个基因）：")
print(expression_df_zscore.iloc[:5, :5])

此外，基因表达数据通常呈现右偏态分布，对其进行对数转换可以使数据更接近正态分布，满足许多统计方法的前提假设。

# 对数转换（通常使用 log2(expression + 1)），+1是为了避免对0取对数
expression_df_log = np.log2(expression_df_filtered + 1)
print("对数转换后的数据范围：")
print(f"最小值：{expression_df_log.values.min():.2f}, 最大值：{expression_df_log.values.max():.2f}")

经过以上步骤，我们得到了一个干净、经过初步处理的基因表达矩阵 expression_df_filtered（原始尺度）和 expression_df_log（对数尺度），可以用于下游分析了。

3. 探索性分析与可视化：看见数据背后的故事

在开始复杂的统计检验之前，通过可视化探索数据是至关重要的。它能帮助我们理解数据的整体结构，检查批次效应，观察样本分组情况。

3.1 样本间相关性分析与热图

我们可以计算样本两两之间的表达谱相关性。高相关性意味着两个样本的基因表达模式相似。

import seaborn as sns
import matplotlib.pyplot as plt

# 计算样本间相关系数矩阵（使用对数转换后的数据更稳定）
sample_correlation = expression_df_log.corr() # 默认是皮尔逊相关系数

# 绘制热图
plt.figure(figsize=(10, 8))
sns.heatmap(sample_correlation,
            annot=True, # 在方格中显示数值
            fmt='.2f', # 数值格式
            cmap='coolwarm', # 颜色映射
            center=0.8, # 颜色中心点
            square=True)
plt.title('样本间基因表达相关性热图')
plt.tight_layout()
plt.show()

从热图中，你可能会发现正常样本（Sample_Normal_*）彼此相关性高，肿瘤样本（Sample_Tumor_*）彼此相关性高，而正常与肿瘤样本间的相关性较低。这初步提示我们的样本分组是合理的。

3.2 主成分分析：降维与批次效应检测

主成分分析（PCA）是一种强大的降维技术，能将成千上万个基因的信息压缩到几个主成分（PC）中，从而在二维或三维空间中直观展示样本间的整体差异。

from sklearn.decomposition import PCA

# 准备数据：转置，使得行是样本，列是基因
X = expression_df_log.T.values # 样本 x 基因
sample_labels = expression_df_log.T.index.tolist()
group_labels = ['Normal' if 'Normal' in s else 'Tumor' for s in sample_labels] # 创建分组标签

# 执行PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# 将结果转换为DataFrame便于绘图
pca_df = pd.DataFrame(data=X_pca, columns=['PC1', 'PC2'])
pca_df['Sample'] = sample_labels
pca_df['Group'] = group_labels

# 绘制PCA散点图
plt.figure(figsize=(10, 8))
scatter = plt.scatter(pca_df['PC1'], pca_df['PC2'],
                      c=[0 if g=='Normal' else 1 for g in pca_df['Group']],
                      cmap='viridis', s=100, alpha=0.7)
plt.colorbar(scatter, label='Group (0=Normal, 1=Tumor)')
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
plt.title('基因表达数据的PCA分析')

# 为每个点添加样本标签
for i, row in pca_df.iterrows():
    plt.annotate(row['Sample'], (row['PC1'], row['PC2']), xytext=(5, 5), textcoords='offset points', fontsize=9)

plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

理想的PCA图中，同一组内的样本应该聚集在一起，不同组间的样本应该明显分离。如果样本按实验批次（而非生物学分组）聚集，则说明存在强烈的批次效应，需要在后续分析中校正。

3.3 单个基因的表达分布：小提琴图与箱线图

假设我们关注一个在肿瘤中可能高表达的基因 Gene_1，我们可以直观地比较它在两组间的表达分布。

# 提取 Gene_1 的表达数据
gene_of_interest = 'Gene_1'
if gene_of_interest in expression_df_log.index:
    gene_data = expression_df_log.loc[gene_of_interest]
    # 转换为适合绘图的长格式DataFrame
    plot_df = pd.DataFrame({
        'Expression': gene_data.values,
        'Sample': gene_data.index,
        'Group': ['Normal' if 'Normal' in s else 'Tumor' for s in gene_data.index]
    })
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 小提琴图
    sns.violinplot(x='Group', y='Expression', data=plot_df, ax=axes[0], inner='box', palette='Set2')
    axes[0].set_title(f'{gene_of_interest} 表达分布 - 小提琴图')
    axes[0].set_ylabel('log2(Expression+1)')
    
    # 箱线图叠加散点图（显示个体样本）
    sns.boxplot(x='Group', y='Expression', data=plot_df, ax=axes[1], width=0.5, palette='Set2', fliersize=0)
    sns.stripplot(x='Group', y='Expression', data=plot_df, ax=axes[1], color='black', alpha=0.7, jitter=True)
    axes[1].set_title(f'{gene_of_interest} 表达分布 - 箱线图')
    axes[1].set_ylabel('log2(Expression+1)')
    
    plt.tight_layout()
    plt.show()
else:
    print(f"基因 {gene_of_interest} 不在过滤后的数据集中。")

这种可视化能清晰展示基因表达的中位数、四分位范围以及整体分布形状，是组间比较的经典方法。

4. 差异表达分析：寻找关键基因

这是生物信息学核心中的核心：找出在不同条件（如正常 vs. 肿瘤）下表达水平有显著差异的基因。这些差异表达基因（DEGs）往往是理解生物学过程或疾病机制的关键。

4.1 使用T检验进行组间比较

对于两组比较，独立样本T检验是一个基础方法。我们将对过滤后的每一个基因，在正常组和肿瘤组之间进行T检验。

from scipy import stats

# 准备分组数据
normal_samples = [col for col in expression_df_log.columns if 'Normal' in col]
tumor_samples = [col for col in expression_df_log.columns if 'Tumor' in col]

deg_results = []
for gene in expression_df_log.index:
    normal_expr = expression_df_log.loc[gene, normal_samples].values
    tumor_expr = expression_df_log.loc[gene, tumor_samples].values
    
    # 执行独立双样本T检验（默认假设方差不相等，使用Welch's t-test）
    t_stat, p_value = stats.ttest_ind(tumor_expr, normal_expr, equal_var=False)
    
    # 计算效应量：Cohen's d
    pooled_std = np.sqrt(((len(normal_expr)-1)*np.var(normal_expr) + (len(tumor_expr)-1)*np.var(tumor_expr)) / (len(normal_expr)+len(tumor_expr)-2))
    cohens_d = (np.mean(tumor_expr) - np.mean(normal_expr)) / pooled_std if pooled_std != 0 else 0
    
    # 计算log2 Fold Change (肿瘤/正常)
    log2fc = np.mean(tumor_expr) - np.mean(normal_expr)
    
    deg_results.append({
        'Gene': gene,
        'log2FC': log2fc,
        'T_statistic': t_stat,
        'P_value': p_value,
        'Cohen_d': cohens_d
    })

# 转换为DataFrame
deg_df = pd.DataFrame(deg_results)
deg_df.set_index('Gene', inplace=True)

# 查看结果
print("差异表达分析结果（前10个基因）：")
print(deg_df.sort_values('P_value').head(10))

4.2 多重检验校正与显著性判断

对成千上万个基因同时进行检验，会极大增加假阳性（I类错误）的风险。因此，必须对P值进行多重检验校正。最常用的方法是 错误发现率（FDR） 校正，即计算q值。

from statsmodels.stats.multitest import multipletests

# 进行FDR校正（Benjamini-Hochberg方法）
reject, pvals_corrected, _, _ = multipletests(deg_df['P_value'], alpha=0.05, method='fdr_bh')
deg_df['P_adj'] = pvals_corrected

# 定义差异表达基因的标准：通常为 |log2FC| > 1 且 P_adj < 0.05
deg_df['Significant'] = (deg_df['P_adj'] < 0.05) & (np.abs(deg_df['log2FC']) > 1)

# 统计显著基因数量
num_sig = deg_df['Significant'].sum()
print(f"在FDR<0.05且|log2FC|>1的标准下，共发现 {num_sig} 个差异表达基因。")

# 查看最显著的上调和下调基因
sig_up = deg_df[(deg_df['Significant']) & (deg_df['log2FC'] > 1)].sort_values('log2FC', ascending=False)
sig_down = deg_df[(deg_df['Significant']) & (deg_df['log2FC'] < -1)].sort_values('log2FC', ascending=True)

print(f"\n表达上调最显著的5个基因：")
print(sig_up[['log2FC', 'P_adj']].head())
print(f"\n表达下调最显著的5个基因：")
print(sig_down[['log2FC', 'P_adj']].head())

4.3 结果可视化：火山图与热图

火山图是展示差异表达分析结果的经典方式，它能同时显示基因表达变化幅度（log2FC）和统计显著性（-log10(P值)）。

plt.figure(figsize=(10, 8))
# 绘制所有基因
plt.scatter(deg_df['log2FC'], -np.log10(deg_df['P_adj']),
            c=['red' if sig else 'gray' for sig in deg_df['Significant']],
            alpha=0.6, s=20, edgecolors='none')

# 添加阈值线
plt.axvline(x=1, color='blue', linestyle='--', linewidth=0.8, alpha=0.7)
plt.axvline(x=-1, color='blue', linestyle='--', linewidth=0.8, alpha=0.7)
plt.axhline(y=-np.log10(0.05), color='green', linestyle='--', linewidth=0.8, alpha=0.7)

plt.xlabel('log2 Fold Change (Tumor / Normal)')
plt.ylabel('-log10(Adjusted P-value)')
plt.title('差异表达基因火山图')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

接下来，我们可以提取显著差异表达基因的表达矩阵，绘制热图来直观展示其在不同样本中的表达模式。

# 提取显著基因的表达数据（Z-score标准化后的数据，便于行间比较）
sig_genes = deg_df[deg_df['Significant']].index.tolist()
if len(sig_genes) > 0:
    sig_expression_data = expression_df_zscore.loc[sig_genes]
    
    # 绘制聚类热图
    plt.figure(figsize=(12, 16))
    # 对行（基因）和列（样本）进行聚类
    g = sns.clustermap(sig_expression_data,
                       cmap='RdBu_r', # 红蓝渐变色，中心为白色
                       center=0,
                       figsize=(12, 16),
                       yticklabels=True if len(sig_genes) < 50 else False, # 基因太多时不显示标签
                       col_cluster=True,
                       row_cluster=True,
                       dendrogram_ratio=(.1, .2),
                       cbar_pos=(.02, .32, .03, .2))
    g.ax_heatmap.set_xlabel('Samples')
    g.ax_heatmap.set_ylabel('Significant DEGs')
    g.ax_heatmap.set_title('差异表达基因聚类热图 (Z-score by Gene)')
    plt.show()
else:
    print("未发现显著差异表达基因，无法绘制热图。")

热图中，红色代表高表达，蓝色代表低表达。你可以看到基因是如何被聚类的（功能可能相关），以及样本是如何被聚类的（反映了生物学分组）。

5. 功能富集分析初探：理解基因列表的生物学意义

找到一堆差异表达基因后，下一个关键问题是：这些基因共同参与了哪些生物学过程？这需要通过功能富集分析来实现。虽然完整的富集分析通常依赖专门的数据库（如GO, KEGG）和工具（如clusterProfiler in R），但我们可以在Python中完成一些基础步骤和可视化。

5.1 准备基因列表与背景集

假设我们有一个感兴趣的基因列表（例如，从差异分析中得到的前50个上调基因），我们想了解它们是否在某些通路中富集。

# 获取前50个上调基因
top_up_genes = sig_up.head(50).index.tolist()
print(f"用于富集分析的上调基因数量：{len(top_up_genes)}")
print(top_up_genes[:10]) # 打印前10个看看

在实际分析中，你需要将这些基因标识符（如Gene_1）映射到标准的基因符号（如TP53, BRCA1），并使用专门的包（如 gseapy）连接到富集分析数据库。这里，我们模拟一个简单的超几何检验过程，并展示如何用Python进行基础计算。

5.2 模拟富集分析结果与可视化

假设我们通过某个在线工具或gseapy包得到了富集分析结果，数据格式通常包含通路名称、富集分数、P值、包含的基因等。

# 模拟富集分析结果（在实际中，这是从分析工具得到的）
enrichment_data = {
    'Pathway': ['Cell Cycle', 'DNA Repair', 'Immune Response', 'Metabolic Process', 'Apoptosis'],
    'P_value': [1.2e-10, 3.5e-7, 0.00015, 0.0023, 0.018],
    'Gene_Count': [28, 19, 35, 42, 12], # 该通路中属于我们基因列表的基因数
    'Total_Genes': [150, 120, 300, 500, 100] # 该通路总基因数
}
enrich_df = pd.DataFrame(enrichment_data)
enrich_df['-log10(P)'] = -np.log10(enrich_df['P_value'])
enrich_df['Gene_Ratio'] = enrich_df['Gene_Count'] / enrich_df['Total_Genes']

# 绘制条形图
plt.figure(figsize=(10, 6))
bars = plt.barh(enrich_df['Pathway'], enrich_df['-log10(P)'], color='skyblue')
plt.xlabel('-log10(P-value)')
plt.title('Top Enriched Pathways for Up-regulated Genes')
# 在条形末端添加基因比例
for i, (ratio, pval) in enumerate(zip(enrich_df['Gene_Ratio'], enrich_df['-log10(P)'])):
    plt.text(pval + 0.05, i, f'{ratio:.2%}', va='center')
plt.tight_layout()
plt.show()

这个简单的条形图展示了哪些通路在基因列表中显著富集，以及富集的程度。在实际项目中，你会使用更专业的工具得到包含数百个通路的结果，并需要对其进行筛选和深入解读。

5.3 后续分析方向建议

完成基础分析后，你可以沿着多个方向深入：

  • 深入挖掘：对关键基因进行生存分析（结合临床数据），或构建共表达网络。
  • 工具进阶：学习使用 DESeq2（R）或 limma（R）进行更稳健的RNA-Seq差异分析，它们专门处理计数数据并考虑离散度。
  • 流程整合：将分析流程脚本化，例如使用 Snakemake 或 Nextflow 创建可重复的分析管道。
  • 数据库利用：熟练使用 GEO、TCGA、GTEx 等公共数据库获取数据。

生物信息学分析就像侦探破案，数据是你的线索，统计和编程是你的工具，而生物学知识则是你的地图。从这份指南起步，保持好奇心，多动手实践，你很快就能从基因数据的海洋中，打捞出属于自己的科学发现。

Logo
DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

cover

巅峰赛上分秘籍：细节决定成败

avatar DAMO开发者矩阵
cover

医疗AI的“破局时刻”：从辅助诊断到手术机器人，真实落地的3大困局

avatar DAMO开发者矩阵
cover

如何用python做一个机器人

avatar DAMO开发者矩阵