数据下载与读取：从GEO/SRA到pandas/numpy
做Slide-seq分析的第一步，就是拿到原始数据。你可能觉得不就是下载嘛，点个链接就行。但真做起来，坑还是挺多的。我在带学生的时候，总有人卡在“为什么下载完打不开”这一步。
这一章我们就把整个过程捋清楚。从GEO/SRA怎么找数据，到用什么工具下载，最后怎么用pandas和numpy把数据读进内存。保证每一步都有代码，有经验，有坑。
先搞清楚数据在哪儿
Slide-seq的原始数据，大部分都存放在GEO（Gene Expression Omnibus）里。每次发文章，作者都会提交一个GSE号。比如经典的Slide-seqV2小鼠大脑数据，GSE号是GSE111111（我随便举的，别当真）。
你到了GEO页面，会看到几个关键区域：
- Supplementary files：这里一般是处理好的表达矩阵、barcodes、features之类，可以直接用pandas读。
- SRA run selector：这里指向SRA数据库，存放原始测序reads（fastq或sra格式）。如果你要做自己的pipeline，必须从这步开始。
我个人习惯是先看Supplementary files。如果人家已经给出了表达矩阵，我直接下载省时间。如果我要自己跑pipeline，或者文章没提供矩阵，才去SRA拉原始数据。
小提示：GEO页面底部通常有个“Series Matrix File(s)”链接，那是整理的样本信息表。我建议顺手下载下来，后面做metadata时很有用。
从SRA下载原始数据
SRA提供了两种原始数据格式：一是.sra（压缩二进制），二是.fastq（解压后的reads）。我推荐直接用prefetch和fastq-dump这套工具。
首先，你得安装SRA Toolkit。我用的是conda安装，稳得很：
conda install -c bioconda sra-tools然后找到一个SRA accession，比如“SRR11111111”。用prefetch下载到本地：
prefetch SRR11111111下载完是个.sra文件。我习惯把它转成fastq，用fastq-dump：
fastq-dump --split-files SRR11111111.sra--split-files是啥意思？如果你的是双端测序，不加这个参数会给你合并成一个文件。我刚开始做的时候不知道，跑完一看只有一个大文件，还以为数据丢了。后来查文档才明白。所以，嗯，这里要注意一下。
踩坑经验：我曾经用默认参数跑了100多G的Slide-seq数据，结果发现都是单端reads，实际上人家是双端的。浪费了好几个小时。后来我每次下载前，都会先去SRA页面看一眼文库类型，再决定加不加--split-files。
Slide-seq数据特有的格式问题
Slide-seq和10X Visium不一样。Visium的barcode是已知的，直接能匹配。但Slide-seq的barcode是随机的bead序列，你得自己去匹配。原始fastq里，reads的结构通常是：
- Read 1（25~30 bp）：包含bead barcode + UMI
- Read 2（50~100 bp）：基因序列
很多工具（比如STARsolo）能直接从fastq处理，不需要你自己拆。但我个人有个小习惯：先用seqtk或者python脚本快速检查一下reads的前几个碱基，确认barcode位置和长度。避免后面跑了半天发现参数设错了。
zcat SRR11111111.fastq.gz | head -8如果看到前25个碱基是ATCG随机组合，大概率就是barcode+UMI。
直接下载处理好的表达矩阵
如果你不想从fastq跑起，GEO上很多Slide-seq文章也直接提供了表达矩阵。常见的格式是：
- counts matrix：通常是.tsv或.csv，行是基因，列是bead barcode
- barcodes.tsv：每个bead的barcode列表
- features.tsv：基因名或ID
- bead metadata：例如每个bead的空间坐标x,y
我一般直接wget下来：
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE111nnn/suppl/GSE111111_RAW.tar
tar -xvf GSE111111_RAW.tar解压后就是几个文本文件。注意，有些文件是gzip压缩的，记得先解压：
gunzip *.gz用pandas和numpy加载数据
到这里才算进入正题。假设你拿到了表达矩阵（矩阵文件），以及对应的barcode和坐标文件。我们用pandas来读。
先读取barcode和基因列表：
import pandas as pd
import numpy as np
barcodes = pd.read_csv('barcodes.tsv', header=None, names=['barcode'])
features = pd.read_csv('features.tsv', header=None, names=['gene_id', 'gene_name'])
然后读表达矩阵。注意矩阵格式：有的是一行一个gene、一列一个barcode的稀疏矩阵（market matrix格式）。Slide-seq常用的是.mtx格式（Market Exchange Format）：
from scipy.io import mmread
matrix = mmread('matrix.mtx').tocsr()
如果你拿到的是稠密的.csv或.tsv，那更简单：
counts = pd.read_csv('counts.tsv', sep='\t', index_col=0)有些时候，矩阵里有很多列是空白的（没有bead附着）。我的做法是先过滤掉all-zero列：
counts = counts.loc[:, counts.sum(axis=0) > 0]完了加载坐标文件：
coords = pd.read_csv('bead_coords.csv')
coords.columns = ['barcode', 'x', 'y']
我个人习惯把barcode设成index，这样后面合并元数据更方便：
counts.columns = barcodes['barcode'].values
coords.set_index('barcode', inplace=True)
数据整合与快速检查
现在有了表达矩阵和坐标，我们可以做一个简单的Anndata对象（当然这一章我们先用numpy/pandas验证一下）。比如检查一下表达矩阵的维度，看看有多少beads和基因：
print(f"Matrix shape: {counts.shape}")
print(f"Number of beads: {counts.shape[1]}")
print(f"Number of genes: {counts.shape[0]}")
你想想看，如果矩阵的列数远小于barcode文件里的bead数量，那说明有很多beads没检测到表达。我遇到过这种情况，后来发现是因为原数据把低质量的beads提前过滤了。这不一定是错误，但需要你在后续分析中注意。
再检查一下坐标范围：
print(coords.describe())如果x,y都集中在0~10000的范围，那基本正常。如果出现负坐标或者NaN，可能数据格式有问题。我曾经下载过一次老数据，坐标列名写成了“x_position”而不是“x”，结果合并时对不上。后来我用列名模糊匹配才找到。所以，读取之前先 print(coords.columns) 是一个好习惯。
常见数据格式对比
为了让你快速区分不同的下载方式，我整理了一个表：
| 数据来源 | 格式 | 读取方式 | 优点 | 缺点 |
|---|---|---|---|---|
| SRA原始fastq | .fastq / .fastq.gz | 需pipeline预处理 | 完全控制预处理 | 慢，需较大计算资源 |
| GEO表达矩阵（tsv） | .tsv / .csv | pd.read_csv | 快速获得结果 | 可能被过滤过 |
| GEO表达矩阵（mtx） | .mtx + barcodes + features | scipy.io.mmread | 稀疏节省内存 | 需额外合并坐标 |
| GEO metadata（坐标csv） | .csv | pd.read_csv | 小巧易读 | 有时列名不规范 |
如果你只想做下游分析，建议直接用GEO提供的矩阵+坐标。如果你要自己跑定量，必须走SRA raw fastq路线。
本章知识：从原始数据到可分析矩阵
上面的图展示了两种常见路线。左侧是完整pipeline，适合需要定制预处理的高端用户；右侧是快速上手路线，我推荐初学者先走右侧。
还有一些实用技巧
最后补充几个我经常用的小技能：
- 内存不足怎么办？ Slide-seq的矩阵可能很大（几十万beads × 两万基因）。我一般用
pd.read_csv(..., chunksize=10000)分块读，或者转成scipy.sparse再操作。 - 坐标和barcode对应不上？ 我曾遇到过barcode文件里有很多“-1”后缀，但坐标文件里没有。用 .str.replace('-1', '') 清洗一下：
barcodes['barcode'] = barcodes['barcode'].str.replace('-1', '')- 文件编码问题：有些GEO数据用UTF-16编码，pandas读出来全是乱码。加个参数：
encoding='utf-16-le'就好了。
嗯，数据下载和读取这块看起来琐碎，但打好底子后面分析才能顺。你按照我上面列的那几个步骤走一遍，基本不会出大问题。