解决99%的GEO芯片平台探针基因注释问题！

已于 2024-12-23 17:43:51 修改
阅读量1.8k 收藏 26
点赞数 26
CC 4.0 BY-SA版权
分类专栏： R函数开发 GEO数据挖掘 文章标签： 数据库
于 2024-11-09 17:07:08 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/jaychouwong/article/details/143647706

获取更好阅读体验：解决99%的GEO芯片平台探针基因注释问题！

1. 前言

上一篇文章 听说你的使用getGEO经常有网络问题？一行代码解决下载GEO表达谱和临床数据 解决了GEO表达谱下载的问题。

然而，从众多客户的需求来看，大家对于后续的探针注释似乎觉得更难。

那么，这次我就教大家如何进行探针的基因注释，基本可以涵盖 99% 平台的情况。

B站已经同步讲解 本篇内容，复制BV号【BV1ifDoY4E1b】至B站可直接观看。

本教程所有代码、资料已打包，公众号【生信摆渡】后台回复【download_GEO】即可获取。

2. GPL文件的下载和读入

虽然有一些自动注释的R包，但大家还是要学习如何从GPL文件自己构建注释表的，毕竟R包更新是不及时的。

以GPL570为例。

首先看网页，包含了头部信息 Data table header descriptions 和一部分的表格 Data table。

头部信息包含了注释表格每一列的解释。在GPL文件中以#标识符开头，因此在读入时需要跳过。

接着展示注释表格的行数和文件大小。

再往下是下载按钮。不过有的平台不提供下载按钮，只能在线查看，而且打开后浏览器就卡住了。

这种情况可以将在线查看的HTML网页下载下来，然后进行读取处理也是可以的，下面会详细介绍。

读取：

  • 首先安装我编写的工具包
remotes::install_github("JiahaoWongg/bioinfobaidu.utils")
  • 下载下来后进行读取
read_GPL <- function(GPL_file){
    nskip = sum(subString(readr::read_lines(GPL_file, n_max = 100), 1) == "#")
    anno_tb = data.table::fread(GPL_file, skip = nskip)
    return(anno_tb)
}

GPL_file = "GPL570.txt"
raw_tb = read_GPL(GPL_file)
hd(raw_tb)
# dim: 54675 × 16 
#          ID GB_ACC SPOT_ID Species Scientific Name Annotation Date
# 1 1007_s_at U48705                    Homo sapiens     Oct 6, 2014
# 2   1053_at M87338                    Homo sapiens     Oct 6, 2014
# 3    117_at X51757                    Homo sapiens     Oct 6, 2014
# 4    121_at X69699                    Homo sapiens     Oct 6, 2014
# 5 1255_g_at L36861                    Homo sapiens     Oct 6, 2014

构建探针注释数据库

  • 提取注释信息。
  • 去除一探针对应多个基因的情况，只取第一个基因。
  • 统一将没有对应基因的探针注释为空 ""
anno_tb = data.frame(ID = raw_tb$ID, SYMBOL = raw_tb$`Gene Symbol`)
anno_tb$SYMBOL = subString(anno_tb$SYMBOL, 1, "[ /|]")
anno_tb$SYMBOL[anno_tb$SYMBOL == "NA"] = ""
anno_tb$SYMBOL[anno_tb$SYMBOL == "---"] = ""
hd(anno_tb)
# dim: 54675 × 2 
#          ID SYMBOL
# 1 1007_s_at   DDR1
# 2   1053_at   RFC2
# 3    117_at  HSPA6
# 4    121_at   PAX8
# 5 1255_g_at GUCA1A

这样注释文件就准备好了，为了方便使用，可以将每一个GPL注释数据库保存在一个列表里，并保存为RData，需要用时载入就行，这就构建了探针注释数据库了。

obj_file = "./anno_obj.RData"
if(!file.exists(obj_file)){
	anno_obj = list()
} else {
	anno_obj = get(load(obj_file))
}	
anno_obj[[GPL]] = anno_tb
save(anno_obj, file = obj_file)

