简介这份资源是面向生物信息学入门者与单细胞研究初学者的数据查找指南配套项目代码聚焦单细胞转录组学中公共数据获取这一关键环节帮助读者解决从海量数据库中精准定位可用数据集的问题。压缩包共3个文件以inscode项目配置、html页面与gitignore为主整体约9KB结构轻量便于快速浏览与二次整理。内容围绕GEO、Single Cell Portal、Human Cell Atlas等常用数据库展开涵盖布尔逻辑检索、数据分类筛选、数据格式识别、质量检查及下载许可与引用规则等要点并附有作业要求说明可帮助读者建立从查找、下载到初步应用的完整认知。目前已有131人学习适合希望掌握单细胞数据检索思路、提升数据处理技能的研究者参考使用。1. 单细胞转录组数据查找指南从项目代码到可复现的数据获取路径做过单细胞的人都有一个共识分析卡住的时候八成不是算法不会写而是数据没找对。你手里可能已经有一份项目代码里面写着read.csv(data/expr_matrix.csv)但那个data/目录是空的——原作者没传或者传了一个你根本下不动的压缩包。单细胞转录组数据查找指南要解决的核心问题就是这个给定一份项目代码怎么反推出它需要什么格式的数据去哪里找怎么验证找对了。适合两类人一是拿到别人项目代码想复现但缺数据的研究生和工程师二是自己做完分析想整理成可共享项目、但不知道数据该怎么组织的人。这篇不聊空泛的数据库列表只讲从代码出发、反向定位数据、下载后校验的完整链路。2. 从项目代码反推数据需求先读懂代码在找什么2.1 项目代码里数据引用的四种典型写法拿到一份单细胞项目代码第一件事不是跑而是搜。搜什么搜所有跟文件读取相关的函数调用。常见的有四类第一类是直接读矩阵文件比如read.csv()、read.table()、read.delim()后面跟的路径就是表达矩阵。第二类是读 RDS 或 h5ad 这类序列化对象比如readRDS(pbmc_processed.rds)或sc.read_h5ad(adata.h5ad)这种通常已经是做完质控和注释的完整对象。第三类是读 10x Genomics 的标准输出目录代码里会出现Read10X(data.dir filtered_feature_bc_matrix)这样的调用。第四类是从公共数据库的 API 拉取比如GEOquery::getGEO()或者scanpy.datasets.pbmc3k()。# 在 R 项目里快速定位所有数据读取调用 # 在项目根目录执行列出所有涉及文件读取的行 grep -rn read\.csv\|read\.table\|readRDS\|Read10X\|read_h5ad\|read_10x \ --include*.R --include*.Rmd --include*.py .这段命令的作用是把项目里所有可能涉及数据读取的代码行列出来-r表示递归搜索子目录-n显示行号--include限定文件类型。跑完之后你会得到一张清单每一行都指向一个数据依赖。参数方面如果你的项目同时有 R 和 Python 代码把--include加上*.py和*.ipynb如果项目很大可以加--exclude-dirrenv跳过依赖包目录。2.2 判断数据是原始测序数据还是处理后的对象搜到读取调用之后下一步是判断这个项目需要的是原始数据还是已经处理好的对象。判断依据看文件扩展名和变量名.fastq.gz是原始测序读长.bam是比对结果barcodes.tsv features.tsv matrix.mtx是 10x 的计数矩阵.rds和.h5ad通常是已经做完质控、降维、聚类的对象。这个判断直接决定你去哪里找数据。如果代码里是Read10X()你需要去 GEO 或 ENA 找 10x 原始输出如果是readRDS()你需要找作者上传的处理后对象或者自己从原始数据跑一遍完整流程。常见做法是先看项目有没有data/目录下的 README 或sessionInfo()输出里面往往会写清楚数据来源和版本。# 在 Python 项目里检查 AnnData 对象的来源线索 import scanpy as sc import h5py # 如果 h5ad 文件已存在先看它的 unstructured 元数据 with h5py.File(adata.h5ad, r) as f: # uns 里常存有数据来源、处理参数等信息 if uns in f: for key in f[uns].keys(): print(key, f[uns][key]) # obs 的列名能告诉你做了哪些注释 if obs in f: print(obs columns:, list(f[obs].keys()))这段代码不加载整个对象只用h5py读元数据速度快且省内存。uns是 AnnData 的 unstructured 字典作者如果规范会把数据来源、参考文献、处理参数放在这里。obs的列名能告诉你这份数据已经做了哪些注释比如有cell_type列说明已经注释过细胞类型有leiden列说明做过聚类。参数上如果文件很大h5py的懒加载特性让你不用把整个矩阵读进内存。2.3 用依赖清单锁定数据版本和参考基因组单细胞项目对版本极其敏感。同一份数据用 GRCh38 和 GRCh37 比对出来的矩阵不一样Seurat v4 和 v5 读同一个 h5 文件可能报不同的错。所以从代码反推数据需求时必须同时锁定三个版本参考基因组版本、细胞注释参考集版本、分析工具版本。# 从项目里提取版本信息 # R 项目看 renv.lock 或 sessionInfo cat renv.lock 2/dev/null | head -50 # Python 项目看 requirements.txt 或 environment.yml cat requirements.txt 2/dev/null || cat environment.yml 2/dev/null # 搜索代码里出现的基因组版本关键词 grep -rn GRCh38\|GRCh37\|mm10\|hg19\|hg38 --include*.R --include*.py .这三条命令分别对应 R 项目的依赖锁文件、Python 项目的依赖清单、以及代码里硬编码的基因组版本。renv.lock是 R 的依赖锁定文件里面记录了每个包的精确版本requirements.txt和environment.yml同理。最后一条 grep 是关键——很多项目不会在 README 里写基因组版本但代码里会硬编码比如reference GRCh38。把这些信息凑齐你才知道该去下载哪个版本的数据。3. 公共数据源的实际检索策略GEO、ArrayExpress 与单细胞专用库3.1 GEO 检索的五个有效过滤条件GEO 是最常用的公共表达数据仓库但直接搜关键词往往返回几百个结果。有效的做法是用过滤条件逐步收窄。我一般会同时用这五个条件研究类型选Expression profiling by high throughput sequencing物种选目标物种样本数设下限单细胞研究通常不少于 3 个样本时间范围限定近五年然后加single cell或scRNA作为关键词。# 用 NCBI E-utilities 命令行检索 GEO 单细胞数据集 # 搜索人类单细胞转录组、样本数大于 5、近三年的数据集 esearch -db gds -query \ single cell RNA seq[Title] AND Homo sapiens[Organism] \ AND 2022:2024[PDAT] AND gse[Entry Type] \ | efetch -format docsum \ | xtract -pattern DocumentSummary \ -element Accession,Title,n_samples \ | head -30esearch是 NCBI 的检索入口-db gds指定搜 GEO DataSets 库。查询语法里[Title]限定标题字段[Organism]限定物种[PDAT]限定发布日期[Entry Type]限定条目类型为 GSE数据集而非 GPL平台或 GSM样本。efetch -format docsum拉取摘要xtract按 XML 结构提取 Accession、Title 和样本数三列。参数上head -30只显示前 30 条实际使用时可以去掉或改成 results.txt存文件。3.2 单细胞专用库的适用场景对比GEO 虽然大但对单细胞数据来说有几个专用库更高效。下表是我实际用过的四个库的对比数据库数据规模主要优势适用场景GEO最大原始数据处理后数据都有找特定研究的数据ArrayExpress中等欧洲研究居多元数据规范补充 GEO 未收录的数据单细胞图谱类库大统一处理的标准化矩阵做跨样本比较细胞注释参考库中等带细胞类型标签做注释迁移选择逻辑是如果你要复现某篇论文优先去 GEO 或 ArrayExpress 找对应 Accession如果你要做大规模跨样本分析用单细胞图谱类库的标准化矩阵省去自己处理的麻烦如果你需要带标签的数据做监督学习去细胞注释参考库。常见做法是先在 GEO 搜找不到再去 ArrayExpress两者都没有再考虑专用库。3.3 下载 10x 数据的目录结构与校验方法从 GEO 下载 10x 数据时经常遇到的一个问题是下载下来一堆文件但不知道哪个是barcodes.tsv、哪个是matrix.mtx。10x 的标准输出是三个文件放在一个目录里目录名通常是filtered_feature_bc_matrix或raw_feature_bc_matrix。# 下载并整理 10x 数据目录 # 假设从 GEO 下载了一个 tar 包 wget -O GSEXXXXX_RAW.tar https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/GSEXXXXX_RAW.tar # 解压 mkdir -p data/GSEXXXXX tar -xvf GSEXXXXX_RAW.tar -C data/GSEXXXXX # 查看解压后的文件结构 ls -lh data/GSEXXXXX/ # 如果文件是分散的按样本整理成 10x 标准目录 for sample in GSM*; do mkdir -p data/GSEXXXXX/$sample mv $sample*_barcodes.tsv.gz data/GSEXXXXX/$sample/barcodes.tsv.gz 2/dev/null mv $sample*_features.tsv.gz data/GSEXXXXX/$sample/features.tsv.gz 2/dev/null mv $sample*_matrix.mtx.gz data/GSEXXXXX/$sample/matrix.mtx.gz 2/dev/null done这段脚本做了三件事下载 tar 包、解压到指定目录、把分散的文件按样本整理成 10x 标准结构。wget -O指定输出文件名tar -xvf解压并显示过程-C指定解压目标目录。后面的 for 循环遍历所有 GSM 开头的文件按样本建子目录把对应的 barcodes、features、matrix 文件移进去并重命名。参数上2/dev/null是忽略文件不存在的报错因为有些样本可能只有部分文件。校验方法是整理完之后用Read10X()试读一下能读出矩阵且维度合理基因数在 20000-40000 之间细胞数在几百到几万之间就说明结构对了。4. 数据下载后的格式转换与质量校验4.1 从 GEO 矩阵到 Seurat 对象的完整转换GEO 上很多单细胞数据是以表达矩阵的形式提供的通常是一个大的 txt 或 csv 文件行是基因、列是细胞。这种格式不能直接喂给 Seurat需要先转成稀疏矩阵再建对象。# 从 GEO 表达矩阵构建 Seurat 对象 library(Seurat) library(Matrix) # 读取矩阵第一列是基因名 expr - read.csv(GSEXXXXX_expression_matrix.csv, row.names 1, check.names FALSE) # 转成稀疏矩阵单细胞数据稀疏度很高用稀疏矩阵省内存 expr_sparse - as(as.matrix(expr), dgCMatrix) # 创建 Seurat 对象设最小细胞数和最小基因数过滤低质量 seu - CreateSeuratObject(counts expr_sparse, project GSEXXXXX, min.cells 3, min.features 200) # 查看对象基本信息 print(seu) # 计算线粒体基因比例用于后续质控 seu[[percent.mt]] - PercentageFeatureSet(seu, pattern ^MT-)read.csv的row.names 1表示第一列作为行名基因名check.names FALSE防止 R 把基因名里的特殊字符自动转换。as.matrix再转dgCMatrix是为了把密集矩阵转成稀疏矩阵单细胞数据的稀疏度通常在 90% 以上用稀疏矩阵能省大量内存。CreateSeuratObject的min.cells 3表示一个基因至少在 3 个细胞里表达才保留min.features 200表示一个细胞至少检测到 200 个基因才保留。PercentageFeatureSet用^MT-正则匹配线粒体基因算出来的比例存在percent.mt列里后续质控要用。4.2 用 DoubletFinder 和 Scrublet 做双细胞检测双细胞doublet是单细胞数据里最容易被忽略的坑。两个细胞被同一个液滴包裹测序结果看起来像一个细胞但表达谱是混合的。如果不检测聚类时会出现一些“四不像”的群。# 用 Scrublet 检测双细胞 import scrublet as scr import scanpy as sc # 读取数据 adata sc.read_10x_mtx(data/filtered_feature_bc_matrix/) # 初始化 Scrubletexpected_doublet_rate 根据细胞数估算 scrub scr.Scrublet(adata.X, expected_doublet_rate0.06) # 运行检测 doublet_scores, predicted_doublets scrub.scrub_doublets( min_counts2, min_cells3, min_gene_variability_pctl85, n_prin_comps30 ) # 把结果存回 AnnData adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets # 查看预测的双细胞比例 print(fPredicted doublet rate: {predicted_doublets.sum() / len(predicted_doublets):.2%})expected_doublet_rate是关键参数10x 官方给的经验值是每 1000 个细胞约 0.8%但实际项目里我一般设 0.06 到 0.08 之间因为很多数据集细胞数在 5000 到 10000 之间。min_counts2表示一个基因至少在 2 个计数里出现才参与模拟min_cells3同理。n_prin_comps30是 PCA 主成分数影响模拟的准确性。跑完之后看预测比例如果超过 15% 就要警惕可能是参数设得太激进或者数据本身质量有问题。4.3 校验数据完整性的三个检查点下载完数据、转完格式之后别急着往下跑分析。先做三个检查第一检查基因名和细胞名有没有重复或缺失第二检查矩阵的稀疏度和总计数是否合理第三检查线粒体基因比例分布是否正常。# 数据完整性检查 # 检查点 1基因名和细胞名 cat(基因数:, nrow(seu), 细胞数:, ncol(seu), \n) cat(重复基因名:, sum(duplicated(rownames(seu))), \n) cat(重复细胞名:, sum(duplicated(colnames(seu))), \n) # 检查点 2矩阵稀疏度和总计数 total_counts - Matrix::colSums(seuassays$RNAcounts) cat(总计数中位数:, median(total_counts), \n) cat(总计数范围:, range(total_counts), \n) # 检查点 3线粒体比例分布 mt_summary - summary(seu$percent.mt) cat(线粒体比例摘要:\n) print(mt_summary) # 如果中位数超过 20%说明数据质量可能有问题第一个检查点看基因数和细胞数是否在合理范围重复名会导致后续分析报错。第二个检查点看总计数单细胞数据每个细胞的总计数中位数通常在 2000 到 5000 之间太低说明测序深度不够太高可能有双细胞。第三个检查点看线粒体比例正常样本的中位数在 5% 到 15% 之间超过 20% 说明细胞状态不好或者数据有污染。这三个检查跑完你才能确认数据是“能用”的。5. 避坑与常见问题数据查找路上的血泪经验5.1 下载的矩阵和代码里的基因名对不上现象代码里用rownames(seu)取基因名去匹配参考集结果匹配率不到 50%。原因通常是基因名版本不一致——有些数据集用 Ensembl ID如ENSG00000141510有些用 Gene Symbol如TP53还有些用带版本号的 Ensembl ID如ENSG00000141510.17。解决方法是先检查基因名格式然后用biomaRt或本地注释包做 ID 转换。# 检查基因名格式并转换 head(rownames(seu)) # 如果全是 ENSG 开头用 biomaRt 转 Symbol library(biomaRt) mart - useMart(ensembl, dataset hsapiens_gene_ensembl) gene_map - getBM(attributes c(ensembl_gene_id, hgnc_symbol), filters ensembl_gene_id, values rownames(seu), mart mart) # 合并转换结果去掉没有 Symbol 的 seu - seu[rownames(seu) %in% gene_map$ensembl_gene_id, ] gene_map - gene_map[match(rownames(seu), gene_map$ensembl_gene_id), ] rownames(seu) - gene_map$hgnc_symbol5.2 内存不够导致读取大矩阵时崩溃现象read.csv()读一个 2GB 的矩阵文件R 直接报cannot allocate vector of size。原因是read.csv先把整个文件读进内存再解析单细胞矩阵动辄几万基因乘几万细胞密集存储会爆内存。解决方法是改用data.table::fread()或Matrix::readMM()前者读取速度快且内存效率高后者直接读 Matrix Market 格式的稀疏矩阵。# 用 data.table 分块读取大矩阵 library(data.table) # fread 比 read.csv 快且省内存 expr_dt - fread(large_matrix.csv, data.table FALSE) # 或者用 Matrix 直接读 mtx 格式 library(Matrix) expr_sparse - readMM(matrix.mtx)5.3 细胞注释参考集版本不匹配现象用参考集做注释迁移结果大部分细胞被标成Unknown。原因是参考集和查询数据的基因集重叠度太低或者参考集是用不同平台如 Smart-seq2 vs 10x生成的。解决方法是先算基因重叠率低于 60% 就换参考集如果平台不同考虑用Seurat的FindTransferAnchors做跨平台整合而不是直接迁移。5.4 下载的 h5ad 文件版本不兼容现象sc.read_h5ad()报错ValueError: Unknown format version。原因是 h5ad 文件是用旧版 AnnData 写的新版读不了。解决方法是先看文件版本然后用对应版本的 AnnData 读或者用h5py手动降级读取。# 检查 h5ad 文件版本 import h5py with h5py.File(old_data.h5ad, r) as f: print(Format version:, f.attrs.get(encoding-type, unknown)) print(AnnData version:, f.attrs.get(version, unknown))5.5 数据下载到一半断连导致文件损坏现象tar -xvf解压时报unexpected EOF。原因是下载不完整文件大小和服务器上的不一致。解决方法是下载后先校验文件大小再用md5sum对比服务器提供的校验值如果有的话。没有校验值就重新下载用wget -c支持断点续传。# 断点续传下载 wget -c -O GSEXXXXX_RAW.tar https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/GSEXXXXX_RAW.tar # 校验文件完整性 md5sum GSEXXXXX_RAW.tar # 对比服务器提供的 md5如果有6. 把数据查找流程固化成可复用的脚本6.1 用配置文件管理数据源和参数每次找数据都手动搜、手动下载、手动整理效率太低。我的做法是写一个配置文件把数据源、下载地址、整理规则、校验参数都写进去然后用一个脚本统一执行。# data_config.yaml datasets: - name: pbmc_10k source: geo accession: GSEXXXXX download_url: https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/ files: - GSEXXXXX_RAW.tar format: 10x_mtx genome: GRCh38 min_cells: 3 min_features: 200 max_mt_percent: 20 - name: brain_atlas source: custom download_url: https://example.org/brain_atlas.h5ad format: h5ad genome: GRCh38这个配置文件把每个数据集的来源、下载地址、格式、基因组版本、质控参数都定义清楚。format字段决定用哪个读取函数genome字段决定后续比对和注释用哪个参考集min_cells和max_mt_percent是质控阈值。有了这个文件换一个项目只需要改配置不用改代码。6.2 一个脚本完成下载、解压、校验、转换# data_pipeline.py import yaml import subprocess import os import scanpy as sc import pandas as pd def load_config(config_path): with open(config_path) as f: return yaml.safe_load(f) def download_dataset(dataset): 下载数据集支持断点续传 url dataset[download_url] for fname in dataset.get(files, []): full_url url fname # -c 断点续传-q 静默模式-P 指定下载目录 subprocess.run([wget, -c, -q, -P, downloads/, full_url], checkTrue) return True def validate_download(dataset): 校验下载文件大小是否合理 for fname in dataset.get(files, []): path os.path.join(downloads, fname) size_mb os.path.getsize(path) / (1024 * 1024) print(f{fname}: {size_mb:.1f} MB) # 如果文件小于 1MB大概率是下载失败 if size_mb 1: raise ValueError(f{fname} 文件过小可能下载失败) return True def convert_to_adata(dataset): 根据格式转换成 AnnData 对象 fmt dataset[format] if fmt 10x_mtx: adata sc.read_10x_mtx(downloads/filtered_feature_bc_matrix/) elif fmt h5ad: adata sc.read_h5ad(downloads/ dataset[files][0]) elif fmt csv: df pd.read_csv(downloads/ dataset[files][0], index_col0) adata sc.AnnData(df.T) else: raise ValueError(f不支持的格式: {fmt}) # 统一基因名格式 adata.var_names_make_unique() return adata def qc_filter(adata, dataset): 按配置参数做质控 sc.pp.filter_cells(adata, min_genesdataset[min_features]) sc.pp.filter_genes(adata, min_cellsdataset[min_cells]) adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], inplaceTrue) adata adata[adata.obs[pct_counts_mt] dataset[max_mt_percent]] return adata if __name__ __main__: config load_config(data_config.yaml) for ds in config[datasets]: print(f处理数据集: {ds[name]}) download_dataset(ds) validate_download(ds) adata convert_to_adata(ds) adata qc_filter(adata, ds) # 保存处理后的对象 adata.write_h5ad(fprocessed/{ds[name]}.h5ad) print(f{ds[name]}: {adata.n_obs} 细胞, {adata.n_vars} 基因)这个脚本把整个流程串起来load_config读 YAML 配置download_dataset用 wget 断点续传下载validate_download检查文件大小convert_to_adata按格式读取并统一基因名qc_filter按配置参数做质控。每个函数只做一件事方便单独调试。参数方面min_genes对应配置里的min_featuresmin_cells对应min_cellsmax_mt_percent对应max_mt_percent。跑完之后每个数据集存成一个 h5ad 文件文件名就是配置里的name。6.3 用校验和日志确保每次跑的结果可追溯数据查找最怕的是“上次跑通了这次跑不通但不知道哪里变了”。解决办法是每次运行都记录日志和校验和。import hashlib import logging from datetime import datetime def setup_logging(name): logging.basicConfig( filenameflogs/{name}_{datetime.now():%Y%m%d_%H%M%S}.log, levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s ) return logging.getLogger(name) def file_checksum(path): 计算文件的 SHA256 校验和 sha256 hashlib.sha256() with open(path, rb) as f: for chunk in iter(lambda: f.read(8192), b): sha256.update(chunk) return sha256.hexdigest() # 在下载和转换之后记录校验和 logger setup_logging(data_pipeline) logger.info(f下载文件校验和: {file_checksum(downloads/GSEXXXXX_RAW.tar)}) logger.info(f处理后对象: {adata.n_obs} 细胞, {adata.n_vars} 基因)setup_logging按时间戳建日志文件file_checksum分块计算 SHA256 避免大文件爆内存。每次跑完把下载文件的校验和、处理后对象的维度记进日志。下次再跑对比日志就能知道是数据变了还是代码变了。这个习惯我坚持了两年帮我省了无数次“重新下载一遍试试”的时间。6.4 一个容易忽略的技巧先下小样本验证流程最后说一个我踩过坑之后养成的习惯不管数据多大先用一个小样本把整个流程跑通。GEO 上很多数据集提供suppl目录下的部分文件或者你可以只下载一个样本的 10x 数据。用这个小样本跑完下载、解压、读取、质控、转换的全流程确认没有报错、参数合理、输出正常再去下完整数据。这样即使流程有问题你浪费的只是几分钟而不是几小时。数据查找这件事快就是慢慢就是快——先把链路走通再批量执行比一上来就下几百 GB 然后发现格式不对要重来划算得多。希望帮到你。本文还有配套的精品资源点击获取
简介这份资源是面向生物信息学入门者与单细胞研究初学者的数据查找指南配套项目代码聚焦单细胞转录组学中公共数据获取这一关键环节帮助读者解决从海量数据库中精准定位可用数据集的问题。压缩包共3个文件以inscode项目配置、html页面与gitignore为主整体约9KB结构轻量便于快速浏览与二次整理。内容围绕GEO、Single Cell Portal、Human Cell Atlas等常用数据库展开涵盖布尔逻辑检索、数据分类筛选、数据格式识别、质量检查及下载许可与引用规则等要点并附有作业要求说明可帮助读者建立从查找、下载到初步应用的完整认知。目前已有131人学习适合希望掌握单细胞数据检索思路、提升数据处理技能的研究者参考使用。1. 单细胞转录组数据查找指南从项目代码到可复现的数据获取路径做过单细胞的人都有一个共识分析卡住的时候八成不是算法不会写而是数据没找对。你手里可能已经有一份项目代码里面写着read.csv(data/expr_matrix.csv)但那个data/目录是空的——原作者没传或者传了一个你根本下不动的压缩包。单细胞转录组数据查找指南要解决的核心问题就是这个给定一份项目代码怎么反推出它需要什么格式的数据去哪里找怎么验证找对了。适合两类人一是拿到别人项目代码想复现但缺数据的研究生和工程师二是自己做完分析想整理成可共享项目、但不知道数据该怎么组织的人。这篇不聊空泛的数据库列表只讲从代码出发、反向定位数据、下载后校验的完整链路。2. 从项目代码反推数据需求先读懂代码在找什么2.1 项目代码里数据引用的四种典型写法拿到一份单细胞项目代码第一件事不是跑而是搜。搜什么搜所有跟文件读取相关的函数调用。常见的有四类第一类是直接读矩阵文件比如read.csv()、read.table()、read.delim()后面跟的路径就是表达矩阵。第二类是读 RDS 或 h5ad 这类序列化对象比如readRDS(pbmc_processed.rds)或sc.read_h5ad(adata.h5ad)这种通常已经是做完质控和注释的完整对象。第三类是读 10x Genomics 的标准输出目录代码里会出现Read10X(data.dir filtered_feature_bc_matrix)这样的调用。第四类是从公共数据库的 API 拉取比如GEOquery::getGEO()或者scanpy.datasets.pbmc3k()。# 在 R 项目里快速定位所有数据读取调用 # 在项目根目录执行列出所有涉及文件读取的行 grep -rn read\.csv\|read\.table\|readRDS\|Read10X\|read_h5ad\|read_10x \ --include*.R --include*.Rmd --include*.py .这段命令的作用是把项目里所有可能涉及数据读取的代码行列出来-r表示递归搜索子目录-n显示行号--include限定文件类型。跑完之后你会得到一张清单每一行都指向一个数据依赖。参数方面如果你的项目同时有 R 和 Python 代码把--include加上*.py和*.ipynb如果项目很大可以加--exclude-dirrenv跳过依赖包目录。2.2 判断数据是原始测序数据还是处理后的对象搜到读取调用之后下一步是判断这个项目需要的是原始数据还是已经处理好的对象。判断依据看文件扩展名和变量名.fastq.gz是原始测序读长.bam是比对结果barcodes.tsv features.tsv matrix.mtx是 10x 的计数矩阵.rds和.h5ad通常是已经做完质控、降维、聚类的对象。这个判断直接决定你去哪里找数据。如果代码里是Read10X()你需要去 GEO 或 ENA 找 10x 原始输出如果是readRDS()你需要找作者上传的处理后对象或者自己从原始数据跑一遍完整流程。常见做法是先看项目有没有data/目录下的 README 或sessionInfo()输出里面往往会写清楚数据来源和版本。# 在 Python 项目里检查 AnnData 对象的来源线索 import scanpy as sc import h5py # 如果 h5ad 文件已存在先看它的 unstructured 元数据 with h5py.File(adata.h5ad, r) as f: # uns 里常存有数据来源、处理参数等信息 if uns in f: for key in f[uns].keys(): print(key, f[uns][key]) # obs 的列名能告诉你做了哪些注释 if obs in f: print(obs columns:, list(f[obs].keys()))这段代码不加载整个对象只用h5py读元数据速度快且省内存。uns是 AnnData 的 unstructured 字典作者如果规范会把数据来源、参考文献、处理参数放在这里。obs的列名能告诉你这份数据已经做了哪些注释比如有cell_type列说明已经注释过细胞类型有leiden列说明做过聚类。参数上如果文件很大h5py的懒加载特性让你不用把整个矩阵读进内存。2.3 用依赖清单锁定数据版本和参考基因组单细胞项目对版本极其敏感。同一份数据用 GRCh38 和 GRCh37 比对出来的矩阵不一样Seurat v4 和 v5 读同一个 h5 文件可能报不同的错。所以从代码反推数据需求时必须同时锁定三个版本参考基因组版本、细胞注释参考集版本、分析工具版本。# 从项目里提取版本信息 # R 项目看 renv.lock 或 sessionInfo cat renv.lock 2/dev/null | head -50 # Python 项目看 requirements.txt 或 environment.yml cat requirements.txt 2/dev/null || cat environment.yml 2/dev/null # 搜索代码里出现的基因组版本关键词 grep -rn GRCh38\|GRCh37\|mm10\|hg19\|hg38 --include*.R --include*.py .这三条命令分别对应 R 项目的依赖锁文件、Python 项目的依赖清单、以及代码里硬编码的基因组版本。renv.lock是 R 的依赖锁定文件里面记录了每个包的精确版本requirements.txt和environment.yml同理。最后一条 grep 是关键——很多项目不会在 README 里写基因组版本但代码里会硬编码比如reference GRCh38。把这些信息凑齐你才知道该去下载哪个版本的数据。3. 公共数据源的实际检索策略GEO、ArrayExpress 与单细胞专用库3.1 GEO 检索的五个有效过滤条件GEO 是最常用的公共表达数据仓库但直接搜关键词往往返回几百个结果。有效的做法是用过滤条件逐步收窄。我一般会同时用这五个条件研究类型选Expression profiling by high throughput sequencing物种选目标物种样本数设下限单细胞研究通常不少于 3 个样本时间范围限定近五年然后加single cell或scRNA作为关键词。# 用 NCBI E-utilities 命令行检索 GEO 单细胞数据集 # 搜索人类单细胞转录组、样本数大于 5、近三年的数据集 esearch -db gds -query \ single cell RNA seq[Title] AND Homo sapiens[Organism] \ AND 2022:2024[PDAT] AND gse[Entry Type] \ | efetch -format docsum \ | xtract -pattern DocumentSummary \ -element Accession,Title,n_samples \ | head -30esearch是 NCBI 的检索入口-db gds指定搜 GEO DataSets 库。查询语法里[Title]限定标题字段[Organism]限定物种[PDAT]限定发布日期[Entry Type]限定条目类型为 GSE数据集而非 GPL平台或 GSM样本。efetch -format docsum拉取摘要xtract按 XML 结构提取 Accession、Title 和样本数三列。参数上head -30只显示前 30 条实际使用时可以去掉或改成 results.txt存文件。3.2 单细胞专用库的适用场景对比GEO 虽然大但对单细胞数据来说有几个专用库更高效。下表是我实际用过的四个库的对比数据库数据规模主要优势适用场景GEO最大原始数据处理后数据都有找特定研究的数据ArrayExpress中等欧洲研究居多元数据规范补充 GEO 未收录的数据单细胞图谱类库大统一处理的标准化矩阵做跨样本比较细胞注释参考库中等带细胞类型标签做注释迁移选择逻辑是如果你要复现某篇论文优先去 GEO 或 ArrayExpress 找对应 Accession如果你要做大规模跨样本分析用单细胞图谱类库的标准化矩阵省去自己处理的麻烦如果你需要带标签的数据做监督学习去细胞注释参考库。常见做法是先在 GEO 搜找不到再去 ArrayExpress两者都没有再考虑专用库。3.3 下载 10x 数据的目录结构与校验方法从 GEO 下载 10x 数据时经常遇到的一个问题是下载下来一堆文件但不知道哪个是barcodes.tsv、哪个是matrix.mtx。10x 的标准输出是三个文件放在一个目录里目录名通常是filtered_feature_bc_matrix或raw_feature_bc_matrix。# 下载并整理 10x 数据目录 # 假设从 GEO 下载了一个 tar 包 wget -O GSEXXXXX_RAW.tar https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/GSEXXXXX_RAW.tar # 解压 mkdir -p data/GSEXXXXX tar -xvf GSEXXXXX_RAW.tar -C data/GSEXXXXX # 查看解压后的文件结构 ls -lh data/GSEXXXXX/ # 如果文件是分散的按样本整理成 10x 标准目录 for sample in GSM*; do mkdir -p data/GSEXXXXX/$sample mv $sample*_barcodes.tsv.gz data/GSEXXXXX/$sample/barcodes.tsv.gz 2/dev/null mv $sample*_features.tsv.gz data/GSEXXXXX/$sample/features.tsv.gz 2/dev/null mv $sample*_matrix.mtx.gz data/GSEXXXXX/$sample/matrix.mtx.gz 2/dev/null done这段脚本做了三件事下载 tar 包、解压到指定目录、把分散的文件按样本整理成 10x 标准结构。wget -O指定输出文件名tar -xvf解压并显示过程-C指定解压目标目录。后面的 for 循环遍历所有 GSM 开头的文件按样本建子目录把对应的 barcodes、features、matrix 文件移进去并重命名。参数上2/dev/null是忽略文件不存在的报错因为有些样本可能只有部分文件。校验方法是整理完之后用Read10X()试读一下能读出矩阵且维度合理基因数在 20000-40000 之间细胞数在几百到几万之间就说明结构对了。4. 数据下载后的格式转换与质量校验4.1 从 GEO 矩阵到 Seurat 对象的完整转换GEO 上很多单细胞数据是以表达矩阵的形式提供的通常是一个大的 txt 或 csv 文件行是基因、列是细胞。这种格式不能直接喂给 Seurat需要先转成稀疏矩阵再建对象。# 从 GEO 表达矩阵构建 Seurat 对象 library(Seurat) library(Matrix) # 读取矩阵第一列是基因名 expr - read.csv(GSEXXXXX_expression_matrix.csv, row.names 1, check.names FALSE) # 转成稀疏矩阵单细胞数据稀疏度很高用稀疏矩阵省内存 expr_sparse - as(as.matrix(expr), dgCMatrix) # 创建 Seurat 对象设最小细胞数和最小基因数过滤低质量 seu - CreateSeuratObject(counts expr_sparse, project GSEXXXXX, min.cells 3, min.features 200) # 查看对象基本信息 print(seu) # 计算线粒体基因比例用于后续质控 seu[[percent.mt]] - PercentageFeatureSet(seu, pattern ^MT-)read.csv的row.names 1表示第一列作为行名基因名check.names FALSE防止 R 把基因名里的特殊字符自动转换。as.matrix再转dgCMatrix是为了把密集矩阵转成稀疏矩阵单细胞数据的稀疏度通常在 90% 以上用稀疏矩阵能省大量内存。CreateSeuratObject的min.cells 3表示一个基因至少在 3 个细胞里表达才保留min.features 200表示一个细胞至少检测到 200 个基因才保留。PercentageFeatureSet用^MT-正则匹配线粒体基因算出来的比例存在percent.mt列里后续质控要用。4.2 用 DoubletFinder 和 Scrublet 做双细胞检测双细胞doublet是单细胞数据里最容易被忽略的坑。两个细胞被同一个液滴包裹测序结果看起来像一个细胞但表达谱是混合的。如果不检测聚类时会出现一些“四不像”的群。# 用 Scrublet 检测双细胞 import scrublet as scr import scanpy as sc # 读取数据 adata sc.read_10x_mtx(data/filtered_feature_bc_matrix/) # 初始化 Scrubletexpected_doublet_rate 根据细胞数估算 scrub scr.Scrublet(adata.X, expected_doublet_rate0.06) # 运行检测 doublet_scores, predicted_doublets scrub.scrub_doublets( min_counts2, min_cells3, min_gene_variability_pctl85, n_prin_comps30 ) # 把结果存回 AnnData adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets # 查看预测的双细胞比例 print(fPredicted doublet rate: {predicted_doublets.sum() / len(predicted_doublets):.2%})expected_doublet_rate是关键参数10x 官方给的经验值是每 1000 个细胞约 0.8%但实际项目里我一般设 0.06 到 0.08 之间因为很多数据集细胞数在 5000 到 10000 之间。min_counts2表示一个基因至少在 2 个计数里出现才参与模拟min_cells3同理。n_prin_comps30是 PCA 主成分数影响模拟的准确性。跑完之后看预测比例如果超过 15% 就要警惕可能是参数设得太激进或者数据本身质量有问题。4.3 校验数据完整性的三个检查点下载完数据、转完格式之后别急着往下跑分析。先做三个检查第一检查基因名和细胞名有没有重复或缺失第二检查矩阵的稀疏度和总计数是否合理第三检查线粒体基因比例分布是否正常。# 数据完整性检查 # 检查点 1基因名和细胞名 cat(基因数:, nrow(seu), 细胞数:, ncol(seu), \n) cat(重复基因名:, sum(duplicated(rownames(seu))), \n) cat(重复细胞名:, sum(duplicated(colnames(seu))), \n) # 检查点 2矩阵稀疏度和总计数 total_counts - Matrix::colSums(seuassays$RNAcounts) cat(总计数中位数:, median(total_counts), \n) cat(总计数范围:, range(total_counts), \n) # 检查点 3线粒体比例分布 mt_summary - summary(seu$percent.mt) cat(线粒体比例摘要:\n) print(mt_summary) # 如果中位数超过 20%说明数据质量可能有问题第一个检查点看基因数和细胞数是否在合理范围重复名会导致后续分析报错。第二个检查点看总计数单细胞数据每个细胞的总计数中位数通常在 2000 到 5000 之间太低说明测序深度不够太高可能有双细胞。第三个检查点看线粒体比例正常样本的中位数在 5% 到 15% 之间超过 20% 说明细胞状态不好或者数据有污染。这三个检查跑完你才能确认数据是“能用”的。5. 避坑与常见问题数据查找路上的血泪经验5.1 下载的矩阵和代码里的基因名对不上现象代码里用rownames(seu)取基因名去匹配参考集结果匹配率不到 50%。原因通常是基因名版本不一致——有些数据集用 Ensembl ID如ENSG00000141510有些用 Gene Symbol如TP53还有些用带版本号的 Ensembl ID如ENSG00000141510.17。解决方法是先检查基因名格式然后用biomaRt或本地注释包做 ID 转换。# 检查基因名格式并转换 head(rownames(seu)) # 如果全是 ENSG 开头用 biomaRt 转 Symbol library(biomaRt) mart - useMart(ensembl, dataset hsapiens_gene_ensembl) gene_map - getBM(attributes c(ensembl_gene_id, hgnc_symbol), filters ensembl_gene_id, values rownames(seu), mart mart) # 合并转换结果去掉没有 Symbol 的 seu - seu[rownames(seu) %in% gene_map$ensembl_gene_id, ] gene_map - gene_map[match(rownames(seu), gene_map$ensembl_gene_id), ] rownames(seu) - gene_map$hgnc_symbol5.2 内存不够导致读取大矩阵时崩溃现象read.csv()读一个 2GB 的矩阵文件R 直接报cannot allocate vector of size。原因是read.csv先把整个文件读进内存再解析单细胞矩阵动辄几万基因乘几万细胞密集存储会爆内存。解决方法是改用data.table::fread()或Matrix::readMM()前者读取速度快且内存效率高后者直接读 Matrix Market 格式的稀疏矩阵。# 用 data.table 分块读取大矩阵 library(data.table) # fread 比 read.csv 快且省内存 expr_dt - fread(large_matrix.csv, data.table FALSE) # 或者用 Matrix 直接读 mtx 格式 library(Matrix) expr_sparse - readMM(matrix.mtx)5.3 细胞注释参考集版本不匹配现象用参考集做注释迁移结果大部分细胞被标成Unknown。原因是参考集和查询数据的基因集重叠度太低或者参考集是用不同平台如 Smart-seq2 vs 10x生成的。解决方法是先算基因重叠率低于 60% 就换参考集如果平台不同考虑用Seurat的FindTransferAnchors做跨平台整合而不是直接迁移。5.4 下载的 h5ad 文件版本不兼容现象sc.read_h5ad()报错ValueError: Unknown format version。原因是 h5ad 文件是用旧版 AnnData 写的新版读不了。解决方法是先看文件版本然后用对应版本的 AnnData 读或者用h5py手动降级读取。# 检查 h5ad 文件版本 import h5py with h5py.File(old_data.h5ad, r) as f: print(Format version:, f.attrs.get(encoding-type, unknown)) print(AnnData version:, f.attrs.get(version, unknown))5.5 数据下载到一半断连导致文件损坏现象tar -xvf解压时报unexpected EOF。原因是下载不完整文件大小和服务器上的不一致。解决方法是下载后先校验文件大小再用md5sum对比服务器提供的校验值如果有的话。没有校验值就重新下载用wget -c支持断点续传。# 断点续传下载 wget -c -O GSEXXXXX_RAW.tar https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/GSEXXXXX_RAW.tar # 校验文件完整性 md5sum GSEXXXXX_RAW.tar # 对比服务器提供的 md5如果有6. 把数据查找流程固化成可复用的脚本6.1 用配置文件管理数据源和参数每次找数据都手动搜、手动下载、手动整理效率太低。我的做法是写一个配置文件把数据源、下载地址、整理规则、校验参数都写进去然后用一个脚本统一执行。# data_config.yaml datasets: - name: pbmc_10k source: geo accession: GSEXXXXX download_url: https://ftp.ncbi.nlm.nih.gov/geo/series/GSEXXnnn/GSEXXXXX/suppl/ files: - GSEXXXXX_RAW.tar format: 10x_mtx genome: GRCh38 min_cells: 3 min_features: 200 max_mt_percent: 20 - name: brain_atlas source: custom download_url: https://example.org/brain_atlas.h5ad format: h5ad genome: GRCh38这个配置文件把每个数据集的来源、下载地址、格式、基因组版本、质控参数都定义清楚。format字段决定用哪个读取函数genome字段决定后续比对和注释用哪个参考集min_cells和max_mt_percent是质控阈值。有了这个文件换一个项目只需要改配置不用改代码。6.2 一个脚本完成下载、解压、校验、转换# data_pipeline.py import yaml import subprocess import os import scanpy as sc import pandas as pd def load_config(config_path): with open(config_path) as f: return yaml.safe_load(f) def download_dataset(dataset): 下载数据集支持断点续传 url dataset[download_url] for fname in dataset.get(files, []): full_url url fname # -c 断点续传-q 静默模式-P 指定下载目录 subprocess.run([wget, -c, -q, -P, downloads/, full_url], checkTrue) return True def validate_download(dataset): 校验下载文件大小是否合理 for fname in dataset.get(files, []): path os.path.join(downloads, fname) size_mb os.path.getsize(path) / (1024 * 1024) print(f{fname}: {size_mb:.1f} MB) # 如果文件小于 1MB大概率是下载失败 if size_mb 1: raise ValueError(f{fname} 文件过小可能下载失败) return True def convert_to_adata(dataset): 根据格式转换成 AnnData 对象 fmt dataset[format] if fmt 10x_mtx: adata sc.read_10x_mtx(downloads/filtered_feature_bc_matrix/) elif fmt h5ad: adata sc.read_h5ad(downloads/ dataset[files][0]) elif fmt csv: df pd.read_csv(downloads/ dataset[files][0], index_col0) adata sc.AnnData(df.T) else: raise ValueError(f不支持的格式: {fmt}) # 统一基因名格式 adata.var_names_make_unique() return adata def qc_filter(adata, dataset): 按配置参数做质控 sc.pp.filter_cells(adata, min_genesdataset[min_features]) sc.pp.filter_genes(adata, min_cellsdataset[min_cells]) adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], inplaceTrue) adata adata[adata.obs[pct_counts_mt] dataset[max_mt_percent]] return adata if __name__ __main__: config load_config(data_config.yaml) for ds in config[datasets]: print(f处理数据集: {ds[name]}) download_dataset(ds) validate_download(ds) adata convert_to_adata(ds) adata qc_filter(adata, ds) # 保存处理后的对象 adata.write_h5ad(fprocessed/{ds[name]}.h5ad) print(f{ds[name]}: {adata.n_obs} 细胞, {adata.n_vars} 基因)这个脚本把整个流程串起来load_config读 YAML 配置download_dataset用 wget 断点续传下载validate_download检查文件大小convert_to_adata按格式读取并统一基因名qc_filter按配置参数做质控。每个函数只做一件事方便单独调试。参数方面min_genes对应配置里的min_featuresmin_cells对应min_cellsmax_mt_percent对应max_mt_percent。跑完之后每个数据集存成一个 h5ad 文件文件名就是配置里的name。6.3 用校验和日志确保每次跑的结果可追溯数据查找最怕的是“上次跑通了这次跑不通但不知道哪里变了”。解决办法是每次运行都记录日志和校验和。import hashlib import logging from datetime import datetime def setup_logging(name): logging.basicConfig( filenameflogs/{name}_{datetime.now():%Y%m%d_%H%M%S}.log, levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s ) return logging.getLogger(name) def file_checksum(path): 计算文件的 SHA256 校验和 sha256 hashlib.sha256() with open(path, rb) as f: for chunk in iter(lambda: f.read(8192), b): sha256.update(chunk) return sha256.hexdigest() # 在下载和转换之后记录校验和 logger setup_logging(data_pipeline) logger.info(f下载文件校验和: {file_checksum(downloads/GSEXXXXX_RAW.tar)}) logger.info(f处理后对象: {adata.n_obs} 细胞, {adata.n_vars} 基因)setup_logging按时间戳建日志文件file_checksum分块计算 SHA256 避免大文件爆内存。每次跑完把下载文件的校验和、处理后对象的维度记进日志。下次再跑对比日志就能知道是数据变了还是代码变了。这个习惯我坚持了两年帮我省了无数次“重新下载一遍试试”的时间。6.4 一个容易忽略的技巧先下小样本验证流程最后说一个我踩过坑之后养成的习惯不管数据多大先用一个小样本把整个流程跑通。GEO 上很多数据集提供suppl目录下的部分文件或者你可以只下载一个样本的 10x 数据。用这个小样本跑完下载、解压、读取、质控、转换的全流程确认没有报错、参数合理、输出正常再去下完整数据。这样即使流程有问题你浪费的只是几分钟而不是几小时。数据查找这件事快就是慢慢就是快——先把链路走通再批量执行比一上来就下几百 GB 然后发现格式不对要重来划算得多。希望帮到你。本文还有配套的精品资源点击获取
