一. 实验目的
学习如何下载甲基化数据及如何获得每个探针的甲基化水平（beta值）。
二. 实验数据及步骤
数据：来自GEO GSE27284_non-normalized_data.txt
57578*31探针表达谱矩阵 总共10个样本，5个肺癌，5个正常，整个表达谱第一列是探针名，每个样本有三列数据，第一列是通道A的表达情况，第二列是通道B的表达情况，第三列是p值从GEO数据下载GPL8490_HumanMethylation27_270596_v.1.2.csv和GSE27284_non-normalized_data.txt
-
去掉头文件，也就是提取GSE27284_non-normalized_data里从第五行到最后一行的数据(也就是探针表达谱数据)保存为zGSE27284.csv,提取平台文件里的探针基因列
image.png
image.png 利用R导入
利用R算出Beta值
探针对应基因
(1) 删除既有高甲基化位点和低甲基化位点的基因
(2) 一一对应的不变
(3) 一个基因对应多个探针，取最大值
(4) 写出表达谱
代码
setwd("D:/实验/实验三")
#beta值
zGSE27284<-read.csv("zGSE27284.csv",header=T) #读入数据
str(zGSE27284)
zGSE27284data<-as.matrix(zGSE27284)
str(zGSE27284data)
M<-zGSE27284data[,seq(3,31,by=3)]#提取第通道B列
U<-zGSE27284data[,seq(2,31,by=3)]#提取第通道A列
M<-as.numeric(M)
U<-as.numeric(U)
a<-dim(zGSE27284data)[2] #表达谱列数
b<-dim(zGSE27284data)[1] #表达谱行数
for(i in a){
for(j in b){
ifelse(M>0,c<-M,c<-0)
}
}
str(c) #把0和M更大的赋给c
C<-matrix(c,27578,10)
beta<-C/(M+U+100)
tumorBeta<-beta[,1:5]
normalBeta<-beta[,6:10]
Nbeta<-rowMeans(normalBeta)
Tbeta<-rowMeans(tumorBeta)
methy<-Tbeta-Nbeta
probeid_geneid<-read.table("GPL8490-65-cg-gid.txt",header=T,sep="\t") #读取平台文件，第一列probeid 第二列genesymbol，第三列geneid
geneid<-probeid_geneid[,3]#提取geneid列
probeid<-probeid_geneid[,2]#probeid
probeid<-rownames(zGSE27284data)
methy1<-cbind(probeid,geneid,methy)
write.table(methy1,file="exp.txt",sep="\t")
#探针对应基因
Methy<-read.table("exp.txt",header=T,sep="\t",row.names=1)
#行名是probeid，第一列是geneid，第二列是甲基化水平，肿瘤-正常
#删除交集
library(dplyr)
high<-subset(Methy,methy>0) #高甲基化 14347
low<-subset(Methy,methy<0) #低甲基化 13231
jiaoji<-intersect(high$geneid,low$geneid) #取交集 5116
jiaojiweizhi<-which(jiaoji%in%geneid)
geneid1<-geneid[-jiaojiweizhi] #删除交集 22462
for(i in 1:nrow(Methy))
{
place<-which(Methy$geneid[i]!=geneid1)
}#返回位置
Methy1<-Methy[place,] #提取这个表达谱
aa<-table(geneid1) #行名为geneid
str(aa)
bb<-as.data.frame(aa)
unique1<-subset(bb,Freq=1)#取单独的一一对应的geneid 7125
noUnique<-subset(bb,Freq>=2)#把大于二的选出来 7983
noUnigene<-noUnique[,1] #geneid
write.table(noUnique,file="nogeneid.txt",sep="\t")
geneid2<-read.table("geneidun.txt",header=T,sep="\t",row.names=1)
geneid3<-geneid2[,1]
geneidfactor<-factor(geneid3)
Methy2<-apply(Methy1,2,function(x) tapply(x,geneidfactor,max)) #7293 #一对多取最大值
Methy2<-Methy1%>%group_by(geneid1)%>%summarise_all(max) #一对多取最大值
Methy2<-as.matrix(Methy2)
write.table(Methy2,file="methy_exp.txt",sep="\t")#写出
结果文件
image.png
