學徒出師后去帝都某公司做單細胞服務了，最近他在解析一個食管癌單細胞文章，標題是：《Integrated single-cell transcriptome analysis reveals heterogeneity of esophageal squamous cell carcinoma microenvironment》，它并沒有給出來表達量矩陣，僅僅是一個fastq文件的公開。不過，他要跟我需要討論的是里面的的一個生存分析，如下所示：
原文的生存分析可以看到它這個生存分析是在0.05這個閾值邊緣瘋狂試探，但是學徒重復它這個分析的時候發(fā)現(xiàn)一定都不顯著，如下所示：
學徒重復發(fā)現(xiàn)不顯著我看到這樣的求助非常興奮，想著可以打假了。所以要來了數(shù)據(jù)和代碼，結果一個簡單的生存分析，學徒寫了兩百多行代碼，讓我頭疼。帶娃間隙勉強打開學徒的代碼試圖debug。
首先是下載芯片的表達量矩陣并且獲取CST1基因表達量
首先下載表達量矩陣
# Sys.setenv("VROOM_CONNECTION_SIZE"=999999)
# 參考：https://cran./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個探針， 它表達量矩陣里面的探針是純粹的數(shù)字，真的很讓人討厭。
所以取表達量矩陣，非常的注意！
然后整理臨床信息
這里需要自己去下載 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數(shù)據(jù)庫的表達量矩陣里面的gsm的id系統(tǒng)哦。
> 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
不過，勉強也算是有規(guī)律的。需要整理一下， 代碼如下所示：
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)
有了如下所示的數(shù)據(jù)：
> 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)
出圖如下所示：
確實是統(tǒng)計學顯著的文章里面是0.039，我做出來是0.041，差異并不大。
但是為什么學徒做出來都大于0.1了呢，我鼓起勇氣把他的200多行代碼看了看，發(fā)現(xiàn)其實問題出在他使用了 limma::normalizeBetweenArrays 這個函數(shù)對表達量矩陣進行簡單的處理。
一般來說，確實是有必要的，而且我們相信 limma::normalizeBetweenArrays 是有助于我們更好的尋找到真實的表達量上下調的差異基因。
但是這次居然如此巧合，僅僅是因為加上了 limma::normalizeBetweenArrays 就使得一個在文章里面有統(tǒng)計學顯著性的生存相關基因變得不顯著了。
大家怎么看這件事呢？使用 limma::normalizeBetweenArrays 這個函數(shù)對表達量矩陣進行簡單的處理有什么問題嗎？
還是說生存分析本來就是如此的玄乎呢？我在生信技能樹多次分享過生存分析的細節(jié)；
生存分析是目前腫瘤等疾病研究領域的點睛之筆！
