表达谱数据中相同基因如何处理

在分析表达谱芯片的时候，我们经常会遇到多个探针对应同一个基因的情况。一般遇到这种情况，最常见的两种处理方法是

1）取平均

2）取表达值高的那个探针

那么今天我们就用R来实现这两种处理方式。至于，如何将探针转换成相应的基因名字，相对来说还是比较容易的。一般的芯片数据都会有一个相应的注释文件，从中可以找到探针对应的基因名字。对于一些Agilent的商用芯片和一些比较特殊的芯片平台，可能找不到探针的注释文件。前面我们也简单介绍过

☞探针注释文件中没有基因名字怎么办？

☞探针注释文件中没有基因名字怎么办？（二）

首先我们先来随便造一个基因名有重复的表达谱数据。

#设置随机过程的seed，保证结果可重复
set.seed(123)
#随机生成一个30行10列的矩阵
expr=matrix(runif(300,5,10),ncol=10)
#列名字为sample1-10
colnames(expr)=paste0("sample",1:10)
#行名从26个大写字母里面有放回的抽取30个字母，作为基因名
genes=sample(LETTERS,30,replace=T)
#合并得到基因名有重复的表达谱矩阵
expr=data.frame(genes,expr)
expr

接下来我们先用第一种方法

1）取平均

#利用aggregate函数，对相同的基因名按列取平均
expr_mean=aggregate(.~genes,mean,data=expr)
expr_mean

会得到如下结果，感兴趣的小伙伴可以随便挑几个check一下

2）对于重复的基因名字，取表达值最大的哪一行

其实aggregate也可以对相同的基因使用max函数取最大值，但是这样处理是有问题的。例如同一个基因出现了三次，那么会有三行数据。如果使用aggregate+max，对于每一个样本，他会从三个值中挑选最大的那个值最为这个样本的表达值，这样做是不科学的。我们先来看看效果

#利用aggregate函数，对相同的基因名按列取取最大值
expr_max=aggregate(.~genes,max,data=expr)
expr_max

原始数据

处理之后的数据

所以这个做法不可取。

对于相同的基因，我们应该挑选行平均值大的那一整行，而不应该打乱。

#计算行平均值，按降序排列
index=order(rowMeans(expr[,-1]),decreasing = T)
#调整表达谱的基因顺序
expr_ordered=expr[index,]
#对于有重复的基因，保留第一次出现的那个，即行平均值大的那个
keep=!duplicated(expr_ordered$genes)
#得到最后处理之后的表达谱矩阵
expr_max=expr_ordered[keep,]
expr_max

最后结果是这样的

表达谱数据中相同基因如何处理