GEO数据库的多个表达量数据集的整合分析方法（表达量芯片和转录组测序）

# 导入数据集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()

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()