- 获得平台探针注释信息
#source("https://bioconductor.org/biocLite.R")
#biocLite("GEOquery")
setwd("/Users/baiyunfan/desktop/GEO")
library(GEOquery)
GPL <- getGEO("GPL17843")
probe<-GPL@dataTable@table
这个probe矩阵里就有详细的探针信息啦～其中第15行是gene_symbol，接下来会用到
- 获取表达谱数据和表型矩阵
expr.df <- read.table(file = "GSE84957_series_matrix.txt", header =TRUE,
comment.char = "!", row.names=1)
Data <- getGEO(filename="GSE84957_series_matrix.txt")
pData <- pData(phenoData(Data))
这里的expr.df是读取我们预先下载好的表达谱矩阵，而getGEO是直接在线读取数据，二者选取一种方法即可～
- 将探针和symbol对应上
symbol<-sapply(1:dim(expr.df)[1],function(i){as.character(probe[which(probe[,1]==rownames(expr.df)[i]),15])})
expr.df<-cbind(symbol,expr.df)
- 由于多个探针对应着一个基因，所以把相同的基因的表达谱数据取平均数
更正：感谢@生信技能树@土豆学生信 的友情指导。经深入学习后发现，多个探针对应同一基因的情况可以取最大值，平均数和中位数，各有优缺点，用什么方法目前尚无定论。平均数是比较保守，本实验暂且用平均值来计算。
uni<-as.character(unique(expr.df[,1]))
uni<-uni[-2]
test<-matrix(data=NA,ncol=19,nrow=1)
test<-as.data.frame(test)
colnames(test)<-colnames(expr.df)
expr.df[,1]<-as.character(expr.df[,1])
expr.df<-expr.df[-which(expr.df[,1]==""),]
for(i in 1:length(uni)){if(length(which(expr.df[,1]==uni[i]))!=1){test<-rbind(test,c(uni[i],colMeans(expr.df[which(expr.df[,1]==uni[i]),-1])))}else{test<-rbind(test,expr.df[which(expr.df[,1]==uni[i]),])}}
test<-test[-1,]
就得到最终矩阵啦～
最后安利一波生信技能树的B站良心课程，免费的，免费的，免费的，质量好，质量好，质量好。
