4. 差异表达数据简介：读取 CSV/TSV 格式

要做差异表达分析，第一步就是把数据拿到手里。你可能会从 RNA-seq 流程拿到一个表格，里面列着基因名、log2FC、p-value、padj 这些字段。格式多半是 CSV 或 TSV。

我个人处理过上百份差异表达结果，踩过不少坑。今天就把读取数据这块讲透。

4.1 差异表达结果长什么样？

先看一个典型的差异表达表格。列名大概是这样：

gene_id,base_mean,log2FC,p_value,padj
ENSG000001,1024.5,2.3,0.0001,0.001
ENSG000002,256.1,-1.8,0.0005,0.008
ENSG000003,512.0,0.4,0.4500,0.780

这是 CSV 格式，逗号分隔。TSV 则是用制表符分隔：

gene_id	base_mean	log2FC	p_value	padj
ENSG000001	1024.5	2.3	0.0001	0.001
ENSG000002	256.1	-1.8	0.0005	0.008

看着差不多，但坑就在分隔符上。你想想看，如果用错了分隔符，pandas 会把整行当成一列，解析全乱。

特性	CSV	TSV
分隔符	逗号 ,	制表符 \t
常见场景	Excel 导出、R 语言	生物信息学 pipe 输出
处理引号	字段可带双引号	通常不带引号
读取基准	pd.read_csv()	pd.read_table() 或 sep='\t'

4.2 用 pandas 读取差异表达表格

我建议统一用 pd.read_csv()，通过 sep 参数指定分隔符。这样代码更一致。

读取 CSV 文件：

import pandas as pd

df = pd.read_csv("deseq2_results.csv")
print(df.head())

读取 TSV 文件：

df = pd.read_csv("deseq2_results.tsv", sep="\t")
print(df.head())

有时候文件名不标准，比如叫 .txt 但内容是制表符分隔。别信拓展名，看文件前几行就知道了。我的做法是先打印原始文本前几行，确认分隔符。

with open("results.txt") as f:
    for i in range(5):
        print(repr(f.readline()))

输出类似 'gene_id\tlog2FC\tp_value\n'，看到 \t 那就是 TSV。

4.3 必看字段：log2FC、p-value、padj

差异表达分析的结果里，一般至少有这几列：

• log2FC：log₂ 倍数变化。正数表示上调，负数表示下调。绝对值越大，表达差异越显著。
• p-value：原始 p 值。来自统计检验（如 Wald 检验、exact test）。
• padj：校正后的 p 值，通常是 BH（Benjamini-Hochberg）FDR 校正。用来控制假阳性率。
• base_mean：归一化后的平均表达量，用于过滤低表达基因。

做可视化时，我最常用的就是 log2FC 和 padj。举个例子，挑出显著差异基因：

# 筛选显著上调和下调的基因
significant = df[(df['padj'] < 0.05) & (abs(df['log2FC']) > 1)]
print(f"显著差异基因数量：{len(significant)}")

4.4 数据质量检查

数据读进来后，别急着画图。先跑几个检查：

1. 检查空值：df.isnull().sum()。如果 padj 或 log2FC 有 NaN，可能是低表达被过滤掉。需要决定是丢弃还是填充。
2. 检查列名：有些工具输出列名大小写不同，如 log2FoldChange 或 log2(Fold_change)。统一标准化列名。
3. 检查索引唯一性：基因 ID 应该是唯一的。重复索引会导致后续合并出错。
4. 检查数值范围：log2FC 通常集中在 -10 到 10，p-value 在 0 到 1 之间。如果出现 -Inf 或 Inf，要处理。

# 常用检查代码
print("缺失值统计：")
print(df.isnull().sum())
print("\n列名列表：")
print(df.columns.tolist())
print("\n是否有重复索引：", df.index.duplicated().any())

4.5 快速可视化概览

在正式绘制火山图或 MA 图之前，可以先做个分布检查，看看数据是否正常。

import matplotlib.pyplot as plt

# log2FC 分布直方图
df['log2FC'].hist(bins=50)
plt.title("Distribution of log2 Fold Change")
plt.xlabel("log2FC")
plt.ylabel("Frequency")
plt.show()

# p-value 的 -log10 转换
import numpy as np
df['-log10_pvalue'] = -np.log10(df['p_value'])
df['-log10_pvalue'].hist(bins=50)
plt.title("Distribution of -log10(p-value)")
plt.show()

这些图能帮你快速发现异常——比如 log2FC 全部集中在零附近，或者 p-value 分布严重偏态。

4.6 保存清理后的数据

检查完，建议把清洗后的数据存一份新文件。我通常这样处理：

# 丢弃含 padj 为空的基因（通常是被过滤掉的低表达基因）
df_clean = df.dropna(subset=['padj', 'log2FC'])

# 只保留分析需要的列
df_clean = df_clean[['gene_id', 'base_mean', 'log2FC', 'p_value', 'padj']]

# 保存为新的 CSV 文件
df_clean.to_csv("results_clean.csv", index=False)
print("已保存清洗后的数据")

我记得刚接触差异表达时，拿到一张几千行的表格就急着画火山图。结果发现 padj 全是 0.99，画出来一个上调的都没有……后来才意识到，有个参数忘了设置，导致统计检验压根没跑成功。

所以啊，读取数据和前期检查这一步，看似基础，实则决定了整张图是否可信。别嫌麻烦，多花五分钟检查，后面省一小时 debug。