介绍～

本文将从最基础的医学生信分析讲起，带领大家一步步入门，最终完成并发表属于自己的第一篇SCI论文。全文内容基于个人经验总结，特别适合零基础的小伙伴学习。一般而言，文章设计需要一个数据集进行初步分析（测试集），再用另一个独立数据集进行验证（验证集）。这里我们主要选取其中一个数据集进行详细讲解。我们以中国脑胶质瘤基因图谱计划（CGGA）中的公开数据集CGGA-325为例，介绍数据预处理与基因差异分析的基本流程。该数据集可通过CGGA官网（https://www.cgga.org.cn）获取。接下来，我们将重点讲解该部分的数据处理与分析方法。

1. 数据的读取与数据合并

下载的数据分为临床数据集与测序数据集，不要对原始数据进行改动，我们直接在pandas中进行操作。首先，导入CGGA_325_clinical.csv临床样本数据集：

import pandas as pd

# 读取训练集 - 使用文件自带的列名，然后选择需要的列
dataset_1 = pd.read_csv("CGGA_325_clinical.csv", 
                       header=0,         # 使用第一行作为列名
                       index_col=None,   # 第一列作为行索引
                       sep=',\s*',       # 修正正则表达式
                       na_values=["?"],  # 将"?"识别为缺失值(NaN)
                       engine='python')  # 使用Python引擎以支持正则分隔符


# 显示前10行
dataset_1

其次，以同样的方式导入CGGA_325_gene.csv基因组数据：

dataset_2 = pd.read_csv("CGGA_325_gene.csv", 
                        header=0,           # 第一行作为列名
                        index_col=0,        # 第一列作为行索引
                        sep=',\s*',         # 修正正则表达式
                        na_values=["?"],    # 将"?"识别为缺失值(NaN)
                        engine='python'     # 使用Python引擎以支持正则分隔符
)

# 显示前10行
dataset_2

最后，将基因组数据集进行转置操作，并对临床数据集与测序数据集进行合并，且命名为CGGA_325_cinical_gene.csv，合并后的数据集结果如下所示：

2. Log₂FC、P Value计算

2.1 导入整理好的新数据集

在这一步导入前面整理好的合并数据集CGGA_325_cinical_gene.csv，以后所有的操作均在这个数据集中完成：

import pandas as pd
CGGA_325_cinical_gene = pd.read_csv("CGGA_325_cinical_gene.csv",
                       header=0,         # 使用第一行作为列名
                       index_col=0,      # 第一列作为行索引
                       sep=',\s*',       # 修正正则表达式
                       na_values=["?"],  # 将"?"识别为缺失值(NaN)
                       engine='python')  # 使用Python引擎以支持正则分隔符


# 显示前10行
CGGA_325_cinical_gene

2.2 选择分组

这一步我们需要分组的选择，本节我们以高低级别胶质瘤进行分组，病理组织与WHO分级见补充部分

CGGA_325_cinical_gene['Histology'].unique()
结果：array(['GBM', 'AA', 'rGBM', 'rA', 'O', 'rAA', 'sGBM', 'A', nan, 'AO',
       'rO'], dtype=object)

补充：

高级别胶质瘤（High-Grade Glioma） 这些类别对应WHO III级或IV级，肿瘤恶性程度较高，生长较快：

  • GBM: 胶质母细胞瘤（Glioblastoma），WHO IV级。
  • AA: 间变性星形细胞瘤（Anaplastic Astrocytoma），WHO III级。
  • rGBM: 复发性胶质母细胞瘤（Recurrent GBM），基于GBM，WHO IV级。
  • rAA: 复发性间变性星形细胞瘤（Recurrent AA），基于AA，WHO III级。
  • sGBM: 继发性胶质母细胞瘤（Secondary GBM），通常从低级别进展而来，但最终为WHO IV级。
  • AO: 间变性少突胶质细胞瘤（Anaplastic Oligodendroglioma），WHO III级。

