热点新闻
【GEO数据库挖掘】二、ID转换
2023-08-06 11:35  浏览:4369  搜索引擎搜索“养老之家”
温馨提示:信息一旦丢失不一定找得到,请务必收藏信息以备急用!本站所有信息均是注册会员发布如遇到侵权请联系文章中的联系方式或客服删除!
联系我时,请说明是在养老之家看到的信息,谢谢。
展会发布 展会网站大全 报名观展合作 软文发布

相关原始代码在这里:

https://github.com/jmzeng1314/5years/blob/master/learn-R/tasks/3-r-20-codes.R

1 查看表达矩阵

https://www.jianshu.com/p/b69f1c84ecaa下载了数据并导入R,我们查看一下新建立的a data的具体内容:

> class(a) [1] "data.frame" > str(a) 'data.frame': 33297 obs. of 7 variables: $ X.ID_REF. : int 7892501 7892502 7892503 7892504 7892505 7892506 7892507 7892508 7892509 7892510 ... $ X.GSM1052615.: num 7.25 6.83 4.4 9.48 4.55 ... $ X.GSM1052616.: num 6.81 6.7 4.51 9.68 4.45 ... $ X.GSM1052617.: num 7.73 7.02 4.88 9.63 5.12 ... $ X.GSM1052618.: num 6.19 6.2 4.36 9.69 4.87 ... $ X.GSM1052619.: num 7.05 6.77 4.18 9.91 5.16 ... $ X.GSM1052620.: num 7.2 6.24 4.73 9.66 3.99 ... ##修改a列名并删除第一列 rownames(a)<- a[,1] a <- a[,-1]


image.png


首先回过头查看下载的数据平台:




image.png


为GPL6244,查询该平台使用的R包。

一般重要的芯片在R的bioconductor里面都是有包的,用一个R包可以批量获取有注释信息的芯片平台,不同GPL对应的R包如下:

http://www.bio-info-trainee.com/1399.html ##查询结果是hugene10sttranscriptcluster包,在R中安装加载 67 GPL6244 Homo sapiens hugene10sttranscriptcluster ##安装并运行包 BiocManager::install('hugene10sttranscriptcluster.db') library(hugene10sttranscriptcluster.db) ##原始代码 ids <- toTable(hugene10sttranscriptclusterSYMBOL) length(unique(ids$symbol)) tail(sort(table(ids$symbol))) table(sort(table(ids$symbol))) plot(table(sort(table(ids$symbol)))) ##运行代码 > length(unique(ids$symbol)) [1] 18857 > tail(sort(table(ids$symbol))) UBTFL1 IGKC CDK11B UBE2D3 LRRFIP1 LINC00965 6 7 8 8 10 26 > table(sort(table(ids$symbol))) 1 2 3 4 5 6 7 8 10 26 18097 597 131 15 7 5 1 2 1 1


image.png

##获取基因ID以及数据库列表ids,就可开始匹配 ##原始代码 exprSet <- a table(rownames(exprSet) %in% ids$probe_id) dim(exprSet) exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,] dim(exprSet) ids=ids[match(rownames(exprSet),ids$probe_id),] head(ids) exprSet[1:5,1:5] if(F){ tmp = by(exprSet,ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))] ) probes = as.character(tmp) dim(exprSet) exprSet=exprSet[rownames(exprSet) %in% probes ,] dim(exprSet) ##运行代码 > table(rownames(exprSet) %in% ids$probe_id) FALSE TRUE 13429 19868 ##有1万9千个探针是对应基因的,我们需要过滤表达矩阵 ##过滤之前是3万个探针,过滤后1万9千个探针 > dim(exprSet) [1] 33297 6 > exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,] > dim(exprSet) [1] 19868 6 ##ids顺序修改为和exprSet顺序相同 > ids=ids[match(rownames(exprSet),ids$probe_id),] > head(ids) probe_id symbol 1 7896746 MTND1P23 2 7896754 SEPTIN7P13 3 7896759 LINC01128 4 7896761 SAMD11 5 7896779 KLHL17 6 7896798 PLEKHN1 ##过滤有多个探针的基因,过滤了1千多个 > exprSet[1:5,1:5] X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. 7896746 10.16530 10.50290 10.13550 10.69020 10.34200 7896754 8.81485 8.99721 9.19503 9.49915 9.86029 7896759 8.75126 8.61650 8.81149 8.32067 8.41445 7896761 8.39069 8.52617 8.43338 9.17284 9.10216 7896779 8.20228 8.30886 8.18518 8.13322 8.06453 > tmp = by(exprSet,ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))] ) > probes = as.character(tmp) > dim(exprSet) [1] 19868 6 > exprSet=exprSet[rownames(exprSet) %in% probes ,] > dim(exprSet) [1] 18857 6

对于过滤这一步的详细讲解:

首先查看ids,选取有多个探针的基因如IGKC




image.png

##查看在ids中有多少个IGKC基因的重复探针,有7个 > table(ids[,2]=='IGKC') FALSE TRUE 19861 7 ##然后把exprSet矩阵恢复一下,运行最早的代码,查看IGKC的小表达矩阵 > exprSet[ids[,2]=='IGKC',] X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. 8043431 6.90341 6.74235 6.94234 6.42112 6.60051 8043459 7.81154 8.04823 7.81564 7.68827 7.79578 8043468 5.92133 6.03255 5.94167 5.96146 5.90188 8043480 8.68239 8.57171 8.70114 7.69504 8.29637 8043484 8.52665 8.45000 8.55041 7.68019 8.13244 8053735 8.67474 8.56958 8.69859 7.72255 8.28682 8161563 8.50079 8.44682 8.56368 7.67050 8.11585 X.GSM1052620. 8043431 6.46602 8043459 7.72170 8043468 6.18804 8043480 7.61538 8043484 7.67172 8053735 7.68026 8161563 7.66365 ##把这些基因的表达矩阵命名为x > x <- exprSet[ids[,2]=='IGKC',] > rownames(x) [1] "8043431" "8043459" "8043468" "8043480" "8043484" "8053735" "8161563" ##计算每一个探针的均值,取最大均值,第6个探针 > rownames(x) [1] "8043431" "8043459" "8043468" "8043480" "8043484" "8053735" "8161563" > rowMeans(x) 8043431 8043459 8043468 8043480 8043484 8053735 8161563 6.679292 7.813527 5.991155 8.260338 8.168568 8.272090 8.160215 > which.max(rowMeans(x)) 8053735 6 ##最后新建一个表达矩阵,从exprSet选取数据,按照ids$symbol分类,选取最大均值 tmp <- by(exprSet,ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))])

小结:

1、表达矩阵下载是Matrix

2、ids数据与表达矩阵对应,列名为probe_id和symbol

补充:

一些GPL的包是找不到的:

gpl <- getGEO('GPL6480', destdir=".") colnames(Table(gpl)) ## [1] 41108 17 head(Table(gpl)[,c(1,6,7)]) ##查看需要的列名 write.csv(Table(gpl)[,c(1,6,7)],"GPL6400.csv")##另存为csv

探针转化就到这里,下一篇我继续了解表达矩阵。

我们下一篇再见!

发布人:cfad****    IP:223.213.40.***     举报/删稿
展会推荐
让朕来说2句
评论
收藏
点赞
转发