rm(list = ls())
options(stringsAsFactors = F)
options(digits = 4) #设置全局的数字有效位数为4
load(file = 'step1-output.Rdata')

#导入的数据中，dat为表达矩阵，group_list为分组信息

#按照group_list分组画箱线图
boxplot(dat[1,]~group_list)
# 定义函数，用于和绘制箱式图
bp=function(g){
  # 高级绘图包，用于绘制自定义图
  # 比如添加P值之类的操作
  library(ggpubr)
  # 准备需要的数据
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  # 添加p值
  p + stat_compare_means()# 添加比较函数，默认使用wilcox
}
# 对4个样本的箱式图绘制
bp(dat[1,])

bp(dat[2,])

bp(dat[3,])

bp(dat[4,])

# 如果没有这个包，需要进行install，应该是biocmanger下载
library(limma)
# 数据预处理
# 分组矩阵准备
design <- model.matrix(~0+factor(group_list))
# 给分组矩阵添加分组名
colnames(design)=levels(factor(group_list))
head(design)

##   Control Vemurafenib
## 1       1           0
## 2       1           0
## 3       1           0
## 4       0           1
## 5       0           1
## 6       0           1

# 表达矩阵
exprSet=dat
# 给分组矩阵添加样本名
rownames(design)=colnames(exprSet)
design

##            Control Vemurafenib
## GSM1052615       1           0
## GSM1052616       1           0
## GSM1052617       1           0
## GSM1052618       0           1
## GSM1052619       0           1
## GSM1052620       0           1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`factor(group_list)`
## [1] "contr.treatment"

# 差异比较矩阵，也就是哪个和那个比的矩阵

contrast.matrix<-makeContrasts("Vemurafenib-Control",
                               levels = design)

# 进行差异比较的主要函数
# 这里需要输入三个文件，表达矩阵，分组矩阵，比较矩阵
deg = function(exprSet,design,contrast.matrix){
  ##线性模型拟合
  fit <- lmFit(exprSet,design)
  ##添加比较矩阵，进行差值计算
  fit2 <- contrasts.fit(fit, contrast.matrix)
  ## 贝叶斯检验
  fit2 <- eBayes(fit2)
  ## 生成检验结果报告
  tempOutput = topTable(fit2, coef=1, n=Inf)
  ## 去除缺失的数据
  nrDEG = na.omit(tempOutput)
  return(nrDEG)
}

# 调用函数
deg = deg(exprSet,design,contrast.matrix)

# 保存比较结果后续使用

save(deg,file = 'deg.Rdata')

# 这里的if函数用于折叠代码，方便阅读，没什么重要意义

nrDEG=deg# 重新赋值
attach(nrDEG)# 方便写代码
plot(logFC,-log10(P.Value))# 对p取对数

# 上面的图不太美观，进行美化
# 有点饿了，想吃烤串，一会儿学习
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value) #对p进行取负对数
# 散点图
ggscatter(df, x = "logFC", y = "v",size=0.5)

# 火山图倒置之后，继续对基因标注
# 对上调和下调的基因进行标注
df$g=ifelse(df$P.Value>0.01,'stable', #判断是否为稳定基因
            ifelse( df$logFC >2,'up',  #判断是否上调
                    ifelse( df$logFC < -2,'down','stable') ))#判断是否下调
# 查看有多少上调和下调基因
table(df$g)

##
##   down stable     up
##     51  18686     84

# 添加name列
df$name=rownames(df)
#散点图
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')

# 添加分组信息
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
          label = "name", repel = T,
          label.select = head(rownames(deg)), #挑选一些基因在图中显示出来
          palette = c("#00AFBB", "#E7B800", "#FC4E07") )

# 保存图片
# ggsave('volcano.png')

# MA图绘制
# MA图用于展示两组变化关系，(x+y)/2作为横轴，(y-x)作为纵轴
# 方便观察两组数值偏离程度（差异程度）
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)

df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
                ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )

##
## 0.001<p<0.01      p<0.001       p>0.01
##         2213         5291        11317

ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
          palette = c("green", "red", "black") )

# 保存结果
# ggsave('MA.png')

load(file = 'step1-output.Rdata')
x=deg$logFC # deg为差异比较的结果
names(x)=rownames(deg) # 将基因名添加
# 对差异排序，并取前后100的基因名，赋值到cg
cg=c(names(head(sort(x),100)),
     names(tail(sort(x),100)))
# 绘制热图
library(pheatmap)
# 对差异基因绘制热图
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)

# 对数据归一化后绘制热图
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
pheatmap(n,show_colnames =F,show_rownames = F)

# 添加分组信息
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,
         show_rownames = F,
        cluster_cols = F,
         annotation_col=ac)