# 银屑病单细胞数据 GSE151177 分析复现
文章标题为 Single-cell transcriptomics applied to emigrating cells from psoriasis elucidate pathogenic versus regulatory immune cell subsets（单细胞转录组学应用于银屑病迁移细胞阐明致病性与调节性免疫细胞亚群），发表在 Journal of Allergy And Clinical Immunology 上，IF14.2（2024）
数据为 GSE151177，文章提供了分析代码，但由于是 2021 年的文章，使用的是 Seurat V3，且代码冗余很多，因此我们使用 Seurat V5 进行复现。步骤和方法尽量保持和原文一致。
# 目录组织和数据下载
依据我的习惯，在项目目录 scPsoriasis 下新建如下目录，并将数据下载到 data/GSE151177/
目录结构
scPsoriasis/ | |
├─data/ | |
│ └─GSE151177/ | |
├─results/ | |
│ └─GSE151177/ | |
└─scripts/ | |
└─GSE151177/ |
解压后发现是 69 个 .gz 文件（23 个样本 * 3），文件名格式形如 GSM4567877_Control01_barcodes.tsv.gz
此外，在其 Github 仓库中还提供了一个表型文件 phenotype_data.txt
# 数据预处理
首先加载包，定义数据、结果和脚本路径，并将路径设置到结果
library(Seurat) | |
library(tidyverse) | |
rm(list = ls()) | |
datadir <- "D:/Study/Project/GDG/scPsoriasis/data/GSE151177/" | |
resultdir <- "D:/Study/Project/GDG/scPsoriasis/results/GSE151177/" | |
scriptdir <- "D:/Study/Project/GDG/scPsoriasis/scripts/GSE151177/" | |
setwd(resultdir) |
# 文件整理
该数据集是标准的 10X 测序输出文件，但我们仍需要先稍做处理
自定义 organize_files 函数，将文件按前缀分文件夹并去除多余的前缀：
# 将文件按照前缀移动到相应的子文件夹，并去除前缀的函数 | |
organize_files <- function(directory) { | |
# 设置工作目录为目标目录 | |
setwd(directory) | |
# 获取当前目录下所有文件 | |
files <- list.files(pattern = "*.gz") | |
# 遍历文件 | |
for (file in files) { | |
# 提取前缀（如 GSM4567877_Control01） | |
prefix <- sub("_(barcodes|features|matrix).*", "", file) | |
# 创建文件夹，检查文件夹是否存在，如果不存在则创建 | |
if (!dir.exists(prefix)) { | |
dir.create(prefix) | |
} | |
# 去掉前缀 | |
new_file <- sub("^.*?(barcodes|features|matrix).*", "\\1", file) | |
# 将文件移动到相应的文件夹并重命名 | |
file.rename(file, file.path(prefix, new_file)) | |
} | |
# 返回处理完成后的文件夹名称 | |
return(list.dirs(path = ".", full.names = FALSE, recursive = FALSE)) | |
} | |
sample_dirs <- organize_files(file.path(datadir, "GSE151177_RAW")) |
该函数会自动识别 .gz 文件并根据样本前缀创建相应的文件夹
执行完成后的 data 目录应如下：
data/ | |
└─GSE151177/ | |
└─GSE151177_RAW/ | |
├─GSM4567877_Control01/ | |
├─GSM4567878_Control02/ | |
├─GSM4567879_Control03/ | |
├─GSM4567880_Control04/ | |
├─GSM4567881_Control05/ | |
├─GSM4567882_Control05F/ | |
├─GSM4567883_Psoriasis01/ | |
├─GSM4567884_Psoriasis02/ | |
├─GSM4567885_Psoriasis02F/ | |
├─GSM4567886_Psoriasis03/ | |
├─GSM4567887_Psoriasis03F/ | |
├─GSM4567888_Psoriasis04/ | |
├─GSM4567889_Psoriasis04F/ | |
├─GSM4567890_Psoriasis05/ | |
├─GSM4567891_Psoriasis06/ | |
├─GSM4567892_Psoriasis06F/ | |
├─GSM4567893_Psoriasis07/ | |
├─GSM4567894_Psoriasis08/ | |
├─GSM4567895_Psoriasis09/ | |
├─GSM4567896_Psoriasis10/ | |
├─GSM4567897_Psoriasis11/ | |
├─GSM4567898_Psoriasis12/ | |
└─GSM4567899_Psoriasis13/ |
每一个文件夹下都是三个文件：
ls | |
目录: 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 |
# 样本分组
接下来，我们读取 phenotype_data.txt 文件，以获取样本的表型信息，帮助我们区分银屑病和对照样本：
phenotype <- read.delim("../phenotype_data.txt") |
根据表型信息，我们将样本分为三组：
- 银屑病（Psoriasis）_试剂盒版本 2.0
- 银屑病（Psoriasis）_试剂盒版本 3.0
- 对照（Control）_试剂盒版本 2.0
psoriasis_v2_samples <- phenotype$ID[phenotype$Category == "Psoriasis" & phenotype$Chemistry == "V2.0"] | |
psoriasis_v3_samples <- phenotype$ID[phenotype$Category == "Psoriasis" & phenotype$Chemistry == "V3.0"] | |
control_v2_samples <- phenotype$ID[phenotype$Category == "Control" & phenotype$Chemistry == "V2.0"] |
# 匹配文件夹
接下来，我们匹配样本目录，以便后续读取数据
自定义 get_phenotype_from_dir 和 match_sample_dirs 函数以根据 phenotype 的样本 ID 匹配对应的文件夹
# 从 sample_dirs 中提取 phenotype 部分的函数 | |
get_phenotype_from_dir <- function(dir) { | |
# 文件夹名的格式是 GSMxxx_phenotype，使用 strsplit 按照 "_" 分割 | |
phenotype_part <- strsplit(dir, "_")[[1]][2] | |
return(phenotype_part) | |
} | |
# 根据 phenotype 的样本 ID 匹配对应的文件夹名的函数 | |
match_sample_dirs <- function(samples, sample_dirs) { | |
# 提取所有文件夹名中的 phenotype 部分 | |
phenotype_dirs <- sapply(sample_dirs, get_phenotype_from_dir) | |
# 根据样本 ID 与 phenotype 部分进行精确匹配 | |
matched_dirs <- sample_dirs[phenotype_dirs %in% samples] | |
return(matched_dirs) | |
} |
匹配对应的文件夹
# 根据 phenotype 的信息动态构建分组 | |
psoriasis_v2_samples <- phenotype$ID[phenotype$Category == "Psoriasis" & phenotype$Chemistry == "V2.0"] | |
psoriasis_v3_samples <- phenotype$ID[phenotype$Category == "Psoriasis" & phenotype$Chemistry == "V3.0"] | |
control_v2_samples <- phenotype$ID[phenotype$Category == "Control" & phenotype$Chemistry == "V2.0"] | |
# 匹配对应的文件夹名 | |
psoriasis_v2_dirs <- match_sample_dirs(psoriasis_v2_samples, sample_dirs) | |
psoriasis_v3_dirs <- match_sample_dirs(psoriasis_v3_samples, sample_dirs) | |
control_v2_dirs <- match_sample_dirs(control_v2_samples, sample_dirs) |
命名文件夹向量，使合并后数据的 orig.ident 易识别
# 命名文件夹向量 | |
names(psoriasis_v2_dirs) <- psoriasis_v2_samples | |
names(psoriasis_v3_dirs) <- psoriasis_v3_samples | |
names(control_v2_dirs) <- control_v2_samples |
# 读取数据并创建 Seurat 对象
使用 Read10X 函数读取数据，并创建 Seurat 对象：
Psoriasis_V2.0 <- Read10X(psoriasis_v2_dirs) %>% | |
CreateSeuratObject(counts = ., min.cells = 3, min.features = 100) | |
Psoriasis_V3.0 <- Read10X(psoriasis_v3_dirs) %>% | |
CreateSeuratObject(counts = ., min.cells = 3, min.features = 100) | |
Control_V2.0 <- Read10X(control_v2_dirs) %>% | |
CreateSeuratObject(counts = ., min.cells = 3, min.features = 100) |
# 数据标准化
在合并样本之前，需要对每个 Seurat 对象进行标准化：
Psoriasis_V2.0 <- NormalizeData(Psoriasis_V2.0) %>% FindVariableFeatures(selection.method = "vst", nfeatures = 2000) | |
Psoriasis_V3.0 <- NormalizeData(Psoriasis_V3.0) %>% FindVariableFeatures(selection.method = "vst", nfeatures = 2000) | |
Control_V2.0 <- NormalizeData(Control_V2.0) %>% FindVariableFeatures(selection.method = "vst", nfeatures = 2000) |
# 数据整合
为了消除批次效应，使用 FindIntegrationAnchors() 函数识别不同数据集中的细胞锚点，然后通过 IntegrateData() 函数整合数据：
immune.anchors <- FindIntegrationAnchors(object.list = list(Psoriasis_V2.0, Psoriasis_V3.0, Control_V2.0), dims = 1:30) | |
Round0 <- IntegrateData(anchorset = immune.anchors, dims = 1:30) |
# 保存整合数据
最后，我们将整合后的数据保存为 Rda 文件，便于后续分析：
save(Round0, file = "Round0_integrated.Rda") |
