如何判断芯片表达矩阵是否需要均一化
以GSE43292为例
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE43292
1.png
library(GEOquery)
setwd("/Users/apple/Desktop/GSE43292/")
##下载数据GSE43292
eSet <- getGEO("GSE43292",
destdir = '.',
getGPL = F)
# 从eSet中提取表达矩阵exprSet
exprSet <- exprs(eSet[[1]])
#如果已经下载好了matrix可直接读取
exprSet <- read.table(file = "./GSE43292_series_matrix.txt", header =TRUE, comment.char = "!", row.names=1)
第二种判断方法表达量在50以内，初步判断数据已做过标准化处理
第三种判断方法
#对得到的表达矩阵操作
ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprSet <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}
##判断是否需要log2转换脚本来自果子学生信简书
原始数据的均一化
#安装affy包
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("affy")
library(affy)
setwd("D:\\GEO数据挖掘与meta分析\\练习\\芯片数据的均一化及探针转换（代码）\\GSE43292_RAW")
rawdata <- ReadAffy() #读取当前文件下的CEL格式文件
#同时第一次还会从bioconductor上下载hugene10stv1用来注释cel文件
length(rawdata)#64
#两种归一化方法，推荐使用rma
eset.mas5 <- mas5(rawdata) #综合考虑pm和mm信号，假阳性高，exp数据没有取对数
eset.rma <- rma(rawdata)#基于robust multi-arrary average(RMA)算法衡量表达量，
#只使用pm信号，exp数据已经进行log2处理
#箱线图查看归一化的效果
boxplot(exprs(rawdata))
boxplot(exprs(eset.rma))
#保存表达矩阵
exprSet <- exprs(eset.rma)
write.table(exprSet, "expr_rma_matrix.txt", quote=F, sep="\t")
可看出与直接用getGEO下载的数据略有差异，后续分析还是用网getGEO下载的数据
探针注释
常用id
#提取soft数据信息
library(GEOquery)
setwd(setwd('D:\\GEO数据挖掘与meta分析\\练习\\芯片数据的均一化及探针转换（代码）'))
#下载GPL注释文件
GPL6244 <-getGEO('GPL6244',destdir =".")
GPL6244_anno <- Table(GPL6244)
#如果已经下载了soft文件
GPL6244_anno<-data.table::fread("./GSE43292_family.soft",skip ="ID")
colnames(GPL6244_anno)
[1] "ID" "GB_LIST" "SPOT_ID"
[4] "seqname" "RANGE_GB" "RANGE_STRAND"
[7] "RANGE_START" "RANGE_STOP" "total_probes"
[10] "gene_assignment" "mrna_assignment" "category"
使用R包提出探针和基因名称对应的关系
R有很多注释包，首先我们要知道什么平台用什么R包。 platformMap.txt文件（果子老师整理）
读入R语言
platformMap <- data.table::fread("platformMap.txt")
为了方便查询，果子老师写了以下的代码，每次只要修改平台的名称即可
index <- "GPL6244"
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")
[1] "hugene10sttranscriptcluster.db"(输出结果）
library(hugene10sttranscriptcluster.db)
probe2symbol_df <- toTable(get("hugene10sttranscriptclusterSYMBOL"))
得到探针与id的对应关系
正式进行探针ID之间的转换
准备表达矩阵
##下载数据GSE43292
eSet <- getGEO("GSE43292",
destdir = '.',
getGPL = F)
# 从eSet中提取表达矩阵exprSet
exprSet <- exprs(eSet[[1]])
#如果已经下载好了matrix可直接读取
exprSet <- read.table(file = "./GSE43292_series_matrix.txt", header =TRUE, comment.char = "!")
exprSet
我们要达到的效果是这个样的
把两个数据框按照探针名称合并， 同时在这里，我们还需要把重复的探针去掉，方法就是计算每一个探针的平均值，留下最大的那个。
先调整两个数据框的探针名字和类型，方便合并
names(exprSet)[1] <- names(probe2symbol_df)[1]
exprSet$probe_id <- as.character(exprSet$probe_id)
然后，
重点就来了，下面的代码恢弘大气，一次性完成这两个事情。
library(dplyr)
exprSet <- exprSet %>%
inner_join(probe2symbol_df,by="probe_id") %>% #合并探针的信息
select(-probe_id) %>% #去掉多余信息
select(symbol, everything()) %>% #重新排列，
mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>% #求出平均数(这边的.真的是画龙点睛)
arrange(desc(rowMean)) %>% #把表达量的平均值按从大到小排序
distinct(symbol,.keep_all = T) %>% # symbol留下第一个
select(-rowMean) %>% #反向选择去除rowMean这一列
tibble::column_to_rownames(colnames(.)[1]) # 把第一列变成行名并删除
write.csv(exprSet, "exprSet.csv", quote = F)
不得感叹一下：大神就是牛啊
参考：公众号：果子学生信
