Step 1: Data Cleaning & Differential Analysis

###GSE33302
setwd("-")
library(data.table)
library(stringr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(dplyr)
library(tidyverse)
library(tinyarray)
library(limma)
library(GEOquery)
f='GSE33302_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE33302', destdir=".",
                 AnnotGPL = F,     
                 getGPL = F)       
  save(gset,file=f)   
}

load('GSE33302_eSet.Rdata')  


a=gset[[1]] 
dat=exprs(a) 
dim(dat) 
pd=pData(a)
pd=pd[,c(2,23)]
head(pd)
names(pd)<-c('sample','group')
pd$group=ifelse(pd$group=='non-sleep deprived mouse control','Control','SD')
group_list=pd[,2]
table(group_list)

####Probe Annotation
library(devtools)
#install_github("jmzeng1314/idmap1")
library(idmap1)
ids=getIDs('gpl1261')
head(ids)
ids=ids[,c(1,2)]
dat=as.data.frame(dat)
dat1=merge(ids,dat,by.x=1,by.y=0)
dat1=dat1[,-1]
dat1=dat1[!duplicated(dat1$symbol),]
rownames(dat1)<-dat1$symbol
dat1=dat1[,-1]
dat1=na.omit(dat1)
boxplot(dat1)
#dat1=log(dat1+1)
save(dat1,file = 'GSE33302.Rdata')
save(dat1,pd,group_list,file = 'GSE33302_dat_group.Rdata')

#######GSE9442
f='GSE9442_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE9442', destdir=".",
                 AnnotGPL = F,    
                 getGPL = F)       
  save(gset,file=f)   ## 保存到本地
}
load('GSE9442_eSet.Rdata')  ## 载入数据
a=gset[[1]] 
dat=exprs(a) 
dim(dat) 
pd=pData(a)
pd=pd[,c(2,1)]
head(pd)
names(pd)<-c('sample','group')
pd=pd[24:47,]
pd$group<-c(group=str_split(pd$group,'_',simplify = T)[,2])
pd$group=ifelse(pd$group=='SleepDep','SD',pd$group)
group_list=pd[,2]
table(group_list)


library(devtools)
#install_github("jmzeng1314/idmap1")
library(idmap1)
ids=getIDs('gpl1261')
head(ids)
ids=ids[,c(1,2)]
dat=as.data.frame(dat)
dat1=merge(ids,dat,by.x=1,by.y=0)
dat1=dat1[,-1]
dat1=dat1[!duplicated(dat1$symbol),]
rownames(dat1)<-dat1$symbol
dat1=dat1[,-1]
dat1=na.omit(dat1)
boxplot(dat1)
#dat1=log(dat1+1)
save(dat1,file = 'GSE9442.Rdata')
dat1=as.data.frame(t(dat1))
dat2=merge(pd,dat1,by.x='sample',by.y=0)
rownames(dat2)<-dat2$sample
dat2=dat2[,-c(1,2)]
dat2=as.data.frame(t(dat2))
save(dat2,pd,group_list,file = 'GSE9442_dat_group.Rdata')

########Data Merge
load("GSE33302_dat_group.Rdata")
exp1<-dat1
rm(dat1,pd,group_list)
load("GSE9442_dat_group.Rdata")
exp2=dat2
rm(dat2,pd,group_list)
exp1 <- exp1[intersect(rownames(exp1),rownames(exp2)),]
exp2 <- exp2[intersect(rownames(exp1),rownames(exp2)),]
boxplot(exp1)
boxplot(exp2)
exp = cbind(exp1,exp2)
boxplot(exp)
load("GSE33302_dat_group.Rdata")
group1=group_list
load("GSE9442_dat_group.Rdata")
group2=group_list
write.csv(pd,"pd.csv")
write.csv(pd1,"pd1.csv")
Group = c(group1,group2)
table(Group)
Group = factor(Group,levels = c("Control","SD"))
save(exp,Group,file = "SD_exp.Rdata")

#使用limma包里的removeBatchEffect()函数
rm(list = ls())
load("SD_exp.Rdata")
#处理批次效应
library(limma)
#?removeBatchEffect()
batch <- c(rep("Control",21),rep("SD",20))
exp2 <- removeBatchEffect(exp, batch)
par(mfrow=c(1,2)) 
boxplot(as.data.frame(exp),main="Original")
boxplot(as.data.frame(exp2),main="Batch corrected")


