我们生信入门马拉松授课群里有个学员分析一个GEO数据时，发现这个数据很特别，他的提问：GSE65031这种作者只提供了合并数据的怎么办？来看下是怎么个合并情况吧！
数据背景
GSE65031 这个数据为小鼠肺组织样本，总共有4个分组为三个时间节点 + 一个对照，每个分组有5个样本，但是测序的时候呢同一个组别中的这五个样本的 RNA 进行了 pool，然后再测序。
具体描述如下：
Balb/cJ mice (5 per group) were intranasally challenged with 8-oxoguanine (1 µM, 60 µl) for three times at days 0, 2 and 4. Control group mice were intranasally challenged with saline (60 µl). At 30, 60 and 120 min after the third challenge (day 4), mice were sacrificed and lungs were processed for RNA extraction. RNAs from individual mice within the same group were pooled and subjected to deep-sequencing analysis in duplicate using NSG on an Illumina HiSeq 1000 sequencing system. After alignment and processing, the resulting RPKM from treatment groups (8-oxoG-challenged) were normalized to control group (saline-challenged).
数据链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65031
可以看到是个非常古老的测序仪了：GPL15103 Illumina HiSeq 1000 (Mus musculus)
GSM1586026_Sample5RPKM.txt.gz 198.4 Kb
GSM1586027_Sample6RPKM.txt.gz 199.7 Kb
GSM1586028_Sample7RPKM.txt.gz 200.5 Kb
GSM1586029_Sample8RPKM.txt.gz 200.6 Kb
这种情况，我们看下人工智能大模型kimi（https://kimi.moonshot.cn/）是否能给我们一些有用的建议呢，问：
有4个分组为三个时间节点 + 一个对照的高通量RNA测序数据，每个分组有5个生物学重复样本，但是测序的时候呢同一个组别中的这五个样本的 RNA 进行了 pool，然后再测序。这样的数据如何分析？
但是感觉 这个时候 kimi 不是很智能，他并没有给出一个很好的处理方式：
对于您提到的高通量RNA测序数据分析，其中包含三个时间节点和一个对照组，每个分组有5个生物学重复样本，且这些样本的RNA在测序前进行了pooling，以下是一些建议的分析步骤：
数据预处理：首先，需要对测序数据进行质量控制和标准化处理。这通常包括去除低质量的reads、接头序列以及进行比对到参考基因组等步骤。标准化是必要的，因为它有助于消除由于测序深度不同导致的表达量差异。 差异表达分析：在数据预处理之后，可以进行差异表达分析，以识别在不同时间点或对照组之间表达显著变化的基因。可以使用如DESeq2、edgeR或limma等R包来进行分析。这些工具能够处理重复样本，并考虑到生物学变异和测序技术变异。 数据可视化：为了更好地理解数据，可以进行主成分分析（PCA）或聚类分析，以可视化样本之间的相似性。此外，热图可以用来展示差异表达基因的表达模式。 功能富集分析：对差异表达的基因进行功能富集分析，如GO分析或KEGG通路分析，可以帮助理解这些基因在生物学过程中的作用。 考虑生物学重复：由于每个组别中的样本RNA进行了pooling，这可能会影响对生物学重复的分析。在分析时，需要考虑到这一点，并可能需要使用更复杂的统计模型来估计生物学变异。 数据解释和验证：最后，对分析结果进行解释，并可能需要通过实验（如qPCR）来验证关键基因的表达变化。 请注意，由于样本RNA进行了pooling，这可能会影响对单个样本之间差异的分析，因此在分析时需要特别注意这一点，并可能需要采用特定的统计方法来处理pooling的影响。
分析策略
其实，这种数据相对于使用差异分析来找不同时间节点的差异基因，我们可以使用更好的方法如 Mfuzz 分析。Mfuzz 分析可以找到 一些与时间变化相关的基因模块 ，且进行 Mfuzz 分析的时候，组内的生物学重复样本也会被处理成均值。
早期分享的关于mfuzz 分析的帖子：
简单的预处理一下数据，就可以使用上面的代码进行分析了：
rm(list = ls())#清空当前的工作环境
options(stringsAsFactors = F)#不以因子变量读取
options(scipen = 20)#不以科学计数法显示
library(data.table)
library(tinyarray)
library(GEOquery)
# 创建目录
acc <- "GSE65031"
dir.create(acc)
setwd(acc)
getwd()
list.files("./")
# 先读取一个样本看看
fs <- list.files('GSE65031_RAW//', full.names = T)
fs
tmp <- fread(fs[1],data.table = F)
head(tmp)
gid <- fread(fs[1],data.table = F)[,2]
head(gid)
# 批量读取并按照列合并
rpkm <- do.call(cbind, lapply(fs, function(x){
fread(x,data.table = F)
}))
head(rpkm)
identical(rpkm[,2], rpkm[,5])
grep("Sample",colnames(rpkm))
rpkm <- rpkm[,c(2,grep("Sample",colnames(rpkm)))]
rownames(rpkm) <- rpkm[,1]
rpkm <- rpkm[-1,] # 去掉第一行
rpkm <- rpkm[, -1] # 去掉第一列
# Klk7
rpkm['Klk7',]
head(rpkm)
# 过滤低表达
floor(0.75*ncol(rpkm))
keep <- rowSums(rpkm>0) >= floor(0.75*ncol(rpkm))
table(keep)
rpkm <- rpkm[keep,]
dim(rpkm)
head(rpkm)
获取样本的分组信息：
##############################
# 获取样本的分组信息
acc
gset <- getGEO(acc, destdir = '.', getGPL = F)
a <- gset[[1]]
a
## 2.样本分组
## 根据生物学背景及研究目的人为分组
## 挑选一些感兴趣的临床表型。
pd <- pData(a)
colnames(pd)
pd$title
group_list <- pd[,c("description","time:ch1")]
group_list$`time:ch1` <- gsub(" ","",group_list$`time:ch1`)
group_list
rownames(group_list) <- group_list$description
group_list
head(rpkm)
colnames(rpkm) <- group_list[colnames(rpkm),2]
head(rpkm)
dim(rpkm)
str(rpkm)
rpkm_exp <- matrix(as.numeric(as.matrix(rpkm)),nrow = nrow(rpkm))
rownames(rpkm_exp) <- rownames(rpkm)
colnames(rpkm_exp) <- colnames(rpkm)
head(rpkm_exp)
表达矩阵如下：
更多讲解，可以观看生信技能树的视频号：
学徒作业
拿到 GSE65031 数据集预处理后的表达矩阵，进行 Mfuzz 分析，并与 文章中的结果进行比较，文献《Whole transcriptome analysis reveals a role for OGG1-initiated DNA repair signaling in airway remodeling》：
