# 单细胞转录组数据分析入门 - 2 - 各种格式的单细胞数据读取
单细胞转录组测序的格式可能有多种,包括单样本,多样本,10X 的标准输出文件,h5,h5ad,txt/csv/tsv
网上的教程鱼龙混杂,这使我在初次接触时晕头转向
# 10X 标准输出
这是目前比较主流的一种数据格式,对于一个样本,其包含:
barcodes.tsv.gz
:记录 cell idfeatures.tsv.gz
:记录 gene idmatrix.mtx.gz
:计数 counts 矩阵
三个文件,且命名必须与上面完全相同
直接使用 Seurat 的 Read10X()
函数读取即可
# 如何读取多个样本
通常,测序不会只测一个样本,那么如何读取多个样本呢?
以 GSE151177 的读取为例:
下载的数据为压缩包 GSE151177_RAW.tar
,解压到 GSE151177_RAW
文件夹中,此时所有样本的三个文件都在这个文件夹下
ls | |
目录: D:\Study\Project\GDG\scPsoriasis\data\GSE151177\GSE151177_RAW | |
Mode LastWriteTime Length Name | |
---- ------------- ------ ---- | |
-a---- 5/24/2020 10:59 AM 13776 GSM4567877_Control01_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567877_Control01_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 16419935 GSM4567877_Control01_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 5829 GSM4567878_Control02_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567878_Control02_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 3573641 GSM4567878_Control02_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 2155 GSM4567879_Control03_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567879_Control03_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 1441746 GSM4567879_Control03_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 8430 GSM4567880_Control04_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567880_Control04_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 6875043 GSM4567880_Control04_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 1669 GSM4567881_Control05_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567881_Control05_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 489487 GSM4567881_Control05_matrix.mtx.gz | |
-a---- 6/14/2020 1:30 AM 1494 GSM4567882_Control05F_barcodes.tsv.gz | |
-a---- 6/14/2020 1:30 AM 304728 GSM4567882_Control05F_features.tsv.gz | |
-a---- 6/14/2020 1:30 AM 147394 GSM4567882_Control05F_matrix.mtx.gz | |
-a---- 6/14/2020 1:53 AM 3889 GSM4567883_Psoriasis01_barcodes.tsv.gz | |
-a---- 6/14/2020 1:53 AM 304728 GSM4567883_Psoriasis01_features.tsv.gz | |
-a---- 6/14/2020 1:53 AM 3289901 GSM4567883_Psoriasis01_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 3970 GSM4567884_Psoriasis02_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567884_Psoriasis02_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 3318396 GSM4567884_Psoriasis02_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 3599 GSM4567885_Psoriasis02F_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567885_Psoriasis02F_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 1848116 GSM4567885_Psoriasis02F_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 6924 GSM4567886_Psoriasis03_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567886_Psoriasis03_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 7062387 GSM4567886_Psoriasis03_matrix.mtx.gz | |
-a---- 5/26/2020 11:44 AM 3037 GSM4567887_Psoriasis03F_barcodes.tsv.gz | |
-a---- 5/26/2020 11:44 AM 304728 GSM4567887_Psoriasis03F_features.tsv.gz | |
-a---- 5/26/2020 11:44 AM 975940 GSM4567887_Psoriasis03F_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 5366 GSM4567888_Psoriasis04_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567888_Psoriasis04_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 5613417 GSM4567888_Psoriasis04_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 7450 GSM4567889_Psoriasis04F_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567889_Psoriasis04F_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 3059527 GSM4567889_Psoriasis04F_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 5525 GSM4567890_Psoriasis05_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567890_Psoriasis05_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 5228209 GSM4567890_Psoriasis05_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 7983 GSM4567891_Psoriasis06_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567891_Psoriasis06_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 7441699 GSM4567891_Psoriasis06_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 4842 GSM4567892_Psoriasis06F_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567892_Psoriasis06F_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 5298708 GSM4567892_Psoriasis06F_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 4452 GSM4567893_Psoriasis07_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567893_Psoriasis07_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 4163977 GSM4567893_Psoriasis07_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 11775 GSM4567894_Psoriasis08_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567894_Psoriasis08_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 11177661 GSM4567894_Psoriasis08_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 15318 GSM4567895_Psoriasis09_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567895_Psoriasis09_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 15250717 GSM4567895_Psoriasis09_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 10836 GSM4567896_Psoriasis10_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567896_Psoriasis10_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 11705956 GSM4567896_Psoriasis10_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 9814 GSM4567897_Psoriasis11_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567897_Psoriasis11_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 10868239 GSM4567897_Psoriasis11_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 23278 GSM4567898_Psoriasis12_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567898_Psoriasis12_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 24708520 GSM4567898_Psoriasis12_matrix.mtx.gz | |
-a---- 5/24/2020 10:59 AM 12039 GSM4567899_Psoriasis13_barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 GSM4567899_Psoriasis13_features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 9919157 GSM4567899_Psoriasis13_matrix.mtx.gz |
可以看到,前缀相同的是同一个样本
根据 Seurat 包的 Read10X()
函数要求,我们需要整理一下文件,为每个样本单独建立一个文件夹,然后将样本对应的 3 个文件重命名一下
在 R 语言中,可以写一个函数来完成这些操作:
# 将文件按照前缀移动到相应的子文件夹,并去除前缀的函数 | |
organize_files <- function(datadir) { | |
# 列出文件夹中的所有文件 | |
files <- list.files(datadir, full.names = TRUE) | |
# 获取每个文件的 GSM 编号 | |
# sample_ids <- unique(sub("_.*", "", basename(files))) | |
# 如果你想要完整的前缀,而不只是 GSM 编号,适用于前缀比较短且规则的情况 | |
sample_ids <- unique(sub("_(barcodes|features|matrix).*", "", basename(files))) | |
# 遍历每个样本,创建对应文件夹并移动文件 | |
for (sample_id in sample_ids) { | |
# 创建样本文件夹 | |
sample_dir <- file.path(datadir, sample_id) | |
if (!dir.exists(sample_dir)) dir.create(sample_dir) | |
# 筛选属于当前样本的文件 | |
sample_files <- files[grepl(sample_id, files)] | |
# 根据文件类型移动到样本文件夹并重命名 | |
for (file in sample_files) { | |
new_name <- if (grepl("barcodes", file)) { | |
"barcodes.tsv.gz" | |
} else if (grepl("features", file)) { | |
"features.tsv.gz" | |
} else if (grepl("matrix", file)) { | |
"matrix.mtx.gz" | |
} else { | |
next # 跳过不匹配的文件 | |
} | |
# 移动并重命名文件 | |
file.rename(file, file.path(sample_dir, new_name)) | |
} | |
} | |
cat("文件已成功整理到各自样本文件夹中。\n") | |
} | |
# 使用示例 | |
datadir <- "D:/Study/Project/GDG/scPsoriasis/data/GSE151177/GSE151177_RAW" | |
sample_dirs <- organize_files(file.path(datadir, "GSE151177_RAW")) |
于是我们就获得了以下文件夹:
ls | |
目录: D:\Study\Project\GDG\scPsoriasis\data\GSE151177\GSE151177_RAW | |
Mode LastWriteTime Length Name | |
---- ------------- ------ ---- | |
d----- 9/14/2024 7:01 PM GSM4567877_Control01 | |
d----- 9/14/2024 7:01 PM GSM4567878_Control02 | |
d----- 9/14/2024 7:01 PM GSM4567879_Control03 | |
d----- 9/14/2024 7:01 PM GSM4567880_Control04 | |
d----- 9/14/2024 7:01 PM GSM4567881_Control05 | |
d----- 9/14/2024 7:01 PM GSM4567882_Control05F | |
d----- 9/14/2024 7:01 PM GSM4567883_Psoriasis01 | |
d----- 9/14/2024 7:01 PM GSM4567884_Psoriasis02 | |
d----- 9/14/2024 7:01 PM GSM4567885_Psoriasis02F | |
d----- 9/14/2024 7:01 PM GSM4567886_Psoriasis03 | |
d----- 9/14/2024 7:01 PM GSM4567887_Psoriasis03F | |
d----- 9/14/2024 7:01 PM GSM4567888_Psoriasis04 | |
d----- 9/14/2024 7:01 PM GSM4567889_Psoriasis04F | |
d----- 9/14/2024 7:01 PM GSM4567890_Psoriasis05 | |
d----- 9/14/2024 7:01 PM GSM4567891_Psoriasis06 | |
d----- 9/14/2024 7:01 PM GSM4567892_Psoriasis06F | |
d----- 9/14/2024 7:01 PM GSM4567893_Psoriasis07 | |
d----- 9/14/2024 7:01 PM GSM4567894_Psoriasis08 | |
d----- 9/14/2024 7:01 PM GSM4567895_Psoriasis09 | |
d----- 9/14/2024 7:01 PM GSM4567896_Psoriasis10 | |
d----- 9/14/2024 7:01 PM GSM4567897_Psoriasis11 | |
d----- 9/14/2024 7:01 PM GSM4567898_Psoriasis12 | |
d----- 9/14/2024 7:01 PM GSM4567899_Psoriasis13 |
每个文件夹中的文件名都是一样的:
ls .\GSM4567877_Control01\ | |
目录: D:\Study\Project\GDG\scPsoriasis\data\GSE151177\GSE151177_RAW\GSM4567877_Control01 | |
Mode LastWriteTime Length Name | |
---- ------------- ------ ---- | |
-a---- 5/24/2020 10:59 AM 13776 barcodes.tsv.gz | |
-a---- 5/24/2020 10:59 AM 304728 features.tsv.gz | |
-a---- 5/24/2020 10:59 AM 16419935 matrix.mtx.gz |
在 Seurat V4 中,我们可以循环读取多个样本,然后 merge:
library(Seurat) | |
datadir <- "D:/Study/Project/GDG/scPsoriasis/data/GSE151177/" | |
setwd(file.path(datadir, "GSE151177_RAW")) | |
sample_dirs <- list.dirs(path = getwd(), full.names = FALSE, recursive = FALSE) | |
seurat_object_list <- lapply(sample_dirs, function(x) { | |
print(x) | |
folder <- file.path(datadair, "GSE151177_RAW/", x) | |
CreateSeuratObject(counts = Read10X(folder), | |
project = x, | |
min.cells = 3, | |
min.features = 100) | |
}) | |
# 合并多个样本 | |
seurat_object <- merge(seurat_object_list[[1]], | |
y = seurat_object_list[-1], | |
add.cell.ids = samples, | |
project = "GSE151177") |
但是 Seurat V5 修改了 Seurat 数据对象的结构,merge 操作并不会将样本的表达量矩阵合并
Read10X()
函数可以接受多个路径,因此在 Seurat V5 中,我们可以直接传入这些子文件夹:
library(Seurat) | |
datadir <- "D:/Study/Project/GDG/scPsoriasis/data/GSE151177/" | |
setwd(file.path(datadir, "GSE151177_RAW")) | |
sample_dirs <- list.dirs(path = getwd(), full.names = FALSE, recursive = FALSE) | |
seurat_object <- Read10X(sample_dirs) | |
seurat_object <- CreateSeuratObject(counts = seurat_object, min.cells = 3, min.features = 100) |
# txt/csv/tsv 数据
使用 data.table 的 fread()
函数读取
以 GSE150703 为例
library(data.table) | |
datadir <- "D:/Study/Project/GDG/scPsoriasis/data/GSE151177/" | |
GSE150703_dir <- file.path(datadir, "GSE150703") | |
GSE150703_file <- file.path(GSE150703_dir, "GSE150703_retina_NORM_OIR_P14_P17_C57_WR_CD73FT_noamg_normalizedUMI_Count_DGEmatrix.txt.gz") | |
data <- fread(GSE150703_file, data.table = F) |
# h5ad
H5AD 是一种基于 HDF5 (Hierarchical Data Format) 的高效文件格式,专为存储和分析单细胞转录组数据设计。通常 h5ad
文件由 Python 的单细胞分析工具 Scanpy 生成,但如果我们想在 R 语言中进行处理和绘图,也不是不可以。
首先需要安装这些包:
if (!requireNamespace("remotes", quietly = TRUE)) { | |
install.packages("remotes") | |
} | |
remotes::install_github("mojaveazure/seurat-disk") | |
install.packages("hdf5r") | |
install.packages("hdf5r.Extra") |
先将 h5ad
文件转化为 h5seurat
对象
h5ad_data <- file.path(datadir, "adata.h5ad") | |
Convert(h5ad_data, | |
dest = 'h5seurat', | |
overwrite = F, # 是否覆盖源 h5ad | |
assay = 'RNA') |
加载 h5seurat
对象,此时不要添加 metadata,因为通常会出错
# 使用 LoadH5Seurat 函数加载 h5seurat 文件,并创建 Seurat 对象 | |
seurat_object <- LoadH5Seurat(file.path(datadir, "adata.h5seurat"), meta.data = F, misc = F) |
提取 metadata 并添加
# 添加 metadata | |
meta_data <- h5Read(file.path(datadir, "adata.h5seurat"), "meta.data") | |
# 将每一列转换为因子 | |
meta_data <- as.data.frame(lapply(meta_data, factor)) | |
seurat_object@meta.data <- meta_data |
# h5
我目前没有遇到,看到其他教程中是使用 Read10X_h5()
函数即可
# 其他技巧
# 如何修改单细胞样本的名字?
我们提取之前的 GSE151177 数据中的 10 个样本(银屑病样本,试剂盒版本 V2.0,共 10 个样本)来示例:
psoriasis_v2_dirs <- c("GSM4567883_Psoriasis01", "GSM4567884_Psoriasis02", "GSM4567886_Psoriasis03", "GSM4567888_Psoriasis04", "GSM4567890_Psoriasis05", "GSM4567891_Psoriasis06", "GSM4567893_Psoriasis07", "GSM4567894_Psoriasis08", "GSM4567895_Psoriasis09", "GSM4567896_Psoriasis10") | |
seurat_object <- Read10X(psoriasis_v2_dirs) %>% | |
CreateSeuratObject(counts = ., min.cells = 3, min.features = 100) |
由此获得了一个 Seurat 对象 seurat_object,seurat_object 对象中生成的元数据字段 orig.ident
(通常存储在 meta.data
中)通常用于跟踪每个细胞样本来自哪个实验条件、患者、时间点或其他实验分组。这对于下游分析非常有用,尤其是当需要区分或比较不同样本时。可以用 seurat_object$orig.ident
来访问。
如果你是按照前文进行读取,你会发现 orig.ident
在 Seurat 对象中被存储为从 1 开始的 ** 因子类型(factor)** 而不是我们期望的样本名。在后续的作图分析中,如果样本信息只是数值,还要对应回原本的样本名去查看,这样并不直观。因此我们需要将 Seurat 对象中的 orig.ident
替换为真实的样本名字。
# 为什么替换会出错?
这里会出现一些顺序的问题,例如当我们的样本数量超过 9 个的时候:
> levels(seurat_object$orig.ident) | |
[1] "1" "10" "2" "3" "4" "5" "6" "7" "8" "9" |
可以发现合并样本的顺序是 "1" "10" "2"...
,在 R 中,因子的默认排序方式是按照字典序(也叫作字母顺序)进行排列的,而不是按照数值的大小进行排序。因为 R 默认将因子的水平(levels)视作字符,而在字典序中,"10" 会被认为比 "2" 小,因为比较的是字符串的每个字符的顺序:
- "1" 和 "10" 的首字符都是 "1",所以相同。
- 但在 "10" 中,第二个字符是 "0",而在 "2" 中是 "2",按照字母表顺序,"0" 比 "2" 更靠前。
如何知道这些因子对应的到底是哪些样本呢?可以通过查看细胞数量来判断:
> as.data.frame(table(seurat_object$orig.ident)) | |
Var1 Freq | |
1 1 660 | |
2 10 1976 | |
3 2 682 | |
4 3 1232 | |
5 4 945 | |
6 5 969 | |
7 6 1432 | |
8 7 768 | |
9 8 2186 | |
10 9 2891 | |
> # 查看每个 psoriasis_v2 样本的细胞数量 | |
> lapply(psoriasis_v2_dirs, function(x){ | |
+ print(x) | |
+ print(dim(CreateSeuratObject(counts = Read10X(x) , | |
+ min.cells = 3, | |
+ min.features = 100))) | |
+ }) | |
[1] "GSM4567883_Psoriasis01" | |
[1] 14913 660 | |
[1] "GSM4567884_Psoriasis02" | |
[1] 15027 682 | |
[1] "GSM4567886_Psoriasis03" | |
[1] 16899 1232 | |
[1] "GSM4567888_Psoriasis04" | |
[1] 15941 945 | |
[1] "GSM4567890_Psoriasis05" | |
[1] 15587 969 | |
[1] "GSM4567891_Psoriasis06" | |
[1] 16629 1432 | |
[1] "GSM4567893_Psoriasis07" | |
[1] 15836 768 | |
[1] "GSM4567894_Psoriasis08" | |
[1] 17437 2186 | |
[1] "GSM4567895_Psoriasis09" | |
[1] 18432 2891 | |
[1] "GSM4567896_Psoriasis10" | |
[1] 17813 1976 |
所以我们可以知道,在读取时, Read10X()
其实是按照我们的 psoriasis_v2_dirs
中的文件夹顺序来读取的,但是他并没有按照 psoriasis_v2_dirs
中的文件夹名去命名我们的样本,而是按顺序给了一个因子,只不过在合并样本的时候,由于因子的默认排序方式是按照字典序进行排列的,因此合并的顺序发生了变化。
知道对应关系后我们就可以替换样本名了,使用 Seurat 的 RenameIdents()
函数:
psoriasis_v2_samples <- c("Psoriasis01", "Psoriasis10", "Psoriasis02", "Psoriasis03", "Psoriasis04", "Psoriasis05", "Psoriasis06", "Psoriasis07", "Psoriasis08", "Psoriasis09") | |
names(psoriasis_v2_samples) <- levels(seurat_object$orig.ident) | |
seurat_object <- RenameIdents(seurat_object, psoriasis_v2_samples) |
# 更简单的方法
上面的思路确实可以解决问题,但还是太麻烦了,有没有更简单又强势的方法推荐一下呢?
有的兄弟,有的。其实在数据读入的时候就可以直接命名样本
通过查看 Read10X()
函数的源代码,我们可以发现其实这个函数有一个查看路径向量 data.dir
是否有命名的操作:
if (is.null(x = names(x = data.dir))) {...} |
所以我们可以直接在一开始使用样本名来命名文件夹向量:
psoriasis_v2_samples <- c("Psoriasis01", "Psoriasis02", "Psoriasis03", "Psoriasis04", "Psoriasis05", "Psoriasis06", "Psoriasis07", "Psoriasis08", "Psoriasis09", "Psoriasis10") | |
names(psoriasis_v2_dirs) <- psoriasis_v2_samples | |
seurat_object <- Read10X(psoriasis_v2_dirs) %>% | |
CreateSeuratObject(counts = ., min.cells = 3, min.features = 100) | |
> levels(seurat_object$orig.ident) | |
[1] "Psoriasis01" "Psoriasis02" "Psoriasis03" "Psoriasis04" "Psoriasis05" "Psoriasis06" "Psoriasis07" "Psoriasis08" "Psoriasis09" | |
[10] "Psoriasis10" |