方法一:下载平台数据以得到对应信息
然后进入官网https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi,在这里我以GSE26682为例,点击红圈中的链接
然后下载soft.gz文件(一般是下面一个)
首先在分析前需要加载几个必要的包,可以通过以下方法下载
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager",ask = F, update = F)
}
BiocManager::install("GEOquery")
BiocManager::install("limma")
BiocManager::install("affy")
只有一个平台
gse <- getGEO("你的GSE序号", destdir = ".",
getGPL = T,
AnnotGPL = T) # 此处更改GSE序号
gse[[1]]
data <- read.csv("表达矩阵的地址") # 此处更改表达矩阵地址
# 表达矩阵也可以通过这个方法获得
# exp <- exprs(gse[[1]])
GPL <- getGEO(filename = "刚刚下载的soft.gz文件的地址") # 此处增加soft.gz文件地址
gpl <- GPL@gpls[[1]]@dataTable@table
data <- aggregate(.~ID_REF, data, mean) # 对数据进行去重,重复的取平均
colnames(gpl)
ids <- gpl[, c(a, b)] # 此处将a和b更改为前面列出的列名中的id和gene symbol的位置
colnames(ids) <- c("ID_REF", "symbol") # 此处将列名与表达矩阵中的列名进行统一,不一定是ID_REF和symbol,需要按照对应的进行修改
re <- merge(data, ids, by.x = "ID_REF", by.y = "ID_REF") # 此处同理,不一定是ID_REF
re <- re[!apply(is.na(re) | re$symbol == "", 1, all), ] # 删除空白值
re <- na.omit(re) # 去除na值
re$symbol <- data.frame(sapply(re$symbol,
function(x)unlist(strsplit(x, "///"))[1]),
stringsAsFactors = F)[, 1]
write.csv(re, file = "给文件取个名.csv") # 保存文件
有两个平台
library(GEOquery)
library(limma)
library(affy)
library(dplyr)
gse <- getGEO("你的GSE序号", destdir = ".",
getGPL = T,
AnnotGPL = T) # 此处更改GSE序号
fdata_1 <- fData(gse[[1]])
fdata_2 <- fData(gse[[2]])
sampleNames_1 <- sampleNames(gse[[1]])
sampleNames_2 <- sampleNames(gse[[2]])
gene_1 <- fdata_1[, c("ID", "gene_symbol")] # 将fdata_1中的 “id”列和“gene_symbol”列进行提取,注意:可能叫ID_REF,不一定就叫“id”
gene_1$gene_symbol <- data.frame(sapply(gene_1$gene_symbol,
function(x)unlist(strsplit(x, "//"))[1]),
stringsAsFactors = F)[, 1]
gene_1 <- na.omit(gene_1)
gene_2$gene_symbol <- data.frame(sapply(gene_2$gene_symbol,
function(x)unlist(strsplit(x, "//"))[1]),
stringsAsFactors = F)[, 1]
fdata_2 <- na.omit(fdata_2)
gene_2 <- fdata_2[, c("ID", "gene_symbol")] # 将fdata_1中的 “id”列和“gene_symbol”列进行提取,注意:可能叫ID_REF,不一定就叫“id”
gene_1 <- gene_1[!apply(is.na(gene_1) | gene_1$gene_symbol == "", 1, all), ]
gene_2 <- gene_2[!apply(is.na(gene_2) | gene_2$gene_symbol == "", 1, all), ]
gene_total <- rbind(gene_1, gene_2)
index_gene_total <- duplicated(gene_total$ID)
gene_total <- gene_total[!index_gene_total, ] # 去重,只保留第一个
data <- read.csv("表达矩阵的地址")
# 表达矩阵也可以通过这个方法获得
# exp <- exprs(gse[[1]])
colnames(gene_total) <- c("ID_REF", "gene_symbol") # 此处将列名与表达矩阵中的列名进行统一,不一定是ID_REF和symbol,需要按照对应的进行修改
re <- merge(data, gene_total, by.x = "ID_REF", by.y = "ID_REF")
write.csv(re_1, file = "给他取个名.csv")
注意事项
gene assignment 和 gene symbol的关系
有时候gene symbol会在gene assignment里面,需要我们自行进行提取,这时候我么可以利用这个结构(这个是需要自己理解一下代码基础用法的),只要改一下第二行最后方框里面的数字就可以提取“//”这个符号前面的还是后面的内容
gene_2$gene_symbol <- data.frame(sapply(gene_2$gene_symbol,
function(x)unlist(strsplit(x, "//"))[1]),
stringsAsFactors = F)[, 1]
gene symbol中存在空格的问题
可以用一下的参数进行更改列名(需要改一下名字等等)
names(fdata_2)[names(fdata_2) == "gene_assignment"] <- "gene_symbol"
使用bitr()函数
这是我在最开始学id转换时候用的最多的,因为它很方便,不用下载什么乱七八糟的文件,对于我家这个很烂的网络很有利,但是有个致命的缺点:不一定能对所有的id进行转换,也就是说存在一个“正确率”的问题,这样的话,最后得出的结论也就存在一定的误差,这样就不太符合科研工作者严谨的态度吧,所以不太建议使用!!!
只有在走投无路时候再使用!
加载包
library(clusterProfiler)
library(org.Hs.eg.db)
看一下这个R包里面提供哪几种数据类型的下载方式
keytypes(org.Hs.eg.db)
开始转换工作
ids = rownames(exp)
ids = as.data.frame(ids)
colnames(ids) = "ID"
result <- bitr(ids$ID, # 数据框
fromType = "PMID", # 你的ID的数据类型
toType = c("SYMBOL","ENSEMBL"), # 转化的数据类型
OrgDb = org.Hs.eg.db) # org.Hs.eg.db——人类