ixxmu
commented
6 months ago GEO数据库的多个表达量数据集的整合分析方法（表达量芯片和转录组测序） by 生信菜鸟团
在进行数据挖掘的时候，我们往往会筛选到不止一个符合我们预期的数据集，这些数据集来源于不同的研究人员。这样得到的这些数据集就会存在我们所谓的批次效应，如不同实验时间、不同实验批次、不同处理方法、不同测序平台等。遇到这一情况，我们该如何选择数据和处理数据呢？如果我们只选择一个数据集进行分析，貌似有点不太能说明其他研究人员的研究结果，但如果把所有符合我们实验目标的数据集都能拿来分析又有点无从下手。这里，我们就来介绍一下面对多个GEO数据集，我们该怎么处理？
首先，我们要明确一点，符合我们实验目标的数据集能搜集多少，尽可能的都用上，因为单独数据集的分析存在部分实验误差，不具有代表性。其次，针对多个数据集，我们可以有两种思路来进行整合分析：一是，合并和去除这些批次效应；二是，对各数据集分别进行处理，然后求交集，获得共有结果。
一、合并并去除批次效应
在GEO数据集合并和去除批次校正方法的方法主要包括ComBat方法（parametric prior method，ComBat_p和non-parametric method，ComBat_n）、代理变量法（Surrogate variable analysis，SVA）、基于比值的方法（Geometric ratio-based method，Ratio_G）、平均中心方法（Mean-centering，PAMR）和距离加权判别（Distance-weighted discrimination，DWD）方法。综合多种指标认为ComBat在精确性、准确性和整体性能方面（precision, accuracy and overall performance）总体优于其他方法。R中的SVA包中有ComBat和ComBat_seq函数可以用来校正批次效应，输入数据为干净的、标准化的表达数据（如FPKM、TPM等），通常是芯片数据。
代码示例如下
# 导入数据集1
data1 <- read.csv("PRJNA423456_FPKM.csv", header = TRUE) #这里的数据集为模拟数据集，大家改成自己的数据集进行分析
# 导入数据集2
data2 <- read.csv("PRJNA777728_FPKM.csv", header = TRUE) #这里的数据集为模拟数据集，大家改成自己的数据集进行分析
# 导入数据集3
data3 <- read.csv("PRJNA6878999_FPKM.csv", header = TRUE) #这里的数据集为模拟数据集，大家改成自己的数据集进行分析
# 导入数据集4
data4 <- read.csv("PRJNA590000_FPKM.csv", header = TRUE) #这里的数据集为模拟数据集，大家改成自己的数据集进行分析
#### 合并数据集
library("sva") #没有软件的话请先安装软件
library("tidyverse") #没有软件的话请先安装软件
merge_eset1 <- inner_join(data1, data2, by = "Gene.ID", suffix = c(".data1", ".data2")) #X是第一列名，需要将基因名放在第一列而不是成为行名
merge_eset1
merge_eset2 <- inner_join(data3, data4, by = "Gene.ID", suffix = c(".data3", ".data4"))
merge_eset <- inner_join(merge_eset1, merge_eset2, by = "Gene.ID", suffix = c(".merge_eset1", ".merge_eset2"))
merge_eset
rownames(merge_eset) <- merge_eset$Gene.ID
merge_eset <- merge_eset[,-1] #把第一列基因名放为行名
dim(merge_eset) #查看数据维度
exp <- as.matrix(merge_eset)
dimnames <- list(rownames(exp),colnames(exp))
data <- matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
dim(data)
class(data) #返回"matrix" "array"
write.csv(data,"mrna_nocombat.csv")
#过滤掉低表达的基因
Expr <- data[rowSums(data)>1,]
Expr=log2(Expr+1)
#PCA分析可视化个数据集在去批次前的聚簇效果
# 这里我们先获取分组信息
group_list <- data.frame(
sample = colnames(exp_all), c(rep("PRJNAX423456_A",3),rep("PRJNA423456_B",3),rep("PRJNA777728_A",3),rep("PRJNA777728_C",3),
rep("PRJNA777728_B",3),rep("PRJNA6878999_A",10),rep("PRJNA6878999_B",9),rep("PRJNA590000_A",3),rep("PRJNA590000_C",3),
rep("PRJNA590000_B",3)))
rownames(group_list) <- group_list$sample
colnames(group_list)[2] <- "dataset"
group <- factor(group_list$dataset)
View(group_list)
##载入样品信息，含批次
data2 <- read.csv("All_data.csv",header = T)
data2[,2] <- as.factor(data2$Type )##type表示的是生物学上的不同处理
data2[,3] <- as.factor(data2$Batch )##batch表示不同的批次信息
row.names(data2) <- data2$Sample
View(data2)
data2
library(FactoMineR)##没有请先安装
library(factoextra)
pdf(file = "PCA_before1.pdf",width = 7,height = 6)
pre.pca <- PCA(t(exp_all),graph = FALSE)
fviz_pca_ind(pre.pca,
geom= "point",
col.ind = data2$Batch,
addEllipses = TRUE,
legend.title="Group"
)
dev.off()
#-------------------进行去批次------------------#
#### 使用sva包计算批次效应
library(sva)
exp_all_combat <- ComBat(exp_all, batch = group_list$dataset) # batch为批次信息
# 查看去除批次效应的结果如何
library(tinyarray)
draw_pca(exp = exp_all_combat, group_list = group)
#做校正的PCA分析
pdf(file = "9_PCA_after.pdf",width = 7,height = 6)
pre.pca <- PCA(t(exp_all_combat),graph = FALSE)
fviz_pca_ind(pre.pca,
geom= "point",
col.ind = data2$Batch,
addEllipses = TRUE,
legend.title="Group")
dev.off()
值得注意的是，去批次效应处理后，个别基因的表达量会出现负数，这个时候用我们常规的edgeR和DEseq2进行差异分析是行不通的。而且上述我们使用FPKM值做的去除批次效应，所以在进行差异分析时，不可使用edgeR和DEseq2，但是可以用limma包进行差异分析。对于WGCNA的分析，似乎不受影响。
二、整合数据及分析
在数据挖掘过程中，我们同时会分析多个数据集的表达谱数据，这样就会都得到多个差异分析列表。那么，怎么样才能挑出一些更重要的或者更有生物学意义的基因进行后续实验呢？常规做法就是将三个数据集的差异基因列表进行overlapping，但这种方法只考虑到了gene出现的次数，并没有考虑到基因在多个差异分列表中排序上的重要性。
Robust Rank Aggregation(RRA)方法则可以对对多个排好序的基因集进行求交集，同时还考虑一下它们的排序情况。总体上来说，就是挑选那些在多个数据集都表现差异的基因，并且每次差异都排名靠前的那些，他们的最终综合排名也会比较靠前。
代码示例如下
setwd("F:\\RRA分析-2024-03-29\\4_RRA分析-2\\Anagen_Catagen")
#1.读入GSE数据集
files=c("PRJNA590000_A_C.DESeq2.csv",
"PRJNA779755_A_C.DESeq2.csv")
upList=list()
downList=list()
allFCList=list()
for(i in 1:length(files)){
inputFile=files[i]
rt=read.csv(inputFile)
header=unlist(strsplit(inputFile,"_"))
downList[[header[1]]]=as.vector(rt[,1])
upList[[header[1]]]=rev(as.vector(rt[,1]))
fcCol=rt[,1:2] colnames(fcCol)=c("Gene",header[[1]])
allFCList[[header[1]]]=fcCol
}
View(allFCList)
#2.合并所有的FC值
mergeLe=function(x,y){
merge(x,y,by="Gene",all=FALSE)}
newTab=Reduce(mergeLe,allFCList)
View(newTab)
rownames(newTab)=newTab[,1]
newTab=newTab[,2:ncol(newTab)]
newTab[is.na(newTab)]=0
dim(newTab)
#3.获取上调基因
library(RobustRankAggreg)
upMatrix = rankMatrix(upList)
View(upMatrix)
dim(upMatrix)
upAR = aggregateRanks(rmat=upMatrix)
View(upAR)
dim(upAR)
colnames(upAR)=c("Name","Pvalue")
upAdj=p.adjust(upAR$Pvalue,method="bonferroni")
#help(p.adjust)
upXls=cbind(upAR,adjPvalue=upAdj)
upFC=newTab[as.vector(upXls[,1]),]
upXls=cbind(upXls,logFC=rowMeans(upFC))
write.table(upXls,file="up.xls",sep="\t",quote=F,row.names=F)
upSig=upXls[(upXls$Pvalue<0.01 & upXls$logFC> log2(2)),]
dim(upSig)
write.table(upSig,file="upSig_2.xls",sep="\t",quote=F,row.names=F)
#3.获取下调基因
downMatrix = rankMatrix(downList)
downAR = aggregateRanks(rmat=downMatrix)
colnames(downAR)=c("Name","Pvalue")
downAdj=p.adjust(downAR$Pvalue,method="BH")
downXls=cbind(downAR,adjPvalue=downAdj)
downFC=newTab[as.vector(downXls[,1]),]
downXls=cbind(downXls,logFC=rowMeans(downFC))
write.table(downXls,file="down.xls",sep="\t",quote=F,row.names=F)
downSig=downXls[(downXls$Pvalue<0.01 & downXls$logFC< -log2(2)),]
dim(downSig)
write.table(downSig,file="downSig_2.xls",sep="\t",quote=F,row.names=F)
#整合上下调基因
allSig = rbind(upSig,downSig)
colnames(allSig)
allSig = allSig[,c("Name","logFC")]
View(allSig)
write.table(allSig,file = 'allSign_2.xls',sep = '\t',quote = F)
#显示排名靠前的前20个上下调基因的结果logFC.tiff
hminput=newTab[c(as.vector(upSig[1:20,1]),as.vector(downSig[c(1:20),1])),]
library(pheatmap)
tiff(file="logFC2_1_2.tiff",width = 15,height = 20,units ="cm",compression="lzw",bg="white",res=400)
pheatmap(hminput,display_numbers = TRUE,fontsize_row=10,fontsize_col=12,
color = colorRampPalette(c("blue","white", "red"))(200),
cluster_cols = FALSE,cluster_rows = FALSE, angle_col = 45)
dev.off()
目前最常用的整合分析方法即为以上所介绍的两种方法，最后给大家放置几篇相关参考文献以供借鉴！
Xing J, Chen M, Han Y. Multiple datasets to explore the tumor microenvironment of cutaneous squamous cell carcinoma. Math Biosci Eng. 2022 Apr 8;19(6):5905-5924. doi: 10.3934/mbe.2022276. IF: 2.6 Q2. PMID: 35603384. https://pubmed.ncbi.nlm.nih.gov/35603384/
Li Z, Li X, Jin M, Liu Y, He Y, Jia N, Cui X, Liu Y, Hu G, Yu Q. Identification of potential biomarkers and their correlation with immune infiltration cells in schizophrenia using combinative bioinformatics strategy. Psychiatry Res. 2022 Aug;314:114658. doi: 10.1016/j.psychres.2022.114658. IF: 11.3 Q1. https://pubmed.ncbi.nlm.nih.gov/35660966/
Mo M, Liu B, Luo Y, Tan JHJ, Zeng X, Zeng X, Huang D, Li C, Liu S, Qiu X. Construction and Comprehensive Analysis of a circRNA-miRNA-mRNA Regulatory Network to Reveal the Pathogenesis of Hepatocellular Carcinoma. Front Mol Biosci. 2022 Jan 24;9:801478. doi: 10.3389/fmolb.2022.801478. IF: 5.0 Q2. https://pubmed.ncbi.nlm.nih.gov/35141281/
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
生物信息学马拉松授课（买一得五） ，你的生物信息学入门课 时隔5年，我们的生信技能树VIP学徒继续招生啦 2024的共享服务器交个朋友福利价仍然是800 千呼万唤始出来的独享生物信息学云服务器 永久免费的生信进阶培训（线下）
https://mp.weixin.qq.com/s/T9ITQwUrZpq4UyTQwMWeiA