在我们生信技能树的马拉松授课群里有个学员遇到一个有意思的事情：
在分析GEO芯片数据时，有两个GEO芯片数据，
实验设计一模一样，而且来自同一个课题组
，只有芯片平台不一样，但是对这两个数据做差异分析后，进行差异基因一致性比较，发现一致性很差。
下面就来看看~
来自同一个课题组的两个实验设计一模一样的数据
数据来自东京大学的外科肿瘤学系课题组，这两个数据分别为：
-
GSE3493
：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE3493
-
GSE35452
：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE35452
数据一：GSE3493
这个数据集包括了46例样本，35个 药物 non-responder 和 11 个 responder。
在直肠癌患者开始术前放疗之前，通过肠镜检查前瞻性地收集了活检标本。相应的肿瘤标本经福尔马林固定并包埋于石蜡中，用于组织学检查；其他标本则用于RNA提取。当相应的标本中至少含有70%的肿瘤细胞时（如参考文献9所述；见补充图S1），则使用这些标本进行RNA提取。标本在液氮中迅速冷冻，并在-80℃下保存，直至使用。所有患者接受了总计50.4Gy的放疗剂量，并在放疗结束后4周进行了标准化的根治性切除术。
数据芯片平台为：
GPL8300 [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array
数据二：GSE35452
与上面同样的实验设计：这个数据集包括46个样本，22 个 药物 non-responder 和 24 个 responder。
数据芯片平台为：
GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
简单的预处理数据并进行差异分析
第一个数据集的处理：GSE3493
rm(list=ls())
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
library(limma)
library(tidyverse)
## 1.获取并且检查表达量矩阵:gse编号需修改
gse_number <- "GSE3493"
dir.create(gse_number)
setwd(gse_number)
getwd()
#gset <- geoChina(gse_number)
gset <- getGEO(gse_number, destdir = '.', getGPL = T)
a <- gset[[1]]
## 2.样本分组:挑选一些感兴趣的临床表型。
pd <- pData(a)
pd$description <- gsub("test-","", pd$description)
pd$description <- gsub("training-","", pd$description)
group_list <- factor(pd$description, levels = c("non-responder","responder"))
table(group_list)
## 3.提取探针水平表达矩阵
dat <- exprs(a) # a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat) # 看一下dat这个矩阵的维度
dat[1:4,1:4] # 查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
## ～～～查看数据是否需要log～～～
range(dat)
# log2转换
dat <- log2(dat+1)
## 4.探针转换为基因symbol
gpl <- getGEO(filename=paste0(a@annotation,".soft.gz"))
gpl_anno <- gpl@dataTable@table
id2name <- gpl_anno[,c("ID" ,"Gene Symbol")]
colnames(id2name) <- c("ID","GENE_SYMBOL")
# 1.过滤掉空的探针
id2name <- na.omit(id2name)
id2name <- id2name[which(id2name$GENE_SYMBOL!=""), ]
# 2.过滤探针一对多
id2name <- id2name[!grepl("\\///",id2name$GENE_SYMBOL), ]
head(id2name)
# 3.多对一取均值
# 合并探针ID 与基因，表达谱对应关系
# 提取表达矩阵
dat <- dat %>%
as.data.frame() %>%
rownames_to_column("ID")
exp <- merge(id2name, dat, by.x="ID", by.y="ID")
# 多对一取均值
exp <- avereps(exp[,-c(1,2)],ID = exp$GENE_SYMBOL) %>%
as.data.frame()
dat <- as.matrix(exp[,pd$geo_accession])
dim(dat)
dat[1:5, 1:6]
save(gse_number, dat, group_list, pd, file = 'step1_output.Rdata')
差异分析为 Responder vs Non-responder
，差异结果如下，得到的差异基因比较少，我们前面也有遇到类似的情况：
有些差异本来就是不应该很明显
，
再小的差异也能被gsea找出来
。这其实也是合理的，因为病人药物治疗响应与不响应本身差异就非常小。
GSE35452 数据集下载CEL文件进行预处理：
## 读取芯片原始cell文件
library(oligo)
library(ggplot2)
# 下载raw文件夹
gse <- "GSE35452"
# https://ftp.ncbi.nlm.nih.gov/geo/series/GSE163nnn/GSE163558/suppl/GSE163558_RAW.tar
# 使用正则表达式替换：\w{3}$ 匹配每个单词最后的三个字符替换为空字符串，即去掉它们
s <- gsub("(\\w{3}$)", "", gse, perl = TRUE)
s
url <- paste0("https://ftp.ncbi.nlm.nih.gov/geo/series/",s,"nnn/",gse,"/suppl/",gse,"_RAW.tar")
url
# file <- paste0(gse,"_RAW.tar")
# file
# downloader::download(url, destfile = file)
# 把cel文件存放在文件夹里面
celFiles <- list.celfiles("GSE35452_RAW/",listGzipped=T, full.name=TRUE)
celFiles
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
exon_data <- oligo::read.celfiles(celFiles)
dim(exprs(exon_data))
exprs(exon_data)[1:4,1:4]
### 标准化(一步完成背景校正、分位数校正标准化和RMA(robust multichip average) 表达整合)
exon_data_rma <- oligo::rma(exon_data)
exp_probe <- Biobase::exprs(exon_data_rma)
exp_probe[1:4,1:4]
colnames(exp_probe) <- stringr::str_split_fixed(colnames(exp_probe), pattern = "_", n=2)[,1]
dat <- exp_probe
dim(dat)
dat[1:4,1:4]
后续与上面同样的代码然后进行
Responder vs Non-responder 差异分析
，得到的结果如下：差异同样的微弱
现在来比较两次差异分析的一致性如何
使用 FC 散点图比较两次差异结果：
rm(list = ls())
library(data.table)
# 读取数据
GSE3493 = fread(file = '../../01-GEO-database/2006-GSE3493-46-放疗响应与否直肠癌/GSE3493/DEG.csv',data.table = F)
GSE35452 = fread(file = '../../01-GEO-database/2013-GSE35452-46-放疗响应与否直肠癌/GSE35452//DEG.csv',data.table = F)
colnames(GSE35452)
rownames(GSE3493) = GSE3493$name
rownames(GSE35452) = GSE35452$name
head(GSE35452)
head(GSE3493)
ids=intersect(rownames(GSE3493),
rownames(GSE35452))
df= data.frame(
GSE3493 = GSE3493[ids,'logFC'],
GSE35452 = GSE35452[ids,'logFC']
)
rownames(df)=ids
library(ggpubr)
p <- ggscatter(df, x = "GSE3493", y = "GSE35452",
color = "black", shape = 21, size = 3, # Points color, shape and size
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.sep = "\n")
)
p
p <- p +
geom_hline(yintercept = -1, linetype = "dashed", color = "blue") +
geom_hline(yintercept = 0, linetype = "dotted", color = "red") +
geom_hline(yintercept = 1, linetype = "longdash", color = "green") +
geom_vline(xintercept = -1, linetype = "dashed", color = "blue") +
geom_vline(xintercept = 0, linetype = "dotted", color = "red") +
geom_vline(xintercept = 1, linetype = "longdash", color = "green")
p
ggsave(filename = 'cor-of-two-logFC.pdf', width = 6, height = 6)
可以看到两次差异分析结果基本上不相关：
这是为什么呢？
两个数据除了芯片平台不一致，其他的基本上都相同，我们来问问人工智能大模型看看能不能给我们一个合理的答案：
问：两个GEO芯片数据除了芯片平台不一致，实验设计一模一样，但是差异分析结果一致性却很差，可能的原因是什么呢？
两个GEO芯片数据在实验设计一致的情况下，差异分析结果一致性却很差，可能的原因包括以下几点：
1.
芯片平台差异
不同芯片平台的探针设计、检测原理和数据处理方式存在显著差异，这会导致基因表达水平的测量结果不一致。例如，Affymetrix、Illumina和Agilent等不同公司的芯片平台对探针强度的算法不同，不建议直接合并不同平台的数据。即使实验设计相同，平台差异也可能导致差异基因分析结果的不一致。
2.
批次效应（Batch Effect）
批次效应是指在实验过程中引入的非生物学的技术或实验变异，可能源自实验的不同运行批次、不同的实验人员、仪器的不同或者实验条件的变化等。批次效应可能导致数据不一致性，掩盖真实的生物学信号，增加噪声，从而影响差异分析结果的一致性。
3.
数据预处理和标准化
不同芯片数据的预处理和标准化方式可能不同。如果数据没有经过适当的标准化处理，或者标准化方法不一致，可能会导致差异分析结果的不一致。例如，某些数据可能已经进行了log2转换和标准化，而另一些数据则没有，这会影响后续的差异分析。
4.
样本数量和质量
样本数量较少时，数据的统计力不足，可能导致差异分析结果的不一致。此外，样本质量的差异（如RNA完整性、样本保存条件等）也会影响基因表达水平的测量结果。
