###差异分析
rm(list=ls())
load('exp_group.Rdata')
Sys.setenv(LANGUAGE='EN')
exp[1:4,1:4]
table(group_list)
先画箱线图看一下~
boxplot(exp[1,]~group_list)
t.test(exp[1,]~group_list)
library(ggpubr)
df <- data.frame(gene=exp[1,],group=group_list)
p <- ggboxplot(df,x='group',y='gene',
               color = 'group',palette = 'jco',
               add='jitter')
#  Add p-value
p + stat_compare_means()
##封装为函数
bp=function(g){         #定义一个函数bp，函数为{}里的内容
  library(ggpubr)
  df=data.frame(gene=g,group=group_list)
  p <- ggboxplot(df, x = "group", y = "gene",
                 color = "group", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
bp(exp[1,]) ## 调用上面定义好的函数，避免同样的绘图代码重复多次敲。
bp(exp[2,])
bp(exp[3,])
bp(exp[4,])

boxplot

limma差异分析

library(limma)
design=model.matrix(~factor( group_list ))#创建一个因子矩阵
#model.matrix创建一个矩阵模型（设计）
fit=lmFit(exp,design)
#lmFit 给定阵列数据为每个基因拟合线性模型
fit=eBayes(fit)
#eBayes 给出了一个微阵列线性模型拟合，通过经验贝叶斯调整标准误差到一个共同的值来计算修正后的t统计量、修正后的f统计量和微分表达式的对数概率。
## 上面是limma包用法的一种方式 
options(digits = 4) #设置全局的数字有效位数为4
#topTable(fit,coef=2,adjust='BH') 
topTable(fit,coef=2,adjust='BH') 
#coef可以是column number，也可以是column name，这样就可以指定你所感兴趣的两两比较的结果
#topTable从线性模型拟合中提取排名靠前的基因表。
bp(exp['241662_x_at',])
bp(exp['207639_at',])
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
save(deg,file='deg.Rdata')

5f184e00-0f2c-4f62-b797-4d2ec12aeea1.png

#step 4 基因id转换
rm(list = ls())
head(deg)
deg[1:5,1:5]#想知道上下调基因属于什么通路
deg$probe_id=rownames(deg) #deg重新取一列命名为probe_id,将deg的行名给新取的这一列作为每一行的内容
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
deg=merge(deg,ids,by='probe_id')#通过merge函数，由于deg和ids都有probe_id这一列，因此通过'probe_id'合并为新的deg
-deg
save(deg,file = 'deg.Rdata')

##火山图
rm(list=ls())
load('deg.Rdata')
head(deg)
deg$change=as.factor(
  ifelse(deg$P.Value>0.01, 'stable', 
         ifelse(deg$logFC>1.5,'UP',
                ifelse(deg$logFC< -1.5,'DOWN','stable'))
  )
)
table(deg$change)
deg$'-log10(pvalue)' <- -log10(deg$P.Value)
head(deg)
ggscatter(deg, x="logFC", y="-log10(pvalue)",
          color="change",size = 0.7,
          label = "symbol", repel = T, #自动调整标签位置
          label.select = c("AGR3","CYP4Z1","PROM1","SHC4"))

save(deg,file = "deg.Rdata")
rm(list = ls())

heatmap

热图

##热图
load("deg.Rdata")
head(deg)
tmp <- deg$logFC; names(tmp) <- deg$probe_id
head(tmp)
cg_up <- names(head(sort(tmp,decreasing = T),100))
cg_down <- names(head(sort(tmp),100))

load("exp_group.Rdata")
n=t(scale(t(exp[c(cg_up,cg_down),])))
n[n>2]=2; n[n< -2]= -2
n[1:4,1:4]
group_dat <- data.frame(group=group_list, row.names = colnames(exp))
library(pheatmap)
pheatmap(n,
         show_rownames = F,
         show_colnames = F,
         #cluster_cols = F,
         annotation_col = group_dat)

heatmap

参考：https://www.jianshu.com/p/207c8b8e6662