这个基因芯片平台换了马甲就不认识了吗?
希望所有学员都可以站在生信技能树的舞台上发光发热!
追随生信技能树的脚步,学习生信已经有半年多了。看了哔哩哔哩上的视频,也跑了健明老师的代码。以为自己起码入门了,但是真正分析感兴趣的数据集时才发现还差得远。
今天尝试分析的是GSE52746,首先在探针ID注释上难住了。这个数据集是GPL17966平台,没查到对应的注释R包。
平台描述是:“This array is identical to GPL570 but the data were analyzed with a custom CDF environment”。我略微思考了一下,是不是可以用GPL570的注释包?从《用R获取芯片探针与基因的对应关系三部曲-bioconductor》,查到其注释包是 “hgu133plus2”。
按照健明老师的代码进行注释:
ids=toTable(hgu133plus2SYMBOL)
head(ids)
dim(ids)
colnames(ids)
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in% rownames(dat),]
dim(ids)
结果只匹配到175个基因!??
probe_id symbol
1 1053_at RFC2
2 117_at HSPA6
3 121_at PAX8
4 1255_g_at GUCA1A
5 1316_at THRA
6 1320_at PTPN21
> dim(ids)
[1] 41922 2
> colnames(ids)
[1] "probe_id" "symbol"
> ids=ids[ids$symbol != '',]
> ids=ids[ids$probe_id %in% rownames(dat),]
> dim(ids)
[1] 175 2
傻眼了,回头再看GPL17966网页介绍如下图。在健明老师的提醒下,注意到这里的ID并不是真正的探针ID,而是相应的Entrez GeneID;
而SPOT_ID就是基因完整的名了。
那么接下来就简单了
x <- read.csv(file="https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL17996&id=15209&db=GeoDb_blob105", comment.char='',
header = T, sep = "\t", skip = 22, na.strings=c("","NA"))
ids <- dplyr::select(x,-2)
colnames(ids) <- c("probe_id","symbol")
ids <- na.omit(ids)
ids <- filter(ids,!ids$symbol=="CONTROL")
ids <- ids[!duplicated(order(ids$symbol)),]
ids <- ids[ids$probe_id %in% rownames(dat),]
dat <- dat[ids$probe_id,]
dim(dat)
rownames(dat) <- ids$symbol
dat[1:4,1:4]dim(dat)
[1] 18908 39
根据GEO数据,目前基于 GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array)的基因芯片有61种,下图截取了部分。
看这些平台说明,都是注释后的。形式不尽相同(如下图),相应处理即可。我们随意截图其中两个看看:
GPL4454
GPL6671
我把3年前的收费视频课程:3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月,免费到B站:
这个课程超级棒,B站免费学习咯:https://m.bilibili.com/video/BV1dy4y1C7jz 配套代码在GitHub哈:https://github.com/jmzeng1314/GSE76275-TNBC TCGA数据库挖掘,代码在:https://github.com/jmzeng1314/TCGA_BRCA GTEx数据库挖掘,代码在:https://github.com/jmzeng1314/gtex_BRCA METABRIC数据库挖掘,代码在:https://github.com/jmzeng1314/METABRIC
然后马上就有了3千多学习量,而且有学员给出来了图文并茂版本万字笔记,让我非常感动!
扫描下面二维码马上就可以学习起来啦,笔记需要至少半个小时来阅读哦!