低级别胶质瘤（Low-Grade Glioma） 这些类别对应WHO I级或II级，肿瘤恶性程度较低，生长较慢：

  • rA: 复发性星形细胞瘤（Recurrent Astrocytoma），基于A（星形细胞瘤），通常WHO II级。
  • O: 少突胶质细胞瘤（Oligodendroglioma），通常WHO II级。
  • A: 星形细胞瘤（Astrocytoma），通常WHO II级。
  • rO: 复发性少突胶质细胞瘤（Recurrent Oligodendroglioma），基于O，通常WHO II级。

2.3 数据编码

# 创建真实副本
df = CGGA_325_cinical_gene.copy()

# 明确指定需要编码的列
categorical_columns = ['Censor','Gender','Grade', 'Histology', 'IDH', 'MGMTp', '1p19q']

# 创建字典保存映射关系
encoding_mappings = {}

for col in categorical_columns:
    # 进行编码并保存映射关系
    encoded, categories = df[col].factorize()
    df[col] = encoded
    encoding_mappings[col] = dict(enumerate(categories))
    
# 查看所有列的编码映射
for col, mapping in encoding_mappings.items():
    print(f"\n{col} 列的编码映射:")
    for code, category in mapping.items():
        print(f"  {code}: {category}")
Censor 列的编码映射:
  0: Alive
  1: Dead

Gender 列的编码映射:
  0: Male
  1: Female

Grade 列的编码映射:
  0: WHO IV
  1: WHO III
  2: WHO II

Histology 列的编码映射:
  0: GBM
  1: AA
  2: rGBM
  3: rA
  4: O
  5: rAA
  6: sGBM
  7: A
  8: AO
  9: rO

IDH 列的编码映射:
  0: Wildtype
  1: Mutant

MGMTp 列的编码映射:
  0: Un-methylated
  1: Methylated

1p19q 列的编码映射:
  0: Non-codel
  1: Codel

2.4 Log₂FC、P Value计算

本节以高低级别胶质瘤进行分组，WHO IV与WHO III期属高级别，WHO I与WHO II属低级别。

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

# 提取基因列（从A1BG开始）
gene_columns = [col for col in df.columns if col not in [
    'Histology', 'Grade', 'Gender', 'Age', 'OS',
    'Censor', 'IDH', '1p19q', 'MGMTp'
]]

print(f"找到 {len(gene_columns)} 个基因进行差异表达分析")

# 根据Grade分组
group1 = df[df['Grade'] == 2]  # 低级别
group2 = df[df['Grade'].isin([0, 1])]  # 高级别

print(f"组1 (低级别) 样本数: {len(group1)}")
print(f"组2 (高级别) 样本数: {len(group2)}")

# 计算每个基因的logFC和p值
results = []

for gene in gene_columns:
    # 获取两组的表达值
    expr1 = group1[gene].values
    expr2 = group2[gene].values

    # 计算均值
    mean1 = np.mean(expr1)
    mean2 = np.mean(expr2)

    # 计算logFC (以2为底)
    if mean1 == 0 or mean2 == 0:
        # 处理零值情况
        logFC = np.nan
    else:
        logFC = np.log2(mean2 / mean1)

    # 进行t检验
    if len(expr1) > 1 and len(expr2) > 1:  # 确保每组至少有2个样本
        t_stat, p_value = stats.ttest_ind(expr1, expr2, equal_var=False)  # Welch's t-test
    else:
        p_value = np.nan

    results.append({
        'Gene': gene,
        'Mean_Group1': mean1,
        'Mean_Group2': mean2,
        'logFC': logFC,
        'p_value': p_value
    })

# 创建结果DataFrame
result_df = pd.DataFrame(results)

# 多重检验校正 (FDR)
result_df['p_adjust'] = multipletests(result_df['p_value'], method='fdr_bh')[1]

# 排序结果（按p值）
result_df = result_df.sort_values('p_value')

# 显示前20个最显著的基因
print("\n前20个最显著的差异表达基因:")
print(result_df.head(20))

# 保存结果到CSV文件
result_df.to_csv('differential_expression_results.csv', index=False)
print("\n结果已保存到 'differential_expression_results.csv'")

我们得到differential_expression_results.csv数据集，结果如下：

3. 火山图绘制

3.1 导入分析结果数据集

from adjustText import adjust_text
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