进行注释

首先需要下载数据集的表达矩阵和样本信息。

使用之前的文章听说你的使用getGEO经常有网络问题？一行代码解决下载GEO表达谱和临床数据 的getGEO2函数可以一键下载。

并且这一次对getGEO2也进行更新了，第一点是自动检测和执行是否需要进行log转化，第二点是使用limma进行样本间标准化。

所以需要后台获取最新代码，当然，这次的代码也整合了新版本的getGEO2。

obj = getGEO2("GSE23586")
# INFO [2024-11-09 16:10:40] Querying dataset: GSE23586 ...
# 	Try time 1 ... succeed.
# - GPL: GPL570 
# - Found 6 samples, 38 metas.
# - Writting sample_anno to 00_GEO_data/GSE23586_sample_anno.xlsx. 
# INFO [2024-11-09 16:10:42] Execute log(expM + 1) ...
# INFO [2024-11-09 16:10:43] Normalize between arrays ...
# - Successed, file save to 00_GEO_data/GSE23586.RData.

GPL = obj$GPL
expM = obj$expM
  • 之后将原始表达矩阵的行名根据注释表替换为基因名就行了。
  • 可能存在多探针对应同一个基因的情况，可以选择的去重方法可以选择取平均mean或直接删掉重复的remove。建议取平均值，虽然计算要慢得多。如果很慢的话选择直接去除也是可以的。
anno_GEO <- function(expM, anno_tb, method = "mean"){
	# check match
	n_not_match = sum(! rownames(expM) %in% anno_tb$ID)
	if(n_not_match > 0){
		cat2("- Warnnig:", paste0(n_not_match, "/", nrow(expM)), "porbes fail to matched!\n")
	} else {
		cat2("- All porbes matched!\n")
	}

	# check SYMBOL
	gene_name = anno_tb$SYMBOL[match(rownames(expM),  anno_tb$ID)]
	n_no_symbol = sum(gene_name == "", na.rm = TRUE)
	if(n_no_symbol > 0){
		cat2("- Warnning:", paste0(n_no_symbol, "/", nrow(expM)), "probes fail to annotated!\n")
	} else {
		cat2("- All porbes annotated!\n")
	}

	expM = data.frame(gene_name, expM)
	expM2 = expM[expM$gene_name != "", ]

	report(paste("Removing duplicated genes by method:", method, "..."))
	if(method == "mean"){
		expM3 = aggregate(expM2[, -1], by = list(factor(expM2[, 1])), mean0)
		rownames(expM3) = expM3[, 1]
		expM3[, 1] = NULL
	} else if(method == "remove"){
		expM3 = expM2[!duplicated(expM2[, 1]), ]
		rownames(expM3) = expM3[, 1]
		expM3[, 1] = NULL
	}
	return(expM3)
}

anno_list = get(load(anno_file))
if(GPL %in% names(anno_list)){
	report("Use local annotation file ...")
	anno_tb = anno_list[[GPL]]
	expM_anno = anno_GEO(expM, anno_tb)
}
# INFO [2024-11-09 16:15:38] Use local annotation file ...
# - All porbes matched!
# - All porbes annotated!
# INFO [2024-11-09 16:15:38] Removing duplicated genes by method: mean ...
  • 并提示有多少探针找不到注释信息或者注释信息为空。
  • 至此，探针注释就完成了
hd(expM_anno)
# dim: 22881 × 6 
#          GSM578529 GSM578530 GSM578531 GSM578532 GSM578533
# A1BG      7.022058  7.809902  6.896205  7.037151  7.071827
# A1BG-AS1  7.168314  5.108887  6.941011  5.997464  4.591907
# A1CF      6.182022  6.001592  6.112716  6.823624  5.634210
# A2M       8.910845  9.665579  9.145603  9.548700  8.718978
# A2M-AS1   6.525060  7.602097  7.376438  7.080512  6.642209

R包自动注释

自动注释就是曾老师编写的AnnoProbe包。

我进行了一点的修改，就是前面所说的一探针对应多给予你的情况，这里貌似没有处理。

