GSEA文件准备及表达相关性分析（R语言）

GSEA文件准备

setwd("F:\\GEO\\GEO芯片数据/")

##下载好的载入
load('GSE35896_eSet.Rdata') 
a=gset[[1]] 
##取出第一个元素赋值给一个对象a
dat=exprs(a) 
#a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数，该函数得到表达矩阵
#现在 得到的dat就是一个表达矩阵，只不过基因的ID是探针名
dim(dat)
#看一下dat这个矩阵的维度
dat[1:5,1:5] 

##以下为GPL570的包
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL) 
head(ids) 

colnames(ids)=c('probe_id','symbol')  
length(unique(ids$symbol)) 
#[1] 18832个独特的基因探针，意味着本来19825个里面有一部分是重复的
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
#每个对象出现的个数
plot(table(sort(table(ids$symbol))))
#画图观察

ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in%  rownames(dat),]
##%in%用于判断是否匹配，然后取匹配的几行，去掉无法匹配的信息。

dat[1:5,1:5]   
dat=dat[ids$probe_id,] 
#取表达矩阵中可以与探针名匹配的那些，去掉无法匹配的表达数据，这时只剩下19825个探针及表达信息，其余已被剔除。

ids$median=apply(dat,1,median) 
#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序
##即先按symbol排序，相同的symbol再按照中位数从大到小排列，方便后续保留第一个值。
##将对应的行赋值为一个新的ids，这样order()就相当于sort()
ids=ids[!duplicated(ids$symbol),]
#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果.最后得到18832个基因。
dat=dat[ids$probe_id,] 
#新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol
#把ids的symbol这一列中的每一行给dat作为dat的行名
head(dat)
##至此我们就得到了该数据集的表达矩阵，最后将结果保存。
dat[1:3,1:3]
which(rownames(dat)=="MYC")

pd=pData(a) 
#通过查看说明书知道取对象a里的临床信息用pData
View(pd)
## 查看一下，挑选一些感兴趣的临床表型，这里我们欲得到其分组title信息。
library(stringr)
#运行一个字符分割包
group_list=str_split(pd$source_name_ch1,' ',simplify = T)[,4]
#抽取title一列，按照空格分割，取第四个元素即Control和Vemurafenib
table(group_list)
#看一下两个分组各有几个
group_list<-as.data.frame(group_list)
group_list$ID<-pd$geo_accession
rownames(group_list)<-group_list$ID
which(group_list$group_list=="primary")
rownames(group_list[1:566,])
Tumordat<-dat[,1:566]
colnames(Tumordat)
#Tumordat<-dat
###得到肿瘤表达矩阵,整理GSEA文件
all_df <- cbind(Name = rownames(Tumordat), DESCRIPTION = "na", Tumordat)
#all_df <- all_df[,-1]
dim(all_df)
all_df[1:3,1:3]
group <- c(rep("low", 31), rep("high", 31))
group <- paste(group, collapse = " ")
group <- c(paste(c(62, 2, 1), collapse = " "), "# low high", group)

write.table(file = "GEO35896.txt", all_df, sep = "\t", col.names = T, row.names = F, quote = F)
write.table(file = "group35896.cls", group, col.names = F, row.names = F, quote = F)

相关性分析

###相关性分析
MYC<-dat["MYC",]
GAPDH<-dat["GAPDH",]
data<-rbind(MYC,GAPDH)
data<-t(data)
data<-as.data.frame(data)

class(data)
###

library(ggplot2)
#求相关系数
r<-cor(data$MYC,data$GAPDH,method = "pearson")
#p值
p<-cor.test(data$MYC,data$GAPDH,alternative="two.side")
##画图
tiff(filename = "35896_cor.tiff")
ggplot(data =data,aes(x=MYC, y=GAPDH)) +
  geom_point(color="#d7191c")+
  geom_smooth(method="lm",color="#1a9641") +
  geom_text(aes(x=4.4, y=10,label=paste("R","=",signif(r,3),
                                        seq="")),
            color="black") +
  geom_text(aes(x=4.4, y=9.7,label="p = 0.105"),color="black")
dev.off()