# 单细胞转录组数据分析入门 - 2 - 各种格式的单细胞数据读取

单细胞转录组测序的格式可能有多种,包括单样本,多样本,10X 的标准输出文件,h5,h5ad,txt/csv/tsv

网上的教程鱼龙混杂,这使我在初次接触时晕头转向

# 10X 标准输出

这是目前比较主流的一种数据格式,对于一个样本,其包含:

  • barcodes.tsv.gz :记录 cell id
  • features.tsv.gz :记录 gene id
  • matrix.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"