LogFC_Pvalue = pd.read_csv("differential_expression_results.csv",
                       header=0,         # 使用第一行作为列名
                       index_col=None,      # 第一列作为行索引
                       sep=',\s*',       # 修正正则表达式
                       na_values=["?"],  # 将"?"识别为缺失值(NaN)
                       engine='python')  # 使用Python引擎以支持正则分隔符


result_df = LogFC_Pvalue

3.2 火山图绘制

import matplotlib.pyplot as plt
import numpy as np
from adjustText import adjust_text  # 确保已安装：pip install adjustText

def plot_volcano(result_df, logFC_threshold=1, p_value_threshold=0.05,
                 max_labels=50, figsize=(12, 8), title='Volcano Plot',
                 custom_genes=None):
    """
    绘制SCI风格火山图

    参数:
    result_df: 包含logFC和p_value的数据框
    logFC_threshold: logFC阈值，默认1（|logFC| > 1为显著）
    p_value_threshold: p值阈值，默认0.05
    max_labels: 最大标注基因数量，默认50
    figsize: 图形大小
    title: 图表标题
    custom_genes: 自定义要标注的基因列表，默认为None
    """
    # 设置SCI风格：无衬线字体、白色背景、网格
    plt.style.use('default')
    plt.rcParams['font.family'] = 'Arial'  # 使用Arial字体
    plt.rcParams['axes.edgecolor'] = 'black'
    plt.rcParams['axes.linewidth'] = 0.8
    plt.rcParams['grid.color'] = '#E6E6E6'
    plt.rcParams['grid.linestyle'] = '--'

    # 创建图形
    fig, ax = plt.subplots(figsize=figsize)
    ax.set_facecolor('white')  # 设置背景为白色

    # 计算-log10(p_value)
    result_df['neg_log10_p'] = -np.log10(result_df['p_value'])

    # 识别上调、下调和不变的基因
    up_regulated = (result_df['logFC'] > logFC_threshold) & (result_df['p_value'] < p_value_threshold)
    down_regulated = (result_df['logFC'] < -logFC_threshold) & (result_df['p_value'] < p_value_threshold)
    not_significant = ~(up_regulated | down_regulated)

    # 绘制点：使用SCI风格配色和透明度
    ax.scatter(result_df['logFC'][not_significant],
               result_df['neg_log10_p'][not_significant],
               c='#B0B0B0', alpha=0.3, s=40, label='Not significant', edgecolors='none')

    ax.scatter(result_df['logFC'][up_regulated],
               result_df['neg_log10_p'][up_regulated],
               c='#D62728', alpha=0.4, s=70, label=f'Up (logFC > {logFC_threshold})', edgecolors='black')

    ax.scatter(result_df['logFC'][down_regulated],
               result_df['neg_log10_p'][down_regulated],
               c='#1F77B4', alpha=0.4, s=70, label=f'Down (logFC < -{logFC_threshold})', edgecolors='black')

    # 添加阈值线：黑色虚线，更细
    ax.axhline(-np.log10(p_value_threshold), color='black', linestyle='--', linewidth=0.8, alpha=0.7)
    ax.axvline(logFC_threshold, color='black', linestyle='--', linewidth=0.8, alpha=0.7)
    ax.axvline(-logFC_threshold, color='black', linestyle='--', linewidth=0.8, alpha=0.7)
    ax.axvline(0, color='black', linestyle='-', linewidth=0.5, alpha=0.5)

    # 标注自定义基因
    if 'Gene' in result_df.columns and custom_genes is not None:
        custom_genes_df = result_df[result_df['Gene'].isin(custom_genes)].copy()
        if not custom_genes_df.empty:
            custom_genes_df['regulation'] = custom_genes_df['logFC'].apply(
                lambda x: 'up' if x > 0 else 'down'
            )
            max_labels = min(max_labels, len(custom_genes_df))
            texts = []
            for i, (idx, row) in enumerate(custom_genes_df.head(max_labels).iterrows()):
                # 根据上下调选择颜色
                if row['regulation'] == 'up':
                    bbox_color = '#FFE4E1'  # 浅红色背景
                    text_color = '#8B0000'  # 深红色文字
                else:
                    bbox_color = '#E6F3FF'  # 浅蓝色背景
                    text_color = '#00008B'  # 深蓝色文字

                text = ax.annotate(row['Gene'],
                                  (row['logFC'], row['neg_log10_p']),
                                  fontsize=9,
                                  fontweight='bold',
                                  color=text_color,
                                  bbox=dict(boxstyle="round,pad=0.3",
                                           facecolor=bbox_color,
                                           edgecolor=text_color,
                                           alpha=0.8,
                                           linewidth=0.5))
                texts.append(text)

            # 自动调整文本位置，避免重叠
            if texts:
                adjust_text(texts,
                            ax=ax,
                            arrowprops=dict(arrowstyle="->", color='gray', lw=0.5, alpha=0.6),
                            expand_text=(1.2, 1.2),
                            force_text=0.5,
                            ensure_inside_axes=True)

    # 设置标签和标题：使用粗体和大字体
    ax.set_xlabel('log2(Fold Change)', fontsize=12, fontweight='bold')
    ax.set_ylabel('-log10(p-value)', fontsize=12, fontweight='bold')
    ax.set_title(title, fontsize=14, fontweight='bold', pad=20)

    # 添加网格：浅灰色虚线
    ax.grid(True, alpha=0.4, linestyle='--')

    # 图例：位置在右上角，边框透明
    ax.legend(loc='upper right', frameon=True, framealpha=0.9, edgecolor='black')

    # 调整布局
    plt.tight_layout()

    return fig, up_regulated, down_regulated

