setwd("D:\\GSE140358衰弱")
library(GEOquery)
getGEO(GEO = "GSE140358",)
gset140358 <- getGEO('GSE140358',GSEMatrix = TRUE,destdir=".",
AnnotGPL = T,
getGPL = T)
save(gset140358,file="GSE140358.Rdata")
#取列表中的元素
a=gset140358[[1]]
dat_a=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
pd_a=pData(a) #用pData来提取临床信息
fdata_a<-fData(a)##提取平台信息
save(a,dat_a,pd_a,fdata_a,file = "exp_and_cli.Rdata")
load("file = exp_and_cli.Rdata")
load("D:\\GSE140358衰弱\\下载GEO数据\\exp_and_cli.Rdata")
##发现GPL16686平台不一样，
fdata_a1<-fdata_a[,c(1,6)]
fdata_a1$GB_ACC[2]
fdata_a2<-fdata_a1[fdata_a1$GB_ACC!="",]
#ID转换 用clusterprofiler转换"REFSEQ，ENTERZID,SYMBOL,ENSEMBL”,
library(clusterProfiler)
fdata_a3<-fdata_a2[!duplicated(fdata_a2$GB_ACC),]
sm<-bitr(fdata_a3$GB_ACC,fromType = "REFSEQ",toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
fdata_a4<-fdata_a3[fdata_a3$GB_ACC%in%sm$REFSEQ,]
KK<-merge(fdata_a4,sm,by.x="GB_ACC",by.y="REFSEQ")
da<-dat_a[rownames(dat_a)%in%KK$ID,]
da<-as.data.frame(da)
da$hj<-rownames(da)
da1<-merge(da,KK,by.x="hj",by.y="ID")
table(duplicated(da1$SYMBOL))
da2<-da1[!duplicated(da1$SYMBOL),]
rownames(da2)<-da2$SYMBOL
data<-da2[,-c(1,27,28)]
save(data,pd_a,file="IDchange.Rdata")
##差异基因
rm(list=ls())
load("D:\\GSE140358衰弱\\2ID转换\\IDchange.Rdata")
library(GDCRNATools)##一站式分析工具
##将样本分组
gg1<-pd_a[,c(2,36)]
gg1$`condition:ch1`<-ifelse(gg1$`condition:ch1`=="Robust","control","frailty")
dif<-gdcDEAnalysis(counts = data,group = gg1$`condition:ch1`,comparison = "frailty-control",method = "limma")
max(dif$logFC)
min(dif$logFC)