然后将注释数据库整理成统一的格式，方便统一调用。

anno_GEO_AnnoProbe <- function(expM, GPL, method = "mean"){
	# idmap包可以被 AnnoProbe 包替代
	# if(sum(c("idmap1", "idmap2", "idmap3") %in% data.frame(installed.packages())$Package) != 3){
	# 	cat2("Please install annotation packages 'idmap1' 'idmap2' 'idmap3' from github.\n")
	# 	stop()
    # }

	# loadp(idmap1, idmap2, idmap3)
	# anno_tb1 = try(suppressMessages(getIDs(GPL)), silent = TRUE)
	# anno_tb2 = try(suppressMessages(get_soft_IDs(GPL)), silent = TRUE)
	# anno_tb3 = try(suppressMessages(get_pipe_IDs(GPL)), silent = TRUE)
	# if(class(anno_tb1) == "data.frame"){
	# 	anno_tb = anno_tb1
	# } else if(class(anno_tb2) == "data.frame"){
	# 	anno_tb = anno_tb2
	# } else if(class(anno_tb3) == "data.frame"){
	# 	anno_tb = anno_tb3
	# } else {
	# 	cat2("Can't find annotation table, please DIY one.\n")
	# 	stop()		
	# }

	if(! "AnnoProbe" %in% data.frame(installed.packages())$Package){
		cat2("Please install annotation packages 'AnnoProbe' from github.\n")
		cat2("\tdevtools::install_github(\"jmzeng1314/AnnoProbe\").\n")
		stop()
    }
	
	report("Use AnnoProbe annotation file ...")
	loadp(AnnoProbe)
	idmap <- function (gpl = "GPL570", type = "bioc", mirror = "tercent"){
	    gpl = toupper(gpl)
	    gpl_anno = paste(gpl, c("bioc", "soft", "pipe"), sep = "_")
	    if (mirror == "tercent") {
	        up = "http://49.235.27.111"
	    }
	    if (!checkGPL(gpl)) {
	        stop("This platform is not in our list, please use our shinyAPP to custom annotate your probe sequences, or ask us to process and then update the R package!")
	    }
	    else {
	        tryCatch(utils::data("exists_anno_list", package = "AnnoProbe"))
	        gpl_anno = gpl_anno[gpl_anno %in% exists_anno_list]
	        if (T) {
	            tpf = paste0(paste(gpl, type, sep = "_"), ".rda")
	            down = paste0("/GEOmirror/GPL/", tpf)
	            download.file(paste0(up, down), tpf, mode = "wb", quiet = TRUE)
	            load(tpf)
	            return(get(paste(gpl, type, sep = "_")))
	        }
	        else {
	            stop("We have that platform, but just offer other type of annotaion.")
	        }
	    }
	}

	anno_tb = try(suppressMessages(idmap(GPL, type = "pipe")), silent = TRUE)
	if(class(anno_tb) == "try-error"){
		cat2("- Can't find annotation table in AnnoProbe.\n")
		stop()
	}
	colnames(anno_tb) = c("ID", "SYMBOL")
	anno_tb = anno_tb[!duplicated(anno_tb$ID), ]
	anno_tb$SYMBOL = subString(anno_tb$SYMBOL, 1, "[ /|]")
	anno_tb$SYMBOL[anno_tb$SYMBOL == "NA"] = ""
	anno_tb$SYMBOL[anno_tb$SYMBOL == "---"] = ""

	expM_anno = anno_GEO(expM, anno_tb)
	return(expM_anno)
}

我们使用数据集GSE106090测试：

GSE = "GSE106090"
obj = getGEO2(GSE, out_dir = GSE)
# INFO [2024-11-09 16:22:56] Querying dataset: GSE106090 ...
# 	Try time 1 ... succeed.
# - GPL: GPL21827 
# - Found 18 samples, 44 metas.
# - Writting sample_anno to GSE106090/GSE106090_sample_anno.xlsx. 
# INFO [2024-11-09 16:22:59] Normalize between arrays ...
# - Successed, file save to GSE106090/GSE106090.RData.

