GSEA文件准备及表达相关性分析(R语言)
程序员文章站
2022-05-19 11:57:43
...
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()
上一篇: R绘制相关性图
下一篇: JDBC与JDBC工具类