一. 實(shí)驗(yàn)?zāi)康?/h2>
學(xué)習(xí)如何下載甲基化數(shù)據(jù)及如何獲得每個(gè)探針的甲基化水平(beta值)。
二. 實(shí)驗(yàn)數(shù)據(jù)及步驟
數(shù)據(jù):來(lái)自GEO GSE27284_non-normalized_data.txt
57578*31探針表達(dá)譜矩陣 總共10個(gè)樣本,5個(gè)肺癌,5個(gè)正常,整個(gè)表達(dá)譜第一列是探針名,每個(gè)樣本有三列數(shù)據(jù),第一列是通道A的表達(dá)情況,第二列是通道B的表達(dá)情況,第三列是p值從GEO數(shù)據(jù)下載GPL8490_HumanMethylation27_270596_v.1.2.csv和GSE27284_non-normalized_data.txt
-
去掉頭文件,也就是提取GSE27284_non-normalized_data里從第五行到最后一行的數(shù)據(jù)(也就是探針表達(dá)譜數(shù)據(jù))保存為zGSE27284.csv,提取平臺(tái)文件里的探針基因列
image.png
image.png 利用R導(dǎo)入
利用R算出Beta值
探針對(duì)應(yīng)基因
(1) 刪除既有高甲基化位點(diǎn)和低甲基化位點(diǎn)的基因
(2) 一一對(duì)應(yīng)的不變
(3) 一個(gè)基因?qū)?yīng)多個(gè)探針,取最大值
(4) 寫(xiě)出表達(dá)譜
代碼
setwd("D:/實(shí)驗(yàn)/實(shí)驗(yàn)三")
#beta值
zGSE27284<-read.csv("zGSE27284.csv",header=T) #讀入數(shù)據(jù)
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] #表達(dá)譜列數(shù)
b<-dim(zGSE27284data)[1] #表達(dá)譜行數(shù)
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") #讀取平臺(tái)文件,第一列probeid 第二列g(shù)enesymbol,第三列g(shù)eneid
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")
#探針對(duì)應(yīng)基因
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,] #提取這個(gè)表達(dá)譜
aa<-table(geneid1) #行名為geneid
str(aa)
bb<-as.data.frame(aa)
unique1<-subset(bb,Freq=1)#取單獨(dú)的一一對(duì)應(yīng)的geneid 7125
noUnique<-subset(bb,Freq>=2)#把大于二的選出來(lái) 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 #一對(duì)多取最大值
Methy2<-Methy1%>%group_by(geneid1)%>%summarise_all(max) #一對(duì)多取最大值
Methy2<-as.matrix(Methy2)
write.table(Methy2,file="methy_exp.txt",sep="\t")#寫(xiě)出
結(jié)果文件

image.png