GPL = obj$GPL
expM = obj$expM
expM_anno = anno_GEO_AnnoProbe(expM, GPL)
# INFO [2024-11-09 16:23:00] Use AnnoProbe annotation file ...
# - Warnnig: 16541/61046 porbes fail to matched!
# - All porbes annotated!
# INFO [2024-11-09 16:23:01] Removing duplicated genes by method: mean ...

hd(expM_anno)
# dim: 31119 × 18 
#           GSM2829419 GSM2829420 GSM2829421 GSM2829422 GSM2829423
# A1BG        5.057959   5.030532   4.999752   5.199815   4.931398
# A2M-AS1     8.921258   9.073332   9.113745   9.010616   9.206460
# A2ML1       5.405560   5.161947   6.053355   5.855274   6.034109
# A2ML1-AS2   2.517003   2.352829   2.348781   2.351604   3.349039
# A2MP1       3.413715   3.385146   3.259971   3.766047   3.531293

也是非常滴好用。

在线获取注释

后来，我又发现，有些平台有提供在线的GPLxxx.annot.gz文件，内容和下载GPLxxx.txt文件内容是一样的（头信息外）。所以在自动注释不能用且还没有制作本地注释文件的时候，可以尝试一下。

感觉提供这个文件的平台并不是太多，但不如让代码去试试~

至于为什么GPL网页注释文件那个下载按钮只能点击而没有的链接，我也是不明白，有固定链接的话就不用手动点击了。

同样，编写一个函数进行注释：

anno_GEO_online <- function(expM, GPL, out_dir = "annot_dir", symbol_name = NULL){
    mkdir(out_dir)
    GPL_num = gsub("GPL", "", GPL)
    master_string = ifelse(nchar(GPL_num) < 4, "", subString(GPL_num, 1:(nchar(GPL_num) - 3)))
    master_dir = paste0("GPL", master_string, "nnn")
    url_prefix = "https://ftp.ncbi.nlm.nih.gov/geo/platforms"
    soft_url = paste0(url_prefix, "/", master_dir, "/", GPL, "/annot/", GPL, ".annot.gz")    
    soft_file = paste0(out_dir, "/", GPL, ".annot.gz")
    report("Download online annotation file ...")
    options(timeout = 100)
    dl_res = try(download.file(soft_url, destfile = soft_file, quiet = TRUE), silent = TRUE)
    if(class(dl_res) == "try-error"){
    	cat2("Can't find online annotation table. Please download annotation table manually:\n")
    	cat2(paste0("\thttps://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=", GPL, "\n"))
    	stop()
    }
    anno_tb = suppressWarnings(data.table::fread(soft_file, data.table = FALSE))
    
    if(is.null(symbol_name)){
        SYMBOL_col_names = c("Gene symbol", "GENE_SYMBOL", "GeneSymbol", "Symbol", "SYMBOL", "gene_assignment")
        symbol_name = SYMBOL_col_names[which(SYMBOL_col_names %in% colnames(anno_tb))]
        if(length(symbol_name) == 0){
	    	cat2("Can't find right 'symbol' column name. If exists exactly, please specific it to 'symbol_name'\n")
	    	stop()        	
        }
    }      
    anno_tb = anno_tb[, c("ID", symbol_name)]  
    colnames(anno_tb) = c("ID", "SYMBOL")
    anno_tb$SYMBOL = subString(anno_tb$SYMBOL, 1, "[ /|]")
    anno_tb$SYMBOL[anno_tb$SYMBOL == "NA"] = ""
    expM_anno = anno_GEO(expM, anno_tb)
    return(expM_anno)
}

我们还是使用数据集GSE23586进行演示：

