GPL18573 如何进行注释?
从官网下载的GPL18573 soft文件是一串gse,如何将id与基因名进行对应?
引用 皆我百晓生 小程序回复内容作答:
在R语言中对GPL18573进行注释,你可以使用Bioconductor中的相关包。首先,你需要安装biomaRt
和limma
包:
install.packages("BiocManager") # 如果未安装BiocManager,请先安装
BiocManager::install("biomaRt")
install.packages("limma")
然后,载入所需的库:
library(biomaRt)
library(limma)
接下来,你需要创建一个biomaRt
对象,与你希望进行注释的数据库连接。在这种情况下,我们将使用Ensembl数据库:
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
然后,你可以使用getBM
函数来获取GPL18573中的所有基因信息以及它们的注释。假设GPL18573中的列名为gene_ids
,你可以使用以下代码获取基因名和相关注释:
gene_ids <- c("ENSG00000157764", "ENSG00000230494")
annotations <- getBM(attributes = c("ensembl_gene_id", "external_gene_name", "description"),
filters = "ensembl_gene_id",
values = gene_ids,
mart = ensembl)
在上面的代码中,我们通过设置attributes
参数来指定我们需要的注释信息,然后使用filters
和values
参数指定我们要注释的基因列表。最后,我们将结果存储在annotations
对象中。
请注意,如果你的GPL18573文件不包含ensembl_gene_id
列,你需要根据R包中提供的注释查询其他列名,如gene_id
或gene_name
等。
希望这可以帮助到你!
去看他的说明文档,官网有对应的教程。
援引讯飞星火:
GPL18573是一个芯片数据文件,包含了基因表达数据。如果您想将id与基因名进行对应,可以按照以下步骤:
参考gpt:
结合自己分析给你如下建议:
GPL18573是一个基于Illumina NextSeq 500平台的人类基因芯片,可以用来检测人类基因的表达水平1。如果你从官网下载的GPL18573 soft文件是一串gse,那么可能是因为你下载的是GEO数据集的soft文件,而不是GPL平台的soft文件2。
GEO数据集的soft文件包含了多个样本的基因表达数据,以及对应的GPL平台信息。每个样本的基因表达数据都是以gse开头的一串数字表示的,这些数字是GEO数据库中分配给每个样本的唯一标识符3。如果你想要将这些gse与基因名进行对应,你需要先找到对应的GPL平台信息,然后使用相应的注释包或工具来进行注释。
一种简单的方法是使用R语言中的getGEO函数来下载并解析GEO数据集的soft文件,然后使用annotate函数来进行注释。下面是一个示例代码:
安装并加载相关的R包
install.packages(“Biobase”) install.packages(“GEOquery”) install.packages(“annotate”) library(Biobase) library(GEOquery) library(annotate)
下载并解析GEO数据集的soft文件
假设你下载的是GSE123456这个数据集
gse <- getGEO(“GSE123456”, GSEMatrix = TRUE)
查看该数据集使用的GPL平台
gpl <- Meta(gse)$platform_id
根据GPL平台选择相应的注释包
假设该数据集使用的是GPL18573这个平台
那么可以使用illuminaHumanv4.db这个注释包
如果没有安装这个注释包,需要先安装
source(“https://bioconductor.org/biocLite.R”) biocLite(“illuminaHumanv4.db”) library(illuminaHumanv4.db)
使用annotate函数来进行注释
将gse转换为表达矩阵
exprs <- exprs(gse[1])
获取探针ID
probe_id <- rownames(exprs)
获取基因名
gene_name <- getSYMBOL(probe_id, “illuminaHumanv4.db”)
将探针ID和基因名合并为一个数据框
anno <- data.frame(probe_id, gene_name)
GEO下载的数据在R中读取和注释
# 读取GEO下载的矩阵和注释
#### 加载包 ----
library(readxl)
library(tidyverse)
library(GEOquery)
# 读取 series_matrix文件 ----
GSE_data <- getGEO(filename = "data\\GSE14520-GPL571_series_matrix.txt.gz", getGPL = F)
pd_GSE_data <- pData(GSE_data) #观察临床信息中的data processing,Microarray suite,MAS 5.0,即标准化方法,如果已经经过了MAS,一般情况下不需要再标准化了,避免矫枉过正
pd_GSE_data$data_processing[1] #查看数据处理的方式
GSE_data_targets <- pd_GSE_data %>% # 提取有用的信息以便后续分析
dplyr::select(sample_id = geo_accession,
sample_name = title,
tissue_type = `Tissue:ch1`)
GSE_data_expr <- exprs(GSE_data) %>%
as.data.frame() # 提取表达矩阵
# 数据集的标准化结果查看------
boxplot(exprs(GSE_data)) ## 箱线图查看数据是否标准化
## PCA 查看不同样本间的区别
PCA_new <- function(expr, ntop = 500, group, show_name = F){
library(ggplot2)
library(ggrepel)
rv <- genefilter::rowVars(expr)
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(expr[select, ]))#最核心的代码
percentVar <- pca$sdev^2/sum(pca$sdev^2)
d <- data.frame(PC1 = pca$x[, 1],
PC2 = pca$x[, 2],
group = group,
name = colnames(expr))
attr(d, "percentVar") <- percentVar[1:2]
if (show_name) {
ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group")) +
geom_point(size = 2) +
xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) +
ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) +
geom_text_repel(aes(label = name),
size = 3,
segment.color = "black",
show.legend = FALSE )
} else {
ggplot(data = d, aes_string(x = "PC1", y = "PC2",color = "group")) +
geom_point(size = 2) +
xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) +
ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance"))
}
}
PCA_new(exprs(GSE_data),
nrow(GSE_data),
group = GSE_data_targets$tissue_type) # 组别以颜色显示
## getGEO 读取注释信息 ----
## AnnotGPL = T,用的是.gz的文件
GPL_data<- getGEO(filename ="data//GPL571.annot.gz", AnnotGPL = T)
GPL_data_11 <- Table(GPL_data)
## AnnotGPL = F,用的是.soft的文件
GPL_data_2 <- getGEO(filename = "data//GPL10687_family.soft.gz", AnnotGPL = F)
GPL_data_21 <- Table(GPL_data_2)
# 进行注释------
GSE_data_expr <- GSE_data_expr %>%
rownames_to_column() %>%
reshape::rename(c(rowname = "prode_id"))
GPL <- GPL_data_11 %>%
dplyr::select(ID,`Gene symbol`) %>%
reshape::rename(c(ID = "prode_id",`Gene symbol` = "GeneSymbol")) %>%
add_count(GeneSymbol,name = "n_gene") %>%
add_count(prode_id,name = "n_prode_id") %>%
dplyr::filter(n_gene < 1123) %>% # 选择能有对应基因的探针
dplyr::select(c(-3,-4))
GSE_expr <- GSE_data_expr %>%
inner_join(GPL, ., by = "prode_id") %>% #p2s和上一步得到的结果再取交际,p2s放在右边是以它为准
dplyr::select(-1) %>% #去除第一列probe_id
as.data.frame() %>% #因为aggregate需要数据框格式
aggregate(.~ GeneSymbol, data = ., mean) %>% #以symbol作为因子水平,把相似的数据放在一起取均值,最大值max,中位值median
column_to_rownames(var = "GeneSymbol")
## 可视化注释后的矩阵,因为注释过程丢失了部分探针
boxplot(GSE_expr)
PCA_new(GSE_expr,
nrow(GSE_expr),
group = GSE_data_targets$tissue_type) # 组别以颜色显示
# 以RData 格式保存信息 ----
GSE14520_GPL571 <- GSE_data
GSE14520_GPL571_expr <- GSE_expr
GPL571 <- GPL_data_11
GSE14520_GPL571_targets <- GSE_data_targets
picDir <- "dealed_data"
dir.create(picDir )
save(GSE14520_GPL571, # 矩阵
GSE14520_GPL571_targets, # metadata
GPL571, #注释信息
file = paste('.',picDir,"GEOqury读取的GSE14520-GPL571_data.RData",sep="/"))
save(GSE14520_GPL571_expr,
file = paste('.',picDir,"GSE14520-GPL571_expr.RData",sep="/"))