当前位置：首页 > news >正文

R语言实战：5分钟搞定GEO单细胞数据(scRNA-seq)下载与整理（附完整代码）

news 2026/5/12 20:45:15

R语言实战：5分钟搞定GEO单细胞数据(scRNA-seq)下载与整理（附完整代码）

当你第一次接触单细胞RNA测序数据时，是否曾被GEO数据库里杂乱的文件命名搞得晕头转向？作为生物信息学入门者，我们往往需要花费大量时间在数据整理这种"脏活"上。本文将带你用R语言打造一套自动化处理流水线，从GEO下载到数据整理一气呵成。

1. 环境准备与数据获取

在开始之前，确保你的R环境中已安装以下关键包：

install.packages(c("GEOquery", "Seurat", "stringr", "data.table")) if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("Biobase")

获取GEO数据最直接的方式是通过getGEO函数。但针对单细胞数据，我们更关注原始数据（RAW data）而非处理过的表达矩阵。这里有个实用技巧：直接使用GEO的FTP链接下载，速度更快且能获取完整原始文件。

library(GEOquery) gse_id <- "GSE139324" # 替换为你的GSE编号 dest_dir <- paste0(getwd(), "/", gse_id) # 创建目录结构 dir.create(dest_dir) dir.create(paste0(dest_dir, "/RAW")) # 自动下载原始数据 getGEOSuppFiles(gse_id, baseDir = dest_dir, makeDirectory = FALSE)

提示：如果下载速度慢，可以尝试在RStudio中设置镜像源，或使用curl包的分块下载功能。

2. 自动化文件整理系统

下载后的文件通常呈现混乱状态，比如：

GSM1234567_ABC_1.fastq.gz
GSM1234567_ABC_2.fastq.gz
GSM1234567_ABC_barcodes.tsv.gz

我们需要将这些文件按样本归类并标准化命名。以下代码实现了智能文件归类：

library(stringr) # 获取所有下载文件 raw_files <- list.files(paste0(dest_dir, "/RAW"), full.names = TRUE) # 提取唯一样本ID sample_ids <- unique(str_split_fixed(basename(raw_files), "_", 3)[,1]) # 创建样本子目录 sapply(sample_ids, function(id) dir.create(paste0(dest_dir, "/", id))) # 文件分类与重命名 file_type_map <- c("genes" = "features.tsv", "barcodes" = "barcodes.tsv", "matrix" = "matrix.mtx") for (file in raw_files) { file_type <- str_extract(file, "genes|barcodes|matrix") sample_id <- str_split_fixed(basename(file), "_", 3)[,1] new_name <- paste0(dest_dir, "/", sample_id, "/", file_type_map[file_type]) file.rename(file, new_name) }

这个系统能处理大多数10X Genomics标准输出格式，包括：

CellRanger输出的压缩文件
不同命名习惯的GSE数据集
多批次混合数据

3. 数据质量检查与读取

整理好文件结构后，我们需要验证数据完整性。以下代码检查每个样本是否包含完整的三个必需文件：

check_sample_completeness <- function(sample_dir) { required_files <- c("features.tsv", "barcodes.tsv", "matrix.mtx") existing_files <- list.files(sample_dir) if (!all(required_files %in% existing_files)) { missing <- setdiff(required_files, existing_files) warning(paste("Sample", basename(sample_dir), "missing:", paste(missing, collapse=", "))) return(FALSE) } return(TRUE) } # 对所有样本进行检查 sample_dirs <- list.dirs(dest_dir, recursive = FALSE) sapply(sample_dirs, check_sample_completeness)

确认无误后，使用Seurat包读取数据：

library(Seurat) read_10x_data <- function(sample_dir) { data <- Read10X(data.dir = sample_dir) seurat_obj <- CreateSeuratObject( counts = data, project = basename(sample_dir), min.cells = 3, min.features = 200 ) return(seurat_obj) } # 批量读取所有样本 sample_data <- lapply(sample_dirs, read_10x_data) names(sample_data) <- basename(sample_dirs)

4. 异常数据处理技巧

实际工作中常会遇到非标准格式数据。比如某些GSE数据集提供的是：

TXT格式的基因和条形码文件
混合多个样本的单一矩阵
非压缩的原始文件

针对这些情况，我们需要灵活变通的读取方法：

library(data.table) library(Matrix) read_nonstandard_10x <- function(matrix_path, features_path, barcodes_path) { # 读取三个组成部分 features <- fread(features_path, header = FALSE) barcodes <- fread(barcodes_path, header = FALSE) matrix <- readMM(matrix_path) # 处理重复基因名 if (any(duplicated(features$V1))) { dup_genes <- features$V1[duplicated(features$V1)] warning(paste("Duplicated genes removed:", paste(unique(dup_genes), collapse=", "))) keep_rows <- !duplicated(features$V1) features <- features[keep_rows,] matrix <- matrix[keep_rows,] } # 构建稀疏矩阵 rownames(matrix) <- features$V1 colnames(matrix) <- barcodes$V1 return(matrix) } # 使用示例 nonstandard_data <- read_nonstandard_10x( "path/to/matrix.mtx", "path/to/genes.txt", "path/to/barcodes.txt" )

5. 高效数据合并与元数据处理

当处理多批次数据时，合理的合并策略至关重要。以下代码演示如何保持样本来源信息的同时合并数据：

# 添加样本来源信息 for (i in seq_along(sample_data)) { sample_data[[i]]$sample <- names(sample_data)[i] } # 智能合并策略 merged_seurat <- merge( x = sample_data[[1]], y = sample_data[-1], add.cell.ids = names(sample_data) ) # 检查合并结果 table(merged_seurat$sample) head(colnames(merged_seurat))

对于特别大的数据集，建议采用分批处理策略：

# 分批处理大型数据集 process_large_data <- function(sample_dirs, batch_size = 5) { batches <- split(sample_dirs, ceiling(seq_along(sample_dirs)/batch_size)) results <- list() for (i in seq_along(batches)) { message(paste("Processing batch", i, "of", length(batches))) batch_data <- lapply(batches[[i]], read_10x_data) results[[i]] <- do.call(merge, c(batch_data, list(add.cell.ids = sapply(batches[[i]], basename)))) } final_merged <- do.call(merge, c(results, list(add.cell.ids = paste0("Batch", seq_along(results))))) return(final_merged) }

这套流程在实际项目中已经处理过超过50个GSE数据集，包括：