library(GEOquery)
library(limma)
library(ggplot2)
library(ggrepel)
library(ggthemes)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

getGEOdata <- function(x, GSE.id=FALSE){
  # Use the getGEO() function to download the GEO data for the id stored in x
  GSEdata <- getGEO(x, GSEMatrix=T, AnnotGPL=FALSE)
  # check the data: printing a summary of the expression values for the first 2 columns
  print(summary(exprs(GSEdata[[1]])[,1:2]))
  # Get the eset object
  eset <- GSEdata[[1]]
  # Save the objects generated for future use in the current working directory
  save(GSEdata, eset, file=paste(x, ".RData", sep=""))
  # check whether we want to return the list object we downloaded on GEO or
  # just the eset object with the getGSEobject argument
  if(GSE.id) return(GSEdata) else return(eset)
}

# Store the dataset ids in a vector GEO_dataset just in case you want to loop through several GEO ids
GEO_dataset <- c("GSE17215") #芯片数据在GEO数据库里的ID

# 获得表达矩阵
eset <- getGEOdata(GEO_dataset[1])
exprSet <- as.data.frame(exprs(eset))
# 提取GSM ID作为样品名
GSMID <- pData(eset)[,c(1,2,8,10,12)]
colnames(exprSet) <- GSMID$geo_accession

# 获得探針组跟gene symbol的对应关系
gpl <- getGEO(eset$platform_id[1])
dim(Table(gpl)) #查看一共有多少行探针组

# 到此為止，所需要的表达矩阵和芯片平台注释已經下載成功

colnames(Table(gpl)) #查看列名
#head(Table(gpl))

# 提取自己想要的信息，此处第1、11、12列，分别是行名ID和其对应的Gene Symbol和ENTREZ_GENE_ID
annotation <- Table(gpl)[,c(1,11,12)]
i <- sapply(annotation, is.factor) # Change factor to character
annotation[i] <- lapply(annotation[i], as.character)

# 有些探针组对应多個基因，使用“///”分割，當然你可以選擇保留，我的習慣是保留第一個，所以這裡對於其他的冗余信息進行clean
annotation[,2] <- data.frame(sapply(annotation[,2], function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]
annotation[,3] <- data.frame(sapply(annotation[,3], function(x) unlist(strsplit(x,'///'))[1]),stringsAsFactors = F)[,1]
# 或者宁缺毋滥，删除对应多个基因的探针组，用下面两行代替上面两行
#annotation <- annotation[!grepl("///", annotation$`Gene Symbol`),]
#dim(annotation)

# 把整理好的探针组ID、gene symbol、ENTREZ ID对应关系保存到文件里
colnames(annotation) <- c("probe_id", "gsym", "entrez_id")
write.csv(annotation, "gplTOgene.csv", row.names = F)

annotation <- read.csv("gplTOgene.csv")
head(annotation)

exprSet$probe_id <- rownames(exprSet)
express <- merge(x = exprSet, y = annotation[,c(1:2)], by = "probe_id", all.x = T) 
#如果想轉換成ENTREZ_GENE_ID，就用下面这行
#express <- merge(x = exprSet, y = annotation[,c(1:3)], by = "probe_id", all.x = T) 
express$probe_id <- NULL

# 有些gene对应多个探针组，此处保留其中表达量最大值
dim(express)
rowMeans <- apply(express, 1, function(x) mean(as.numeric(x), na.rm = T))
express <- express[order(rowMeans, decreasing =T),]
express <- express[!duplicated(express[, dim(express)[2]]),]  #express的最后一列为gene name
express <- na.omit(express)
dim(express)

rownames(express) <- express[,dim(express)[2]]
express <- express[,-(dim(express)[2])]
head(express)
# 至此，表达矩阵（express）已构建好，保存到文件
write.csv(express, "easy_input_expr.csv", quote = F)

# 获得样品信息
pdata <- pData(eset)
i <- sapply(pdata, is.factor)
pdata[i] <- lapply(pdata[i], as.character)

# 按照实验设计，构建分组矩阵
group_list <- subset(pdata, select = title) # Sample的分组信息
group_list$condition = rep(c("a","b"), each = 3)
group_list

# 至此，分组矩阵（design）已构建好
write.table(group_list, "easy_input_pheno.txt", quote = F, sep = "\t", row.names = F)

express <- read.csv("easy_input_expr.csv", row.names = 1)

# Have to check whether this matrix is log2 based
# usually, log2 value should be in between 0 - 18
range(express)
# if yes no need to log transfer, if above this range, have to do log transfer.
express <- log2(express)

par(mfrow=c(2,1))
boxplot(express,outline = F) #均一化前

# Normalise data quantiles
library(preprocessCore)
express.norm <- normalize.quantiles(as.matrix(express))

colnames(express.norm) <- colnames(express)
rownames(express.norm) <- rownames(express)
boxplot(express.norm, outline=F) #均一化后

#不需要再做均一化处理，就只运行这一行：
#express.norm <- express

group_list <- read.table("easy_input_pheno.txt", sep = "\t", header = T)
group_list

# Correlation between the samples
heatmap(cor(express.norm), scale = "none")
round(cor(express.norm), 2)

# 通過correlation，可能会發現某个sample与同一組其他sample之间距離很遠，需要删除。
# 此处删除GSM431123
group_list <- group_list[c(-3),] #删除第3个sample
express.norm <- express.norm[,c(-3)] #删除第3列

# Correlation between the samples again
heatmap(cor(express.norm),scale="none")
round(cor(express.norm),2)

design <- model.matrix(~ 0 + factor(group_list$condition))
colnames(design) <- levels(factor(group_list$condition))
rownames(design) <- colnames(express.norm)
design

# 构建差异比较矩阵
contrast.matrix <- makeContrasts(b-a, levels = design)
# 至此，差异表达矩阵已构建好

fit <- lmFit(express.norm,design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2)
# 得到两两差异表达的结果
# b vs. a
x <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")

#把全部基因的limma分析结果保存到文件
write.csv(x, "easy_input_limma.csv", quote = F)