哪有什么对错呢 - Githubissues

哪有什么对错呢 by 生信技能树

学徒出师后去帝都某公司做单细胞服务了，最近他在解析一个食管癌单细胞文章，标题是：《Integrated single-cell transcriptome analysis reveals heterogeneity of esophageal squamous cell carcinoma microenvironment》，它并没有给出来表达量矩阵，仅仅是一个fastq文件的公开。不过，他要跟我需要讨论的是里面的的一个生存分析，如下所示：

原文的生存分析

可以看到它这个生存分析是在0.05这个阈值边缘疯狂试探，但是学徒重复它这个分析的时候发现一定都不显著，如下所示：

学徒重复发现不显著

我看到这样的求助非常兴奋，想着可以打假了。所以要来了数据和代码，结果一个简单的生存分析，学徒写了两百多行代码，让我头疼。带娃间隙勉强打开学徒的代码试图debug。

首先是下载芯片的表达量矩阵并且获取CST1基因表达量

首先下载表达量矩阵

# Sys.setenv("VROOM_CONNECTION_SIZE"=999999)
# 参考：https://cran.r-project.org/web/packages/vroom/readme/README.html 
library(AnnoProbe)
library(GEOquery) 
dir.create("GSE53624",recursive = T)
eset <- getGEO("GSE53624", destdir = "./GSE53624", getGPL = F)
expr <- as.data.frame(exprs(eset[[1]]))
head(expr[,1:4])

可以看到，这个探针蛮奇怪的，居然是顺序编号：

> head(expr[,1:4])
   GSM1297076 GSM1297077 GSM1297078 GSM1297079
1   14.196541  14.151714  13.796948  13.802610
2    3.195847   3.042514   3.211573   2.995495
24  15.261637  15.830739  15.311610  15.527160
25   3.157660   3.378073   3.165554   3.634285
26   5.277165   5.297271   5.193853   5.467351
27   8.545228   8.327291   8.527834   8.590668

接着获取探针的注释信息：

library(AnnoProbe)
# 前面的gse里面有提到它芯片对应的gpl平台哦
ids=idmap('GPL18109','pipe')
head(ids)

很容易找到自己的感兴趣的基因：

> ids[ids$symbol=='CST1',]
      probe_id symbol
80298    69686   CST1

CST1 = expr[as.character(ids[ids$symbol=='CST1',1]),]

可以看到是 69686 这个探针，注意哦，并不是说是第69686个探针，其实是第80298个探针，它表达量矩阵里面的探针是纯粹的数字，真的很让人讨厌。

所以取表达量矩阵，非常的注意！

然后整理临床信息

这里需要自己去下载 GSE53624_clinical_data_of_patients_orignial_set.xls这个Excel文件哦。

clinical <- xlsx::read.xlsx("./GSE53624_clinical_data_of_patients_orignial_set.xlsx",sheetIndex = 1) 
table(clinical$Death.at.FU) 
clinical$Death.at.FU <- gsub("no","0",
                             gsub("yes","1",clinical$Death.at.FU)) 
clinical_data <- data.frame(OS.time=as.numeric(clinical$Survival.time.months.),
                            OS=as.numeric(clinical$Death.at.FU),
                            sample=clinical$Patient.ID)
head(clinical_data)

可以看到这个时候的临床信息里面的每个病人的id并不是我们下载的geo数据库的表达量矩阵里面的gsm的id系统哦。

> head(clinical_data)
    OS.time OS sample
1 48.766667  1    ec4
2  9.766667  1    ec6
3  5.833333  1    ec7
4 72.533333  0    ec9
5 72.633333  0   ec10
6 35.033333  1   ec11

不过，勉强也算是有规律的。需要整理一下，代码如下所示：

phenotype <- pData(eset[[1]])
phe1 <- data.frame(sample = rownames(phenotype),
                   title = phenotype$title) 
phe1$tissue <- str_split(phe1$title," ",simplify = T)[,1]
phe1$patient <- str_split(phe1$title," ",simplify = T)[,5]
head(phe1)
phe1=phe1[phe1$tissue == 'cancer',]
phe1$patient=paste0('ec',phe1$patient)
identical(phe1$patient,clinical_data$sample)

clinical_data=clinical_data[match(phe1$patient,clinical_data$sample),]
CST1=CST1[match(phe1$sample,colnames(expr))]

library(survival)
CST1=as.numeric(CST1)
clinical_data$gene = ifelse( CST1  > median( CST1 ),'high','low')
head(clinical_data)

有了如下所示的数据：

> head(clinical_data)
    OS.time OS sample gene
65 60.30000  0  ec224  low
66 27.56667  1  ec225 high
67 34.66667  1  ec226 high
68 60.96667  0  ec227  low
81 15.43333  1  ec251 high
82 61.33333  0  ec253 high

接下来就是纯粹的生存分析啦：


sfit1=survfit(Surv(OS.time, OS)~gene, data=clinical_data)
ggsurvplot(sfit1,pval =TRUE, data = clinical_data, risk.table = TRUE)

出图如下所示：

确实是统计学显著的

文章里面是0.039，我做出来是0.041，差异并不大。

但是为什么学徒做出来都大于0.1了呢，我鼓起勇气把他的200多行代码看了看，发现其实问题出在他使用了 limma::normalizeBetweenArrays 这个函数对表达量矩阵进行简单的处理。

一般来说，确实是有必要的，而且我们相信 limma::normalizeBetweenArrays 是有助于我们更好的寻找到真实的表达量上下调的差异基因。

但是这次居然如此巧合，仅仅是因为加上了 limma::normalizeBetweenArrays 就使得一个在文章里面有统计学显著性的生存相关基因变得不显著了。

大家怎么看这件事呢？使用 limma::normalizeBetweenArrays 这个函数对表达量矩阵进行简单的处理有什么问题吗？

还是说生存分析本来就是如此的玄乎呢？我在生信技能树多次分享过生存分析的细节；

生存分析是目前肿瘤等疾病研究领域的点睛之笔！

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
生信入门课-2022年6月场，你的生物信息学第一课

ixxmu / mp_duty

哪有什么对错呢 #2227

哪有什么对错呢 by 生信技能树

首先是下载芯片的表达量矩阵并且获取CST1基因表达量

然后整理临床信息

文末友情宣传