GSE = "GSE23586"
obj = download_GEO(GSE, out_dir = "GSE23586", anno_file = anno_file)
# INFO [2024-11-09 16:31:44] Querying dataset: GSE23586 ...
# 	Try time 1 ... succeed.
# - GPL: GPL570 
# - Found 6 samples, 38 metas.
# - Writting sample_anno to GSE23586/GSE23586_sample_anno.xlsx. 
# INFO [2024-11-09 16:31:46] Execute log(expM + 1) ...
# INFO [2024-11-09 16:31:46] Normalize between arrays ...
# - Successed, file save to GSE23586/GSE23586.RData.

GPL = obj$GPL
expM = obj$expM
expM_anno = anno_GEO_online(expM, GPL)
# INFO [2024-11-09 16:31:48] Download online annotation file ...
# - All porbes matched!
# - Warnning: 9557/54675 probes fail to annotated!
# INFO [2024-11-09 16:31:54] Removing duplicated genes by method: mean ...

hd(expM_anno)
# dim: 21754 × 6 
#          GSM578529 GSM578530 GSM578531 GSM578532 GSM578533
# A1BG      7.022058  7.809902  6.896205  7.037151  7.071827
# A1BG-AS1  7.168314  5.108887  6.941011  5.997464  4.591907
# A1CF      6.182022  6.001592  6.112716  6.823624  5.634210
# A2M       8.910845  9.665579  9.145603  9.548700  8.718978
# A2M-AS1   6.525060  7.602097  7.376438  7.080512  6.642209

也肥肠滴好用~

整合代码

不用担心代码很多，我最终把这三种注释方式整合为一个函数download_GEO了，包含表达谱的下载和注释，同时输出样本注释文件为csv格式。

如果你能正常使用xlsx包的话，还会输出xlsx格式的文件，好处是表格的列宽是自适应的，不需要你再因为字挡住而拖来拖去的了（能用代码解决绝不用手。。。

最后，再示范一下，非常简单：

GSE = "GSE23586"
obj = download_GEO(GSE, anno_file = anno_file)
# Step1: Download GEO data ...
# INFO [2024-11-09 16:39:44] Querying dataset: GSE23586 ...
# 	Try time 1 ... succeed.
# - GPL: GPL570 
# - Found 6 samples, 38 metas.
# - Writting sample_anno to 00_GEO_data/GSE23586_sample_anno.xlsx. 
# INFO [2024-11-09 16:39:47] Execute log(expM + 1) ...
# INFO [2024-11-09 16:39:47] Normalize between arrays ...
# - Successed, file save to 00_GEO_data/GSE23586.RData.

# Step2: Annotate probe ...

# INFO [2024-11-09 16:39:49] Use local annotation file ...
# - All porbes matched!
# - All porbes annotated!
# INFO [2024-11-09 16:39:49] Removing duplicated genes by method: mean ...
# INFO [2024-11-09 16:39:52] Done.

names(obj)
# [1] "GPL"         "expM"        "sample_anno"
  • 结果就是一个列表，包含GPL索引、注释后的表达谱和样本注释信息。
  • 其中anno_file是我整理的注释文件，提供的压缩包里有。当然不提供这个文件也可以，就要碰碰运气看另外两种自动注释的方法能不能成功了。

这个数据集是可以的。

obj = download_GEO(GSE)

我相信，应该能够解决99%的注释问题，本地注释文件和代码我会一直更新，公众号关键词回复的链接也会同步更新。

学习更多生信技巧，请持续关注工种浩【生信摆渡】。