# 使用示例（保持不变）
my_custom_genes = ['HSPA7']  # 替换为您的基因列表
plt_fig, up_regulated, down_regulated = plot_volcano(result_df,
                   logFC_threshold=2.0,
                   p_value_threshold=0.05,
                   max_labels=10,
                   title='CGGA-325 Volcano Plot: Low-Grade vs High-Grade',
                   custom_genes=my_custom_genes)

# 保存图片（推荐PNG或PDF格式以保持质量）
plt_fig.savefig('volcano_plot_sci_style.png', dpi=300, bbox_inches='tight')
plt_fig.show()

# 输出统计信息
print(f"总基因数: {len(result_df)}")
print(f"显著上调基因数 (logFC > 2.0, p < 0.05): {up_regulated.sum()}")
print(f"显著下调基因数 (logFC < -2.0, p < 0.05): {down_regulated.sum()}")

绘图示例

def plot_volcano(result_df, logFC_threshold=1, p_value_threshold=0.05,
                 max_labels=50, figsize=(12, 8), title='Volcano Plot'):
    """
    绘制火山图

    参数:
    result_df: 包含logFC和p_value的数据框
    logFC_threshold: logFC阈值，默认1（|logFC| > 1为显著）
    p_value_threshold: p值阈值，默认0.05
    max_labels: 最大标注基因数量，默认50
    figsize: 图形大小
    title: 图表标题
    """
    # 创建图形
    plt.figure(figsize=figsize)

    # 计算-log10(p_value)
    result_df['neg_log10_p'] = -np.log10(result_df['p_value'])

    # 识别上调、下调和不变的基因
    up_regulated = (result_df['logFC'] > logFC_threshold) & (result_df['p_value'] < p_value_threshold)
    down_regulated = (result_df['logFC'] < -logFC_threshold) & (result_df['p_value'] < p_value_threshold)
    not_significant = ~(up_regulated | down_regulated)

    # 绘制点
    plt.scatter(result_df['logFC'][not_significant],
                result_df['neg_log10_p'][not_significant],
                c='gray', alpha=0.2, s=50, label='Not significant', edgecolors='None')

    plt.scatter(result_df['logFC'][up_regulated],
                result_df['neg_log10_p'][up_regulated],
                c='red', alpha=0.5, s=70, label=f'Up-regulated (logFC > {logFC_threshold})', edgecolors='black')

    plt.scatter(result_df['logFC'][down_regulated],
                result_df['neg_log10_p'][down_regulated],
                c='blue', alpha=0.5, s=70, label=f'Down-regulated (logFC < -{logFC_threshold})', edgecolors='black')

    # 添加阈值线
    plt.axhline(-np.log10(p_value_threshold), color='black', linestyle='--', alpha=0.8)
    plt.axvline(logFC_threshold, color='black', linestyle='--', alpha=0.8)
    plt.axvline(-logFC_threshold, color='black', linestyle='--', alpha=0.8)
    plt.axvline(0, color='black', linestyle='-', alpha=0.5)

    # 标注显著基因
    if 'Gene' in result_df.columns:
        # 筛选显著基因
        significant_genes = result_df[
            ((result_df['logFC'] > logFC_threshold) | (result_df['logFC'] < -logFC_threshold)) &
            (result_df['p_value'] < p_value_threshold)
        ].copy()

        # 添加一列标识上/下调
        significant_genes['regulation'] = significant_genes['logFC'].apply(
            lambda x: 'up' if x > 0 else 'down'
        )

        # 按显著性排序
        significant_genes = significant_genes.sort_values('p_value')

        # 限制标注数量
        max_labels = min(max_labels, len(significant_genes))
        if max_labels > 0:
            texts = []

            # 标注基因
            for i, (idx, row) in enumerate(significant_genes.head(max_labels).iterrows()):
                # 根据上/下调选择颜色
                if row['regulation'] == 'up':
                    bbox_color = "#FFCCCB"  # 浅红色
                    text_color = "darkred"
                else:
                    bbox_color = "#ADD8E6"  # 浅蓝色
                    text_color = "darkblue"

                text = plt.annotate(row['Gene'],
                                   (row['logFC'], row['neg_log10_p']),
                                   fontsize=8,
                                   alpha=0.9,
                                   color=text_color,
                                   bbox=dict(boxstyle="round,pad=0.2",
                                            facecolor=bbox_color,
                                            edgecolor=text_color,
                                            alpha=0.8))
                texts.append(text)

            # 自动调整文本位置
            if texts:
                adjust_text(texts,
                            arrowprops=dict(arrowstyle="->", color='gray', lw=0.5, alpha=0.7, shrinkA=10, shrinkB=5),
                            expand_text=(1.2, 1.2),
                            expand_points=(1.3, 1.3),
                            force_text=0.5,
                            force_points=0.5,
                            ensure_inside_axes=True)  # 确保文本在坐标轴内

    # 设置标签和标题
    plt.xlabel('log2(Fold Change)')
    plt.ylabel('-log10(p-value)')
    plt.title(title)
    plt.legend()

    # 添加网格
    plt.grid(True, alpha=0.3)

    # 调整布局
    plt.tight_layout()

    return plt, up_regulated, down_regulated

