【GEO数据库挖掘】二、ID转换

相关原始代码在这里:

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

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

我们下一篇再见!

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,080评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,422评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,630评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,554评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,662评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,856评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,014评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,752评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,212评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,541评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,687评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,347评论 4 331
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,973评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,777评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,006评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,406评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,576评论 2 349

推荐阅读更多精彩内容