学习GEO

exprSet_t <- t(exprSet)
as.data.frame()

group <- rep(c("A","B"),each=3)
groupdata <- cbind(group=group,dataframe)

data <- read.table("data/GSE9103_series_matrix.txt",sep = "\t",comment.char="!",stringsAsFactors=F,header=T)

rownames(data) <- data[,1] #命名行名
data <- data[,-1]  #删除第一列

ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
## 开始判断
if (LogC) {
ex[which(ex <= 0)] <- NaN
## 取log2
exprSet <- log2(ex)
print("log2 transform finished")
}else{
  print("log2 transform not needed")
}

library(limma)
boxplot(data1,outline=FALSE, notch=T, las=2)
### 该函数默认使用quntile 矫正差异
data1=normalizeBetweenArrays(data1)
boxplot(data1,outline=FALSE, notch=T, las=2)
## 这步把矩阵转换为数据框很重要
data1 <- as.data.frame(data1)

platformMap <- data.table::fread("platformMap.txt",data.table = F)
index <- "GPL570"
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")
if(!requireNamespace("hgu133plus2.db")){
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
BiocManager::install("hgu133plus2.db",update = F,ask = F)
}

library(hgu133plus2.db)
bbq <- toTable(get("hgu133plus2SYMBOL"))
## 探针有多少个？
length(unique(a1$probe_id))
## 这么多行中，基因名称有重复的么？
length(unique(a1$symbol))

library(dplyr)
library(tibble)

data2 <- data1 %>%
## 行名转列名,因为只有变成数据框的列,才可以用inner_join
rownames_to_column("probe_id") %>%
## 合并探针的信息
inner_join(bbq,by="probe_id") %>%
## 去掉多余信息
select(-probe_id) %>%  
## 重新排列
select(symbol,everything()) %>%  
## rowMeans求出行的平均数(这边的.代表上面传入的数据)
## .[,-1]表示去掉出入数据的第一列，然后求行的平均值
mutate(rowMean =rowMeans(.[,-1])) %>%
## 把表达量的平均值按从大到小排序
arrange(desc(rowMean)) %>%
## 去重，symbol留下第一个
distinct(symbol,.keep_all = T) %>%
## 反向选择去除rowMean这一列
select(-rowMean) %>%
## 列名转行名
column_to_rownames("symbol")
### 保存数据
save(exprSet,file = "output/exprSet_rmdup.Rdata")

library(limma)
load(file = "output/exprSet_rmdup.Rdata")
### 创建分组;这一步根据样本来就行，原则就是: 跟样本匹配，取决于样本的排序
group <- c(rep("con",3),rep("treat",3))
##group <- c("con","con","treat","con","treat","treat")
### 分组变成向量，并且限定leves的顺序;levels里面，把对照组放在前面
group <- factor(group,levels = c("con","treat"),ordered = F)

### 主成分分析PCA：提前预测结果
### 行是样本列是基因
res.pca <- prcomp(t(exprSet), scale = TRUE)
library(factoextra)
fviz_pca_ind(res.pca,col.ind = group)

### 构建比较矩阵
design <- model.matrix(~group)
### 比较矩阵命名
colnames(design) <- levels(group)
design

### 2.线性模型拟合
fit <- lmFit(exprSet,design)
### 3.贝叶斯检验
fit2 <- eBayes(fit)
### 4.输出差异分析结果,其中coef的数目不能操过design的列数
### 此处的2代表的是design中第二列和第一列的比较
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)
### 这个数据很重要需要保存一下
save(allDiff,file = "output/allDiff.Rdata")

###################################################################################
### 定义差异基因：差异倍数2倍，矫正后的p值小于0.05

library(dplyr)
diffgene <- allDiff %>%
  filter(adj.P.Val < 0.05) %>%
  filter(abs(logFC) >1)

### 如果出现行名丢失的情况，需要先把行名变成列，处理完毕后再把列变成行名
### 这个工作是由tibble这个包里面的rownames_to_column()和column_to_rownames()完成的
library(tibble)
diffgene <- allDiff %>%
  rownames_to_column() %>%
  filter(adj.P.Val < 0.05) %>%
  filter(abs(logFC) >1) %>%
  column_to_rownames()

### 可选方案:使用subset直接获取,&是and的意思
diffgene <- subset(allDiff,abs(logFC) >1 & adj.P.Val < 0.05)
test <- allDiff[allDiff$adj.P.Val < 0.05 & abs(allDiff$logFC)>1,]
### 该数据也需要保存，此处一次性保存两个数据，如果是多个，一次写入变量名称即可。
save(diffgene,group,file = "output/diffgene.Rdata")
### 到此差异基因的分析就结束了
####################################################################################
####################################################################################

## 1.把现在数据调整成可以作图的格式
### 这个技能是data wrangling部分重点掌握的技能
### 复习一下流程：输入数据是表达量，经过三步
### 1.探针ID转换，2.行列转置，3，添加分组信息。最终获得的是数据框

### 行列转置
exprSet <- as.data.frame(t(exprSet))
### 添加分组信息
dd <- cbind(group=group,exprSet)
### 截取部分展示,这就是清洁数据
test = dd[,1:10]

## 2.作图展示
library(ggplot2)
ggplot(data = dd,aes(x=group,y=CD36,fill=group))+
  geom_boxplot()+
  geom_point()+
  theme_bw()

## 3.steal plot
my_comparisons <- list(
  c("treat", "con")
)
library(ggpubr)
ggboxplot(
  dd, x = "group", y = "CD36",
  color = "group", palette = c("#00AFBB", "#E7B800"),
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")

## 改写成函数
diffplot <- function(gene){
  my_comparisons <- list(
    c("treat", "con")
  )
  library(ggpubr)
  ggboxplot(
    dd, x = "group", y = gene,
    color = "group", palette = c("#00AFBB", "#E7B800"),
    add = "jitter"
  )+
    stat_compare_means(comparisons = my_comparisons, method = "t.test")
}

diffplot("CD36")
diffplot("MOXD1")

## 4.多个基因作图查看
## 先把基因提取出来
genelist <- rownames(diffgene)[1:6]
## 再提取表达量，使用名称选取行
data <- dd[,c("group",genelist)]
## 用pivot_longer调整数据，数据变长，增加的是行
library(tidyr)
data <- data %>%
  pivot_longer(cols=-1,
               names_to= "gene",
               values_to = "expression")
## 多基因作图
## 作图
ggplot(data = data,aes(x=gene,y=expression,fill=group))+
  geom_boxplot()+
  geom_jitter()+
  theme_bw()+
  stat_compare_means(aes(group=group), label = "p.signif", method = "t.test")

## 尝试更清晰的展示
ggplot(data = data,aes(x=group,y=expression,fill=group))+
  geom_boxplot()+
  geom_jitter()+
  theme_bw()+
  facet_grid(.~gene)+
  stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "t.test")

## 图片导出
library(export)
## 导成PPT可编辑的格式
graph2ppt(file="output/diffgenboxplot.pptx")
## 其他自己想要的格式
graph2pdf(file="output/diffgenboxplot.pdf")
graph2tif(file="output/diffgenboxplot.tif")
## 导成AI可以编辑的状态
graph2eps(file="output/diffgenboxplot.eps")