[TOC]
0.背景知识
trans_array这个函数的作用是讲探针表达矩阵转换为基因表达矩阵。
因为探针和基因不是一一对应的关系，会存在多个探针对应同一个基因的情况。例如AAMDC这个基因有三个探针：
rm(list = ls())
load("exp_ids.Rdata")
ids0 = ids[ids$symbol=="AAMDC",]
ids0
## probe_id symbol
## 24890 220268_at AAMDC
## 26046 221599_at AAMDC
## 26047 221600_s_at AAMDC
首先，多个探针对应同一个基因的问题是必须要处理的，因为最终的分析是以基因为单位，一个基因在一个样本里只能有一个表达量，在一次差异分析中只能有一个logFC，不能有多个啊。
我只设置了一种转换方式，就是按照symbol去重复，只保留每个基因的第一个探针。
其实全部的处理的方法有三种：
1.随机去重，保留任何一个探针都可以
2.保留行和/行平均值最大的探针
3.取多个探针的平均值
采取这三种方法都是可以的，没有标准答案哦！
1.随机去重
因为这些探针也没什么主次和顺序，所以直接用R语言里的去重复函数搞一下就可以。trans_array这个函数采用的就是直接去重，这个例子里就会保留第一个探针。
可以验证一下啊:
exp0 = exp[ids0$probe_id,]
exp0
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777
## 220268_at 3.801755 3.951920 4.004135 3.943231 4.213515
## 221599_at 8.429328 8.396561 8.369032 7.495108 7.612667
## 221600_s_at 7.922799 7.758738 7.765348 7.195104 7.262884
## GSM175778
## 220268_at 4.234231
## 221599_at 7.463824
## 221600_s_at 7.327311
library(tinyarray)
exp1 = trans_array(exp,ids)
exp1[1:4,1:4]
## GSM175773 GSM175774 GSM175775 GSM175776
## DDR1 8.126595 8.375204 8.216957 8.909610
## RFC2 6.028831 6.254986 6.517528 7.239964
## HSPA6 6.971390 6.586626 6.358798 5.501578
## PAX8 7.517333 7.501483 8.153091 8.537925
可以看到已经变成基因为行名的表达矩阵了。把AAMDC对应的表达量提取出来，和第一个探针的表达量放一起：
exp1["AAMDC",]
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 3.801755 3.951920 4.004135 3.943231 4.213515 4.234231
exp["220268_at",]
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 3.801755 3.951920 4.004135 3.943231 4.213515 4.234231
一模一样的。
2.保留行和/行平均值最大的探针
这个也是常见的去重策略之一。因为我不想这样搞(此处省略一些理由)，所以就没在函数里设置直接的参数。
但是我收到学生许愿了！主打一个有求必应，所以就来写了。
办法就是
load("exp_ids.Rdata")
ids = ids[ids$probe_id %in% rownames(exp),]
exp = exp[ids$probe_id,]
identical(ids$probe_id,rownames(exp)) #他俩的顺序已经相同了
## [1] TRUE
ids = ids[order(rowSums(exp),decreasing = T),] #ids2可以使用rowSums求和的顺序来排序，从大到小排，最大的在上面
library(tinyarray)
exp2 = trans_array(exp,ids)
%in% 和order都是R语言基础里面顶呱呱的函数。可以自行查阅帮助文档
让我们来检查一下：
rowSums(exp0)
## 220268_at 221599_at 221600_s_at
## 24.14879 47.76652 45.23218
（行平均值就是行和➗样本数量，样本数量是固定的，所以行和最大就是行平均值最大）
可见行和最大的探针是221599_at，按照这种方法，这个探针的表达量应该被作为基因表达量。check之：
exp2["AAMDC",]
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 8.429328 8.396561 8.369032 7.495108 7.612667 7.463824
exp0["221599_at",]
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 8.429328 8.396561 8.369032 7.495108 7.612667 7.463824
同样是一模一样哦
3.求平均值
这个方法不需要tinyarray
load("exp_ids.Rdata")
exp3 = exp[ids$probe_id,]
rownames(exp3) = ids$symbol
exp3[1:4,1:4]
## GSM175773 GSM175774 GSM175775 GSM175776
## DDR1 8.126595 8.375204 8.216957 8.909610
## RFC2 6.028831 6.254986 6.517528 7.239964
## HSPA6 6.971390 6.586626 6.358798 5.501578
## PAX8 7.517333 7.501483 8.153091 8.537925
exp3 = limma::avereps(exp3)
这种方法计算出来的是平均值，check之：
exp3["AAMDC",]
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 6.717961 6.702406 6.712838 6.211148 6.363022 6.341789
colMeans(exp0)
## GSM175773 GSM175774 GSM175775 GSM175776 GSM175777 GSM175778
## 6.717961 6.702406 6.712838 6.211148 6.363022 6.341789
网友评论