甲基化数据预处理（R语言）

一. 实验目的

学习如何下载甲基化数据及如何获得每个探针的甲基化水平（beta值）。

二. 实验数据及步骤

1. 数据：来自GEO GSE27284_non-normalized_data.txt
   57578*31探针表达谱矩阵 总共10个样本，5个肺癌，5个正常，整个表达谱第一列是探针名，每个样本有三列数据，第一列是通道A的表达情况，第二列是通道B的表达情况，第三列是p值

2. 从GEO数据下载GPL8490_HumanMethylation27_270596_v.1.2.csv和GSE27284_non-normalized_data.txt

3. 去掉头文件，也就是提取GSE27284_non-normalized_data里从第五行到最后一行的数据(也就是探针表达谱数据)保存为zGSE27284.csv,提取平台文件里的探针基因列

   
   image.png
   image.png

4. 利用R导入

5. 利用R算出Beta值

6. 探针对应基因
   (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