#三个包同时加载
library(AnnoProbe)
library(GEOmirror)
library(GEOquery) 

#下载获取GSE13507数据 
gset=AnnoProbe::geoChina('GSE13255')
gset


#检查表达式集 
eSet=gset[[1]] # 提取表达矩阵 
probes_expr <- exprs(eSet)
dim(probes_expr) 
#如果样本已经进行了log2处理，这里就不需要再进行归一化
probes_expr=log2(probes_expr+1) #表达矩阵行log2+1归一化处理

## 对表达芯片的探针进行基因注释
(gpl=eSet@annotation) 
checkGPL(gpl) 
printGPLInfo(gpl)
probe2gene=idmap(gpl) #获得探针注释
genes_expr <- filterEM(probes_expr,probe2gene )
class(genes_expr)
write.table(genes_expr,file="genes_expr.txt",quote=F,sep = "\t")

保存后，对数据进行整理，只保留想要分析的样本数据。

如果样本已经进行了log2处理，这里就不需要再进行归一化

报错

没有的话就自己整理
下载平台注释文件
整理成如下的文件，然后自己替换探针名

自行替换，然后再分组做差异分析

anno=read.table("anno.txt",sep="\t",header=T,check.names=F)
probes_expr<-as.data.frame(probes_expr)
ID<-as.data.frame(rownames(probes_expr))
probes_expr<-cbind(ID,probes_expr)
colnames(probes_expr)[1]<-'ID'

common<-merge(anno,probes_expr,by='ID')
common<-common[,-1]
write.table(common,file="genes_expr.txt",quote=F,sep = "\t")

差异分析

## 分组
group_list=factor(c(rep('Normal',91),rep('Cancer',137))) #自己定义分组和数量
table(group_list)

#limma包差异分析
library(limma)
library(data.table)
rt<-fread("genes_expr.txt",sep="\t")
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data<-as.data.frame(data)
nrow(data)
data=data[rowMeans(data)>0.5,]


design=model.matrix(~factor(group_list))
design
fit=lmFit(data,design)
fit=eBayes(fit)
DEG=topTable(fit,coef=2,n=Inf)
head(DEG)
write.table(DEG,file="DEG.txt",quote=F,sep = "\t")

绘图

##############绘制火山图
need_deg=data.frame(symbols=rownames(DEG), logFC=DEG$logFC, p=DEG$P.Value)

deg_volcano(need_deg,1) #第一种图
deg_volcano(need_deg,1,logFC_thred = 0) #设置logFC范围
deg_volcano(need_deg,2) # 第二种图


#################热图
deg_heatmap(DEG,data,group_list)
deg_heatmap(DEG,data,group_list,5) #显示前5对，数字可以自己定义


###########箱图
check_diff_genes('TLR7',data,group_list)
check_diff_genes('TLR7',data,group_list)+
  stat_compare_means(method = "t.test") #加个P值，手到擒来

临床数据下载

临床数据从sangerbox里面进行下载