📚 课程目录第2章 L1000数据格式与读取方法

2 L1000 数据格式与读取方法归纳总结

来，咱们直接切入正题。L1000 项目的原始数据，说白了就是一堆基因表达值，但它不是普通的那种高通量芯片结果。它的核心是用了一组“条形码基因”(landmark genes) 来推算全基因组表达谱。

我记得刚接触这个数据时，第一反应就是：这文件后缀怎么这么奇怪？.gct、.grp、.nctx... 别急，我一个个给你拆开讲。

2.1 核心数据标准：GCT 与 GRP

L1000 数据的基石，是两个文件格式：GCT （Gene Cluster Text）和 GRP （Gene Set File）。你想想看，没有这两个文件，后面所有分析都是空中楼阁。

格式	全称	用途	我遇到过的大坑
GCT	Gene Cluster Text	存储基因表达矩阵（行是基因，列是样本/处理条件）	表头格式必须严格带 `#1.3`，少个换行符都读不了
GRP	Gene Set File	存储基因列表（一行一个基因符号）	有些文件末尾有空行，解析时会多出一个空字符串

2.2 GCT 文件结构

GCT 文件有固定的 5 行头，然后才是数据。给我记住这个规矩，不然解析会出错。

#1.3
数量 (行)  数量 (列)
基因 ID1  基因符号1  样本1  样本2  样本3  ...
基因 ID2  基因符号2  val    val    val    ...
...

举个例子你就明白了。我手头有个简单版：

#1.3
3   4
ID_A   GENE1  12.3  14.5  9.8  10.1
ID_B   GENE2  5.6   7.2   6.1  5.9
ID_C   GENE3  20.1  19.8  21.5  22.0

重点：第二行是矩阵大小，但 GCT 的 #1.3 版本中，这个行数是除表头外的数据行数。我记得有一次用 pandas 直接读，参数没设置好，把表头也当数据解析了，结果矩阵维度全乱。嗯，这里要注意。

2.3 GRP 文件：简单但容易翻车

GRP 文件就是纯文本，每行一个基因 ID。看起来简单，但有些文件末尾有个额外的换行或空行。

避坑提醒：我曾经用 readlines 读完 GRP 文件后直接遍历，结果最后一行是空字符串，导致下游匹配基因时总差一个。解决方案很简单：line.strip() 后过滤空字符串。

# 示例 GRP 内容
AFFX-TrpnX-3_at
AFFX-ThrX-5_at
AFFX-ThrX-3_at
AFFX-DapX-5_at

2.4 读取方法：Python 实操

我建议你直接用 cmapPy 这个库，它是 LINCS 项目官方的 Python 工具包，专门解析 GCT 和 GRP。省事，也规范。

我的习惯：如果只是快速瞄一眼数据，我会用 pandas 临时顶一下；但一旦做正式分析，必须上 cmapPy。

2.4.1 安装 cmapPy

pip install cmappy

2.4.2 读取 GCT 文件

from cmapPy.pandasGEXpress import parse

# 加载 GCT 文件
gctoo = parse("your_file.gct")

# 查看数据框
print(gctoo.data_df.head())

# 获取基因列表
genes = gctoo.data_df.index.tolist()
print(f"基因数量: {len(genes)}")

# 获取样本列表
samples = gctoo.data_df.columns.tolist()
print(f"样本数量: {len(samples)}")

2.4.3 读取 GRP 文件

from cmapPy.pandasGEXpress import parse_grp

grp_genes = parse_grp("your_gene_set.grp")
print(f"基因集合大小: {len(grp_genes)}")
print(grp_genes[:10])  # 只看前 10 个

2.4.4 手写解析（不依赖 cmapPy）

有些环境不让装 cmapPy？别慌，其实核心逻辑不复杂。但我要警告你：纯手写容易漏掉 GCT 表头格式的细节。

def read_gct_simple(filepath):
    with open(filepath, "r") as f:
        lines = f.readlines()
    
    # 第一行校验
    assert lines[0].strip() == "#1.3", "不是标准 GCT 文件头"
    
    # 第二行：行数和列数
    n_rows, n_cols = map(int, lines[1].strip().split())
    
    # 从第三行开始才是列名（样本名）
    headers = lines[2].strip().split()
    
    # 数据从第四行开始
    data_lines = lines[3:]
    
    # 我用 pandas 构建矩阵
    import pandas as pd
    data = []
    for line in data_lines:
        parts = line.strip().split()
        # 前两列是基因 ID 和基因符号，后面是表达值
        data.append(parts)
    
    df = pd.DataFrame(data, columns=headers)
    df.set_index(headers[0], inplace=True)
    return df

# 调用
gex_df = read_gct_simple("sample.gct")

核心提醒：手写版本我没做类型转换，读出来的表达值全是字符串。你记得加上 df.iloc[:, 2:] = df.iloc[:, 2:].astype(float)。

2.5 整体数据读取流程

L1000 的数据读取不是一步到位的。我个人习惯是先读 GRP 获取 landmark 基因列表，再读 GCT 过滤出这些基因的表达矩阵。下面这张图能帮你理清顺序：

说白了，读取 L1000 数据，重点就两块：格式认知和工具选择。GCT 的 5 行头结构、GRP 的纯列表结构，都是固定的。我个人强烈推荐用 cmapPy，它能自动处理版本差异和基因名匹配。

实用技巧：有些 GCT 文件是压缩包（.gctx），用 cmapPy 的 parse 函数一样能读，它会自动解压。省去你手动 gunzip 的麻烦。

嗯，关于 L1000 的数据格式和读取方法，我心里有数的几块就这些。实际项目里，格式合规性是第一个拦路虎，把这个啃下来，后面解析就不太会出幺蛾子了。