针对GEO，arrayexpress数据库寻找基因注释包-补充
weixin_55143227的博客
10-13 1733
寻找基因注释小技巧
GEO数据挖掘补充(四)——探针注释
MJ0417的博客
10-28 2436
来自生信技能树课程。
GEO生信数据挖掘（二）下载基因芯片平台文件及注释
热门推荐
zzh1464501547的博客
09-27 1万+
检索到目标数据集后，开始数据挖掘，本文以阿尔兹海默症数据集GSE1297为例上节我们，但是矩阵行名称是芯片探针需要转换为基因名。
GEO平台注释表读取器：从GEO读取并解析平台注释表（* .annot）。-matlab开发
05-30
该函数读取并解析NCBI Gene Expression Omnibus（GEO）的平台注释表（*.annot）。 平台注释表包含有关平台的信息。 输入： 文件名：*.annot 文件的文件名输出： annotation: 平台注解信息 platform_table: 平台表 例子： 1) 获取*.annot 文件例如来自ftp://ftp.ncbi.nlm.nih.gov/geo/platforms/GPL2nnn/GPL2883/annot/GPL2883.annot.gz 解压文件，你会得到：GPL2883.annot 2) 把函数和 annot 文件放在你的工作目录中，或者 addpath 3）运行功能[注释，平台表]=read_annot_file('GPL2883.annot');
GEO数据下载和整理
最新发布
2403_89168692的博客
05-10 1016
本文介绍了从GEO数据库下载基因表达和临床数据的流程，包括数据筛选、注释信息获取及整理、数据导出等步骤。首先，通过关键词检索和筛选条件（如物种、数据类型等）获取所需数据集，重点关注样本量、实验设计等质量指标。其次，通过TXT、Soft格式或R包获取注释信息，并使用R语言进行数据整理，包括表达矩阵的log转换、探针ID与基因symbol的匹配、重复基因处理等。最后，整理临床信息并导出表达矩阵和临床数据，为后续分析提供基础。
GEO数据库的基因注释+代码
XGJCNJ的博客
04-17 1836
先下载。
生信&数据挖掘——人工神经网络篇(2)GEO数据下载及注释
heart_6662的博客
01-20 6173
基因组学是一门研究基因组的科学，它涉及到基因组的结构，功能，表达，进化和疾病相关性等方面。基因组学的研究需要大量的基因数据，而这些数据可以从公共数据库中获得。其中最常用的数据库是GEO（Gene Expression Omnibus），它提供了大量的基因表达数据，可以用于基因组学研究。GEO数据库中的数据可以通过FTP或Web界面下载，下载的数据可以是原始数据，也可以是已经注释的数据。。
GEO的表达矩阵的探针ID转换成基因名称教程
qq_40073899的博客
07-12 2187
在基因表达数据分析中，了解特定的基因探针平台（GPL）的注释信息是非常重要的。`GEO`，或者称为`Gene Expression Omnibus`，是一个存储高通量基因表达数据（如微阵列和次世代测序数据）的公共数据库，由美国国立生物技术信息中心（NCBI）管理。每个数据集通常与一个或多个探针平台相关联，这些平台在GEO中被指定为GPL（Gene Expression Omnibus Platform）。
基于R语言的基因表达芯片注释流程.pdf
07-26
在生物信息学领域，基因表达芯片注释是一个重要的过程，用于确定芯片上每个探针序列所对应的基因信息。R语言作为一种开源且功能强大的统计分析语言，在处理基因表达数据方面表现出色。本文介绍了一种基于R语言的基因...
GEO芯片数据处理、多芯片数据合并、差异分析、核心基因提取、核心基因与免疫细胞相关性
2302_80650915的博客
01-20 5775
下载探针矩阵probe.txt和平台注释ann.txt查看分组信息，复制对照组样本名称s1.txt和实验组样本名称s2.txt。
GEO生信数据挖掘（三）芯片探针ID与基因名映射处理
zzh1464501547的博客
09-27 8389
检索到目标数据集后，开始数据挖掘，本文以阿尔兹海默症数据集GSE1297为例上节我们发现，存在一个芯片探针ID匹配到多个基因的情况，本节来介绍处理方案。
Affy芯片ATH1-121501探针注释的处理
R语言与生物信息学
12-04 5360
BioC有Affy芯片ATH1-121501（GPL198平台）的注释库文件（R包），芯片分析后期处理经常要使用，导出方法如下。 先安装BioC的ath1121501.db包： library("BiocInstaller") biocLite("ath1121501.db") 使用plyr揉数据包整理探针id和AGI： library(ath1121
【结合文献】——Affymatrix芯片数据预处理
weixin_40640700的博客
03-13 8602
【理论分析】 （待补充） 【实战】 数据来源：GSE98793 芯片平台：GPL570（[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array） 登录GEO数据库，在检索框中输入GSE98793，打开数据界面。 下载最下方的原始数据。 下载后大概有857MB大小。注意，使用浏览器下载时可能会出现数据下载不完整的...
GPL页面表格奇奇怪怪的更新
weixin_46587777的博客
05-22 1281
GPL页面表格奇奇怪怪的更新 0.因数据库和R包更新导致的问题 从GPL页面上下载注释表格，是获取探针注释的方法之一，通常是找不到注释包的第二选择。 更新之前的下载方式是(2022年之前我们的直播课堂上是这样讲的，2022年以及之后的课上已经更新好了)： library(GEOquery) a = getGEO("GPL570",destdir = ".") b = a@dataTable@table ids2 = b[,c("ID","Gene Symbol")] 更新之前：运行getGEO这句代码，函
This platform is not in our list, please use our shinyAPP to custom annotate your probe sequences,
shengxingaoshou的博客
12-17 253
遇到这种报错怎么办呢？
零基础入门转录组下游分析——数据处理（GEO数据库——芯片数据）
呆猪儿的博客
10-22 3390
零基础入门转录组下游分析—数据处理（GEO芯片数据），GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。
Supplementing your R package with a Shiny app
四纪启悟 · 数据科学
04-23 716
(This article was first published on  Dean Attali's R Blog, and kindly contributed to R-bloggers)      The R community is generally very fond of open-source-ness and the idea of releasing all code
Affymetrix公司芯片类型和对应平台
weixin_43250801的博客
12-12 1681
GPL-aff芯片
R语言 | GEO表达矩阵的数据清洗与预处理
tianyuu1的博客
12-03 5577
表达量矩阵的数据清洗应该在注释完成之后进行，并且下列操作最好按顺序进行如下图的表格所示，同一个探针ID对应的gene有多个，用///分隔着，而我们想获得一个探针ID只对应一个基因symbol的表格。表达矩阵注释过后，通常会有一些基因名是重复 Gene.symbol：是需要去重的所在例名data：是表达矩阵删除Symbol行
解释以下代码：library(GEOquery) library(limma) # load series and platform data from GEO gset <- getGEO("GSE62203", GSEMatrix =TRUE, AnnotGPL=TRUE) if (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx <- 1 gset <- gset[[idx]] # make proper column names to match toptable fvarLabels(gset) <- make.names(fvarLabels(gset)) # group membership for all samples gsms <- "00001111" sml <- strsplit(gsms, split="")[[1]] # log2 transformation ex <- exprs(gset) 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) if (LogC) { ex[which(ex <= 0)] <- NaN exprs(gset) <- log2(ex) } # assign samples to groups and set up design matrix counts<-exprs(gset) counts<-as.data.frame(counts) df<-read.csv("./gene_info.csv") #counts$gene_name=df[match(rownames(counts),df$ID),"symbol"] counts$ID=rownames(counts) rownames(counts)<-NULL counts=left_join(counts,df,"ID") str(counts);counts<-na.omit(counts);str(counts) final.counts<-counts[,c("symbol","GSM1522139","GSM1522140","GSM1522141","GSM1522142","GSM1522143","GSM1522144","GSM1522145","GSM1522146")]
03-08
<think>嗯，用户让我解释一段R代码，看起来是关于GEO数据集处理和分析的。首先，我需要仔细阅读代码，理解每一步的作用。用户可能是在学习生物信息学分析，特别是使用GEOquery和limma包进行基因表达分析。他们的需求可能不仅仅是代码解释，还想了解背后的原理，比如为什么进行log2转换，如何分组样本等。 代码开头加载了GEOquery和limma库，这两个包常用于处理GEO数据和差异表达分析。然后使用getGEO函数获取数据集GSE62203，这里可能需要解释GSEMatrix和AnnotGPL参数的作用。接下来检查数据集是否包含多个平台，选择GPL570，这可能是因为该平台是常用的Affymetrix芯片。 处理列名部分，用户可能对fvarLabels和make.names的作用不太清楚，需要说明这是为了后续分析匹配列名。分组信息gsms设置为"00001111"，这应该表示前四个样本是对照组，后四个是处理组。但代码里分割后得到8个字符，而后面样本有9个，这里可能存在不一致，需要指出潜在的问题。 log2转换的部分，通过分位数判断是否需要转换，用户可能不理解为什么这样做。需要解释基因表达数据通常需要对数转换以满足正态分布假设。然后处理counts数据，合并gene_info.csv中的基因信息，使用left_join，这里需要注意重复基因名或缺失值的问题，na.omit可能会删除部分数据，影响结果。 最后选择特定列生成final.counts，可能用于后续差异分析。用户可能想知道为什么要选这些样本，是否与实验设计有关。此外，代码中注释掉的一行给gene_name赋值，可能涉及探针ID转换，用户可能对这部分有疑问，需要说明其作用。 总结用户的需求，他们可能希望了解整个流程的步骤、每个函数的作用，以及潜在的问题和改进建议。需要指出可能的错误，比如样本数量不匹配，处理缺失值的方法是否合适，以及探针注释的正确性。确保解释清晰，帮助用户不仅理解代码，还能应用到其他项目中。</think>这个代码是用于处理基因表达数据并进行差异表达分析的典型流程。我为您逐步解释： 1. **数据加载与预处理**： ```R library(GEOquery) # 用于从GEO数据库下载数据 library(limma) # 用于差异表达分析 # 下载GSE62203数据集 gset <- getGEO("GSE62203", GSEMatrix=TRUE, AnnotGPL=TRUE) ``` - 通过`GEOquery`下载基因表达数据集 - `GSEMatrix=TRUE`表示返回表达矩阵格式 - `AnnotGPL=TRUE`会自动下载平台注释信息 2. **平台选择**： ```R if (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx <- 1 gset <- gset[[idx]] ``` - 选择GPL570平台（Affymetrix Human Genome U133 Plus 2.0 Array） 3. **数据标准化**： ```R 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) if (LogC) { ex[which(ex <= 0)] <- NaN exprs(gset) <- log2(ex) } ``` - 通过分位数判断数据是否需要log2转换 - 如果数据范围较大（第99分位数>100或全距>50），则进行log2转换 - 将<=0的值设为NaN（避免log(0)） 4. **样本分组**： ```R gsms <- "00001111" # 前4个样本为对照组，后4个为实验组 sml <- strsplit(gsms, split="")[[1]] ``` - 定义8个样本的分组（但后续实际使用9个样本，这里可能存在矛盾） 5. **基因注释整合**： ```R counts<-exprs(gset) df<-read.csv("./gene_info.csv") counts$ID=rownames(counts) counts=left_join(counts,df,"ID") ``` - 将表达矩阵与本地gene_info.csv中的基因符号信息合并 - `left_join`根据探针ID匹配基因符号 - `na.omit`删除缺失值 6. **最终数据整理**： ```R final.counts<-counts[,c("symbol","GSM1522139",...,"GSM1522146")] ``` - 提取需要的样本列（共9个样本）和基因符号列 需要注意的潜在问题： 1. 样本数量不一致：分组定义8个样本但最终选取9个样本 2. 被注释掉的`counts$gene_name=...`行可能需要取消注释 3. `na.omit`可能意外删除有效数据 4. 未显示差异表达分析的核心代码（可能limma分析在后续步骤） 典型后续分析步骤可能包括： 1. 创建设计矩阵：`design <- model.matrix(~0 + group)` 2. 线性模型拟合：`fit <- lmFit(final.counts, design)` 3. 计算差异表达：`eb <- eBayes(fit)` 4. 结果提取：`topTable(eb, coef=2, adjust="BH")` 建议改进点： 1. 使用`dplyr`包进行更清晰的数据操作 2. 添加样本名称检查确保分组正确 3. 增加探针到基因符号的重复处理（多个探针对应同一基因的情况） 4. 添加数据标准化后的质量评估图（如PCA、箱线图）
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值