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

siwei • • 杂文

相关原始代码在这里：

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

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

我们下一篇再见！

版权声明：
作者：siwei
链接：https://www.techfm.club/p/45897.html
来源：TechFM
文章版权归作者所有，未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>