GEO数据作者不上传具体的样本对应信息是故意为之吗？

生信技能树

发布于 2025-05-09 10:40:42

2140

文章被收录于专栏：生信技能树生信技能树

收到技能树的生信马拉松学员的一个反馈，GEO中的一个公共数据如何找到对应的具体样本，如下：

老师们好，请教一下这个单细胞数据集的怎么读入呢？需要挑选具体的样本读入，给的GSE的那三个文件是整体的没有分组信息。GSM的RAW只给了barcodes的csv.

来看看这个数据吧！

此外，最新一期生信入门已经在5月5号开课，还没上车的可以看一看瞄一瞄：生信入门&数据挖掘线上直播课5月班，完全适合0基础小白。

数据背景

数据连接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE227691

去GEO数据库页面简单了解一下，总共有12个样本，三个分组，每组四个生物学重复：

没有COPD的吸烟者：吸烟但未发展为慢性阻塞性肺病的个体；
轻度COPD（COPD I级）：根据疾病严重程度分级，属于轻度慢性阻塞性肺病的患者；
中度COPD（COPD II级）：根据疾病严重程度分级，属于中度慢性阻塞性肺病的患者。

GSM7105661 wo COPD 1, scRNAseq
GSM7105662 wo COPD 2, scRNAseq
GSM7105663 wo COPD 3, scRNAseq
GSM7105664 wo COPD 4, scRNAseq
GSM7105665 mild COPD 1, scRNAseq
GSM7105666 mild COPD 2, scRNAseq
GSM7105667 mild COPD 3, scRNAseq
GSM7105668 mild COPD 4, scRNAseq
GSM7105669 moderate COPD 1, scRNAseq
GSM7105670 moderate COPD 2, scRNAseq
GSM7105671 moderate COPD 3, scRNAseq
GSM7105672 moderate COPD 4, scRNAseq

数据对应的文献已经发表了，于2023年5月5号发表在Front Immunol杂志上，标题为《Heme oxygenase-1 determines the cell fate of ferroptotic death of alveolar macrophages in COPD》。

作者提供的数据有下面的这些文件：

GSE227691_RAW.tar 160.0 Kb (http)(custom) TAR (of CSV)
GSE227691_barcodes.tsv.gz 29.2 Mb (ftp)(http) TSV
GSE227691_features.tsv.gz 311.3 Kb (ftp)(http) TSV
GSE227691_matrix.mtx.gz 225.1 Mb (ftp)(http) MTX

学员的问题：合并在一起的那三个文件读取进来，有没有办法区分不同的样本并与如GSM7105661样本编号对应起来，可以进行特定样本挑选分析？

数据读取看下

将上面的附件下载下来，整成下面的格式：

├── GSE227691
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz

读取：

###
### Create: Jianming Zeng
### Date:  2023-12-31  
### Email: jmzeng1314@163.com
### Blog: http://www.bio-info-trainee.com/
### Forum:  http://www.biotrainee.com/thread-1376-1-1.html
### CAFS/SUSTC/Eli Lilly/University of Macau
### Update Log: 2023-12-31   First version 
### Update Log: 2024-12-09   by juan zhang (492482942@qq.com)
### 

rm(list=ls())
options(stringsAsFactors = F)
library(ggsci)
library(dplyr) 
library(future)
library(Seurat)
library(clustree)
library(cowplot)
library(data.table)
library(ggplot2)
library(patchwork)
library(stringr)
library(qs)
library(Matrix)

# 创建目录
getwd()
gse <- "GSE227691"
dir.create(gse)

# 方式一：标准文件夹
###### step1: 导入数据 ######   
counts <- Read10X("./GSE227691/", gene.column = 2)
dim(counts)
sce.all <- CreateSeuratObject(counts, min.cells = 3, min.features = 200)
sce.all

# 查看特征
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:10])

读取进来后发现每个细胞barcode后面有不同的尾巴：

对这个尾巴进行切割就可以区分不同的样本：

head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
temp <- str_split(colnames(sce.all),n=2,pattern = "-",simplify = T)
head(temp)
table(temp[,2])
id <- as.data.frame(table(temp[,2]))
id

sce.all$orig.ident <- paste0("P",temp[,2])
table(sce.all$orig.ident)

这个方法确实可以区分不同样本，但是没办法知道哪个样本与具体的样本对应！

然后我看到了上面的每个样本单独的barcode，有希望！

├── GSM7105661_2018jz60_barcode_filtered.csv.gz
├── GSM7105662_2018jz47_barcode_filtered.csv.gz
├── GSM7105663_2018jz36_barcode_filtered.csv.gz
├── GSM7105664_C4_barcode_filtered.csv.gz
├── GSM7105665_2018jz13_barcode_filtered.csv.gz
├── GSM7105666_2018jz15_barcode_filtered.csv.gz
├── GSM7105667_2018jz61_barcode_filtered.csv.gz
├── GSM7105668_2018jz63_barcode_filtered.csv.gz
├── GSM7105669_2018jz14_barcode_filtered.csv.gz
├── GSM7105670_A4-SI-GA-A4-01_barcode_filtered.csv.gz
├── GSM7105671_2018jz39_barcode_filtered.csv.gz
├── GSM7105672_2018jz62_barcode_filtered.csv.gz

打开一看，希望破灭了：

每个样本的尾巴都是1，这跟上面的数据对应不上啊！

我又灵思一动：我可以对应细胞数啊！

统计作者给的每个样本的barcode数：

GSM7105661_2018jz60 5321
GSM7105662_2018jz47 1723
GSM7105663_2018jz36 1931
GSM7105664_C4 4188
GSM7105665_2018jz13 2937
GSM7105666_2018jz15 4334
GSM7105667_2018jz61 4091
GSM7105668_2018jz63 2051
GSM7105669_2018jz14 2890
GSM7105670_A4-SI-GA-A4-01 1865
GSM7105671_2018jz39 1360
GSM7105672_2018jz62 1228

统计代码：

ls *gz |perl -ne 'chomp;print "zless $_ | wc -l \n";' |sh >temp1
ls *gz |perl -ne 'chomp;/(.*)_barcode/;print "$1\n";' >temp2
paste temp2 temp1

前面读取的数据细胞数：

后面试了一些过滤参数，发现不能很好的对应上，也不能瞎猜！

去文献中找过滤参数，发现没有数据分析细节：

我又想了一个办法

做相关性计算，想起来之前老板就用生信分析方法找到了样本的分组信息：做生信压根就不怕造假或者出错

做一个伪bulk，然后计算不同样本的相关性：

# 数据经过了简单的过滤，此处省略其步骤
# 使用经典的AggregateExpression函数得到伪bulk 表达谱：
head(sce.all.filt@meta.data)

cor_data <- AggregateExpression(sce.all.filt, group.by = "orig.ident", 
                                assays = "RNA",return.seurat = TRUE)
cor_data <- cor_data[["RNA"]]$data
range(cor_data)
head(cor_data)

pheatmap::pheatmap(cor(cor_data))

dat <- cor_data
cg=names(tail(sort(apply(dat,1,sd)),1000))
exprSet=dat[cg,] # dat
pheatmap::pheatmap(cor(exprSet)) 
# 组内的样本的相似性应该是要高于组间的！
colD=data.frame(Group=group_list)
rownames(colD)=colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),
                   #annotation_col = colD,
                   show_rownames = F)
#～～～样品相关性热图p4～～～
p4 <- pheatmap::pheatmap(cor(exprSet),
                         # annotation_col = colD,
                         display_numbers = T,
                         show_rownames = T,
                         show_colnames =T)
p4

但是这也只能判断一下大致的分组：