合并重复探针
合并探针的原因
为了避免非特异性结合等干扰因素影响实验结果,芯片厂商往往采取多个探针检测同一基因表达的策略,从而导致注释探针后发现许多探针被注释为同一个基因。但在后续的分析中,程序往往不能接受表达矩阵中存在多个探针对应同一基因。因此,在进行后续分析之前,我们需要选取一个标准,对被注释为同一基因的探针进行合并。唯一要注意的是,要在过滤后再合并重复探针。
合并重复探针的方法
1 使用R语言合并重复探针
这个方法要求有一定的R语言功底。一般来说,我们会选择用探针的最大值或是中位数来合并重复探针。意思是只留下统计量最大的探针,其他都筛掉。这种方法能够一定程度上增加差异基因的数目,但也容易造成假阳性的结果。
使用ALL包的ALL数据作为例子。为了省时间和突出重点,在下文的例子中就先不对探针进行过滤。按照正常的步骤,应该先对探针进行过滤,然后才能进行合并重复探针的操作。
先注释探针。每个步骤的目的请见注释探针。
library(ALL)
data(ALL)
library(hgu95av2.db)
probe2symbol <- toTable(hgu95av2SYMBOL)
exprSet <- as.data.frame(exprs(ALL))
exprSet$probe_id <- rownames(exprSet)
exprSet_symbol <- merge(probe2symbol, exprSet, by = "probe_id")
exprSet_symbol$probe_id <- NULL
使用最大值合并探针。注意,请先运行命令length(which(is.na(exprSet_symbol)))
,若返回值大于0,说明存在NA值。需要先用impute包的命令 exprSet_symbol<- impute.knn(exprSet_symbol)$data
填充NA值,才能进行这一步。
exprSet_symbol <- aggregate(x = exprSet_symbol[,-1],by = list(exprSet_symbol$Symbol), FUN = max)
得到的exprSet_symbol即为合并探针之后的表达矩阵,可以用于后续分析。稍微解释一下上述代码。x = exprSet_symbol[,-1]
是进行整合的数据部分。by = list(exprSet_symbol$Symbol)
指对基因名进行整合。FUN = max
指使用最大值进行整合。
使用中位数合并探针得到的差异基因数量可能少于使用最大值合并探针。这方法更加保守,但更加不容易出错。代码如下,只是更改了FUN
参数。
exprSet_symbol <- aggregate(x = exprSet_symbol[,-1],by = list(exprSet_symbol$Symbol), FUN = median)
2 使用 findLargest 函数合并重复探针
findLargest
函数是整合与genefilter包中的函数之一,功能是根据最大的统计量合并重复探针。findLargest
函数的特点是无需事先对表达矩阵进行注释,函数会自动识别重复探针并返回统计量最大的探针的名称。findLargest
函数比起方法1还有个优点,那就是能够保留ExpressionSet
的结构。既可以进行后续的针对ExpressionSet
的分析,也可以提取表达矩阵。
不过,findLargest
函数需要引用芯片注释包,因此如果bioconductor中没有对应的芯片注释包的话,是不能使用findLargest
函数的。
使用上文的ALL数据作为例子。使用统计量IQR最大值作为合并探针的标准。统计量IQR根据每个探针在所有样本中的表达值计算。
library(ALL)
data(ALL)
library(hgu95av2.db)
library(genefilter)
IQRS <- apply(ALL, 1, IQR) #计算统计量IQR
probenames <- findLargest(rownames(ALL), testStat = IQRS, data = "hgu95av2.db")
#返回IQR值最大的探针名称
ALL <- ALL[probenames, ] #重复探针中只留下统计量IQR最大的探针
稍微解释一下findLargest
函数的用法。rownames(ALL)
表示探针的名称,testStat = IQRS
表示使用IQRS作为探针统计数据,data = "hgu95av2.db"
表示使用hgu95av2.db包识别重复探针。如果想用最大值进行合并的话,将IQRS <- apply(ALL, 1, IQR)
替换为MAXS <- apply(ALL, 1, max)
,probenames <- findLargest(rownames(ALL), testStat = IQRS, data = "hgu95av2.db")
替换为probenames <- findLargest(rownames(ALL), testStat = MAXS, data = "hgu95av2.db")
即可。
参考: