- 獲得平臺探針注釋信息
#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)
- 由于多個探針對應著一個基因，所以把相同的基因的表達譜數據取平均數
更正：感謝@生信技能樹@土豆學生信 的友情指導。經深入學習后發(fā)現，多個探針對應同一基因的情況可以取最大值，平均數和中位數，各有優(yōu)缺點，用什么方法目前尚無定論。平均數是比較保守，本實驗暫且用平均值來計算。
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站良心課程，免費的，免費的，免費的，質量好，質量好，質量好。
