代碼都是直接給出來(lái)了
需要修改的地方我進(jìn)行了標(biāo)記
只要修改一下都能直接用
方法一:下載平臺(tái)數(shù)據(jù)以得到對(duì)應(yīng)信息
首先進(jìn)入官網(wǎng)https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi活鹰,在這里我以GSE26682為例,點(diǎn)擊紅圈中的鏈接
在這里插入圖片描述
在這里插入圖片描述
首先在分析前需要加載幾個(gè)必要的包,可以通過(guò)以下方法下載
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager",ask = F, update = F)
}
BiocManager::install("GEOquery")
BiocManager::install("limma")
BiocManager::install("affy")
只有一個(gè)平臺(tái)
gse <- getGEO("你的GSE序號(hào)", destdir = ".",
getGPL = T,
AnnotGPL = T) # 此處更改GSE序號(hào)
gse[[1]]
data <- read.csv("表達(dá)矩陣的地址") # 此處更改表達(dá)矩陣地址
# 表達(dá)矩陣也可以通過(guò)這個(gè)方法獲得
# exp <- exprs(gse[[1]])
GPL <- getGEO(filename = "剛剛下載的soft.gz文件的地址") # 此處增加soft.gz文件地址
gpl <- GPL@gpls[[1]]@dataTable@table
data <- aggregate(.~ID_REF, data, mean) # 對(duì)數(shù)據(jù)進(jìn)行去重,重復(fù)的取平均
colnames(gpl)
ids <- gpl[, c(a, b)] # 此處將a和b更改為前面列出的列名中的id和gene symbol的位置
colnames(ids) <- c("ID_REF", "symbol") # 此處將列名與表達(dá)矩陣中的列名進(jìn)行統(tǒng)一,不一定是ID_REF和symbol周伦,需要按照對(duì)應(yīng)的進(jìn)行修改
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, "http:///"))[1]),
stringsAsFactors = F)[, 1]
write.csv(re, file = "給文件取個(gè)名.csv") # 保存文件
有兩個(gè)平臺(tái)
library(GEOquery)
library(limma)
library(affy)
library(dplyr)
gse <- getGEO("你的GSE序號(hào)", destdir = ".",
getGPL = T,
AnnotGPL = T) # 此處更改GSE序號(hào)
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”列進(jìn)行提取,注意:可能叫ID_REF娃胆,不一定就叫“id”
gene_1$gene_symbol <- data.frame(sapply(gene_1$gene_symbol,
function(x)unlist(strsplit(x, "http://"))[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, "http://"))[1]),
stringsAsFactors = F)[, 1]
fdata_2 <- na.omit(fdata_2)
gene_2 <- fdata_2[, c("ID", "gene_symbol")] # 將fdata_1中的 “id”列和“gene_symbol”列進(jìn)行提取,注意:可能叫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, ] # 去重里烦,只保留第一個(gè)
data <- read.csv("表達(dá)矩陣的地址")
# 表達(dá)矩陣也可以通過(guò)這個(gè)方法獲得
# exp <- exprs(gse[[1]])
colnames(gene_total) <- c("ID_REF", "gene_symbol") # 此處將列名與表達(dá)矩陣中的列名進(jìn)行統(tǒng)一,不一定是ID_REF和symbol禁谦,需要按照對(duì)應(yīng)的進(jìn)行修改
re <- merge(data, gene_total, by.x = "ID_REF", by.y = "ID_REF")
write.csv(re_1, file = "給他取個(gè)名.csv")
注意事項(xiàng)
gene assignment 和 gene symbol的關(guān)系
有時(shí)候gene symbol會(huì)在gene assignment里面胁黑,需要我們自行進(jìn)行提取,這時(shí)候我么可以利用這個(gè)結(jié)構(gòu)(這個(gè)是需要自己理解一下代碼基礎(chǔ)用法的)州泊,只要改一下第二行最后方框里面的數(shù)字就可以提取“//”這個(gè)符號(hào)前面的還是后面的內(nèi)容
gene_2$gene_symbol <- data.frame(sapply(gene_2$gene_symbol,
function(x)unlist(strsplit(x, "http://"))[1]),
stringsAsFactors = F)[, 1]
gene symbol中存在空格的問(wèn)題
可以用一下的參數(shù)進(jìn)行更改列名(需要改一下名字等等)
names(fdata_2)[names(fdata_2) == "gene_assignment"] <- "gene_symbol"
使用bitr()函數(shù)
這是我在最開(kāi)始學(xué)id轉(zhuǎn)換時(shí)候用的最多的丧蘸,因?yàn)樗芊奖悖挥孟螺d什么亂七八糟的文件遥皂,對(duì)于我家這個(gè)很爛的網(wǎng)絡(luò)很有利力喷,但是有個(gè)致命的缺點(diǎn):不一定能對(duì)所有的id進(jìn)行轉(zhuǎn)換,也就是說(shuō)存在一個(gè)“正確率”的問(wèn)題演训,這樣的話弟孟,最后得出的結(jié)論也就存在一定的誤差,這樣就不太符合科研工作者嚴(yán)謹(jǐn)?shù)膽B(tài)度吧样悟,所以==不太建議使用==7髂肌!窟她!
只有在走投無(wú)路時(shí)候再使用陈症!
加載包
library(clusterProfiler)
library(org.Hs.eg.db)
看一下這個(gè)R包里面提供哪幾種數(shù)據(jù)類型的下載方式
keytypes(org.Hs.eg.db)
開(kāi)始轉(zhuǎn)換工作
ids = rownames(exp)
ids = as.data.frame(ids)
colnames(ids) = "ID"
result <- bitr(ids$ID, # 數(shù)據(jù)框
fromType = "PMID", # 你的ID的數(shù)據(jù)類型
toType = c("SYMBOL","ENSEMBL"), # 轉(zhuǎn)化的數(shù)據(jù)類型
OrgDb = org.Hs.eg.db) # org.Hs.eg.db——人類