rm(list = ls())
library(sva)
library(limma)
group_list <- c(rep("Control",9),rep('SD',8),rep("Control",12),rep('SD',12))
gse<-c(rep('GSE33302',17),rep('GSE9442',24))
table(group_list,gse)
load("SD_exp.Rdata")
dat<-exp
dat[1:4,1:4]
betch<-c(rep('GSE33302',17),rep('GSE9442',24))
design=model.matrix(~group_list)
exp2<- removeBatchEffect(dat,batch = betch,design = design)
dim(exp2) 
boxplot(exp2)
exp2=normalizeBetweenArrays(exp2)
save(exp2,group_list,file = "SD_exp_combined.Rdata")


load("GSE33302_dat_group.Rdata")
pd1=pd
rm(exp,pd,group_list)
load("GSE9442_dat_group.Rdata")
pd2=pd
rm(exp,pd,group_list)
load("SD_exp.Rdata")
library(sva)
#?ComBat
batch <-data.frame(sample = c(pd1$sample,pd2$sample),
                   batch = c(pd1$group,pd2$group))
mod = model.matrix(~Group)
exp2 = ComBat(dat=exp, batch=batch$batch, 
              mod=NULL, par.prior=TRUE, ref.batch="Control")
par(mfrow=c(1,2))
library("FactoMineR")
library("factoextra")
pca.plot = function(dat,col){
  
  df.pca <- PCA(t(dat), graph = FALSE)
  fviz_pca_ind(df.pca,
               geom.ind = "point",
               col.ind = col ,
               addEllipses = TRUE,
               legend.title = "Groups"
  )
}
pca.plot(exp,factor(batch$batch))
pca.plot(exp2,factor(batch$batch))


####Difference Analysis
rm(list=ls())
load("SD_exp_combined.Rdata")
library(limma)
library(dplyr)
library(ggplot2)
design<-model.matrix(~0+factor(group_list1$Group))
colnames(design)<-levels(factor(group_list1$Group))
rownames(design)<-colnames(exp2)
contrast.matrix<-makeContrasts(SD-Control,levels = design ) 
fit<-lmFit(exp2,design)
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
DEG<-topTable(fit2,coef = 1,n=Inf)
DEG$regulate<-ifelse(DEG$adj.P.Val<=0.05,
                     ifelse(DEG$logFC>0, 'Up', 'Down'), 'Stable')
table(DEG$regulate)
write.csv(data.frame(gene_symbol=rownames(DEG),DEG),file = 'DEG.csv')
library(ggplot2)
pdf('DEG-volcano.pdf')
ggplot(DEG,aes(x=logFC,y=-log10(P.Value)))+
  geom_point(alpha=0.6,size=3.5,aes(color=regulate))+
  ylab('-log10(P.value)')+
  scale_color_manual(values=c('blue','grey','red'))+
  theme_bw()
while (!is.null(dev.list())) dev.off()

DEG<-read.csv('DEG.csv',sep=',',header = T)
DEG$regulate=ifelse(DEG$regulate=='Stable',NA,DEG$regulate)
DEG=na.omit(DEG)

#####Convert mouse and human gene IDs (R 4.3.2)
musGenes <- DEG
musGenes=musGenes[,2]
#BiocManager::install('biomaRt')
library(biomaRt)
#require("biomaRt")
human <- useMart('ensembl', dataset = 'hsapiens_gene_ensembl', host = 'https://dec2021.archive.ensembl.org/') #需要加host,不然报错
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl",host = 'https://dec2021.archive.ensembl.org/')

genes = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", 
               values = musGenes, 
               mart = mouse, 
               attributesL = c("hgnc_symbol"), 
               martL = human, uniqueRows=T)
genes=genes[,2]
write.csv(genes,file='SD-human.csv')


####Take the intersection of the two sets.
library(ggvenn)
library(ggplot2)
library(dplyr)
SD<-read.csv('SD-human.csv',sep=',',header = T)
autophagy<-read.csv('自噬相关基因集.csv',sep=',',header=T)
autophagy <- read.csv("自噬相关基因集.csv",sep=',',header=T)
datalist<-list('sleep deprivation'=SD$x,
               'autophagy'=autophagy$`Gene Symbol`)
#绘图
opar<-par(family='Roboto Condensed')
biocolor<-c('green','red')
p=ggvenn(datalist,
         fill_color = biocolor,
         fill_alpha = 0.5,
         stroke_linetype = 'longdash',
         set_name_size = 3,
         text_size = 4)
p
pdf(file='交集韦恩图.pdf',width = 8,height = 6)
print(p)
dev.off()
gene_list<-list(SD$x,autophagy$`Gene Symbol`)
common_gene<-Reduce(intersect,gene_list)
write.table(common_gene,file='gene.txt',sep='\t')