可以质疑来自哈佛医学院的顶刊Cell结果吗？

Duty
Bloger
Page content

可以质疑来自哈佛医学院的顶刊Cell结果吗？ by 生信技能树

公众号写作是一个逐步积累的过程，曾老板为了让我快速成长起来，每天都给我发了好多的资料，这不就在周一扔给了我一篇Cell文献，并迅速开了一个短会教我怎么提炼文章中的点。等我梳理完数据，今天要写的时候，脑子跟抽风了一样，一下子全忘光了。顶着被骂的风险，厚脸皮又去问了一遍。。。（下次一定长记性开完会就要及时记下笔记）。

最早看到这篇文章在某公众号：发现Cell文章“造假”，怎么办？？。他质疑了作者的单细胞数据分析结果有问题，下面来看看这篇 Cell文章中有什么古怪呢！！！

文献信息如下：

文献于2023年3月30号发表在顶刊杂志Cell上，标题为：《Cytotoxic CD4+ T cells eliminate senescent cells by targeting cytomegalovirus antigen》，文章的主要关键结论有以下几点：

  • 衰老成纤维细胞不随着年龄的增长而增加
  • 衰老皮肤中 CD4 CTL 与衰老成纤维细胞的积累呈负相关
  • 衰老的成纤维细胞高表达 HLA-II、NKG2D 配体和 HCMV-glycoprotein B
  • CD4 CTLs 以 HLA-II 依赖的方式消除 HCMV-gB+ 衰老成纤维细胞


通讯作者为来自马萨诸塞州总医院及哈佛医学院的 Shadmehr Demehri 教授：

https://www.cancerresearch.org/cri-funded-scientists/shadmehr-demehri-md-phd
https://www.cancerresearch.org/cri-funded-scientists/shadmehr-demehri-md-phd

质疑点一：bulk RNA-Seq 差异分析只有两个样本却使用了DESEq2分析？

文章中使用bulk转录组对年轻的皮肤和年老的皮肤进行测序，并使用DESEq2进行差异分析，检查了免疫相关基因在差异结果中的差异情况：CD69, CD96, CD276, CXCL9, KLRD1, IL2RB, IL17RC, IL36G, MB21D1, S100A9, and VTCN1, were among the genes significantly altered in the old versus young skin。

文献中的差异火山图如下：（如果是 1 vs 1的差异 出了差异分析方法有问题，火山图也不长这个样子，单样本的火山图只有一条v线）

image-20250228212350386

作者将数据上传到了GEO：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE191055

样本如下：只有两个样本是 皮肤组织 GSM5738155  和 GSM5738156

GSM5738155 Old skin
GSM5738156 Young skin
GSM5738157 P27_5: Senescent-induced human dermal fibroblast
GSM5738158 P27_6: Senescent-induced human dermal fibroblast
GSM5738159 P27_7: Senescent-induced human dermal fibroblast
GSM5738160 P27_8: Senescent-induced human dermal fibroblast
GSM5738161 P4_1: Non-senescent-induced human dermal fibroblast
GSM5738162 P4_2: Non-senescent-induced human dermal fibroblast
GSM5738163 P4_3: Non-senescent-induced human dermal fibroblast
GSM5738164 P4_4: Non-senescent-induced human dermal fibroblast

我们可以简单下载并进行处理：

rm(list = ls())#清空当前的工作环境
options(stringsAsFactors = F)#不以因子变量读取
options(scipen = 20)#不以科学计数法显示
library(data.table)
library(tinyarray)

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

# load counts table from GEO
#urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE224401&format=file&file=GSE224401_raw_counts_GRCh38.p13_NCBI.tsv.gz"
gse <- "GSE191055"
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&"
path <- paste0(urld, "acc=", gse, "&format=file&file=",gse,"_raw_counts_GRCh38.p13_NCBI.tsv.gz");
file <- paste0(gse,"_raw_counts_GRCh38.p13_NCBI.tsv.gz")
path
file
download.file(url = path, destfile = file)

tbl <- as.matrix(data.table::fread(file, header=T, colClasses="integer"), rownames=1)
dim(tbl)
ensembl_matrix=as.data.frame(tbl) 

library(AnnoProbe)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
e2s<-AnnotationDbi::select(org.Hs.eg.db,
                         keys= rownames(ensembl_matrix),
                         columns="SYMBOL",
                         keytype = "ENTREZID")
head(e2s)
ids = na.omit(e2s)
ids=ids[!duplicated(ids$SYMBOL),]
ids=ids[!duplicated(ids$ENTREZID),]
head(ids)
symbol_matrix= ensembl_matrix[match(ids$ENTREZID,rownames(ensembl_matrix)),]
rownames(symbol_matrix) = ids$SYMBOL
kp <- rowSums(symbol_matrix) > 0
table(kp)
symbol_matrix <- symbol_matrix[kp, ]
symbol_matrix[1:4,1:4]

library(AnnoProbe)
library(GEOquery)
gset = getGEO(gse, destdir = '.', getGPL = F) 
pd = pData(gset[[1]]) 
com <- intersect(rownames(pd), colnames(symbol_matrix))
symbol_matrix <- symbol_matrix[, com]
pd=pd[com,]
pd$title 
group_list= str_split(pd$title, pattern = "_", n=2, simplify = T)[,1]
group_list
table(group_list)
dat=log10(edgeR::cpm(symbol_matrix)+1)
save(symbol_matrix,dat,group_list,file = 'step1-output.Rdata')

那么DESEq2能进行单样本的差异分析吗？看官网：https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#can-i-use-deseq2-to-analyze-a-dataset-without-replicates

我们实际上是可以使用 edgeR来做单样本 vs  单样本的分析，见我们以前的稿子：没有生物学重复的转录组差异分析如何挑选基因呢：变化倍数与P值选谁？

疑问点二：明明是 4 vs 4的8个样本，结果展示中少了一个？

Fig3的热图：热图显示了衰老成纤维细胞与正常成纤维细胞中免疫细胞激活配体基因的表达情况（基于RNA测序数据）。蓝色表示下调（log2FC < -0.4，p < 0.05），红色表示上调（log2FC > 0.4，p < 0.05），对比的是衰老成纤维细胞（n = 3）与正常成纤维细胞（n = 4）。

这里少了一个衰老成纤维细胞（文章没有交代）：

如果简单的分析，使用三张图check一下，确实可以发现衰老成纤维P27组中有一个比较异常的样本：

疑问点三：作者还进行了单细胞数据分析，是整个文章的核心：CD4T CTL细胞

为了表征人类皮肤中的CD4细胞毒性T细胞（CD4 CTLs），作者对老年皮肤样本中的淋巴细胞进行了scRNA-seq（图S3A）。在CD4+ T细胞簇中，作者鉴定出一群表达穿孔素（PRF1）的CD4 CTLs，这些细胞高度上调了驻留记忆标记物CD69。此外，CD4 CTLs 还表达了与细胞毒性相关的基因，包括 RUNX3、RAB27A 和KLRK1（编码自然杀伤细胞2D受体[NKG2D]）。这些发现表明，CD4 CTLs 是清除老年皮肤中衰老细胞的效应细胞。

问：那究竟有没有这样的 突然冒出来的 一群CD4 CTL细胞，还高表达 PRF1 呢？

学徒作业

拿到这篇文章中的单细胞数据：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE221232，分析看看吧~ 欢迎留言交流

GSM6856506 old1
GSM6856507 old2
GSM6856508 old3

文末友情宣传


原文链接