# 使用示例
# 假设你已经有了result_df，包含'logFC', 'p_value', 'Gene'列
# 绘制火山图
plt_fig, up_regulated, down_regulated = plot_volcano(result_df,
                   logFC_threshold=2.5,
                   p_value_threshold=0.05,
                   max_labels=10,
                   title='Volcano Plot: Low-Grade vs High-Grade')

# 保存图片
plt_fig.savefig('volcano_plot.png', dpi=300, bbox_inches='tight')

# 显示图片
plt_fig.show()

# 输出统计信息
print(f"总基因数: {len(result_df)}")
print(f"显著上调基因数 (logFC > 2.5, p < 0.05): {up_regulated.sum()}")
print(f"显著下调基因数 (logFC < -2.5, p < 0.05): {down_regulated.sum()}")
print(f"最显著上调基因: {result_df[up_regulated].nsmallest(154, 'p_value')['Gene'].tolist() if up_regulated.sum() > 0 else '无'}")
print(f"最显著下调基因: {result_df[down_regulated].nsmallest(1, 'p_value')['Gene'].tolist() if down_regulated.sum() > 0 else '无'}")

本节我们将得到高级别与低级别脑胶质瘤的差异基因可视化

# 数据库
Logo
DAMO开发者矩阵

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

更多推荐

cover

ICAI动态认知验证：从任务成功率到认知过程可检验性

avatar DAMO开发者矩阵
cover

游戏引擎是假的，真机器人反而学得更好？Agent 拿「假物理」当裁判，D4RL 分数涨 48%

avatar DAMO开发者矩阵
cover

第302篇 策略梯度——从REINFORCE到现代方法

avatar DAMO开发者矩阵