相关原始代码在这里：
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
探针转化就到这里，下一篇我继续了解表达矩阵。
我们下一篇再见！
网友评论