根据方差截断筛选矩阵的R代码

Question

请参见下面的编辑 我想使用R过滤矩阵（基因表达数据），并保留只有高方差值的行（基因/探针）。例如，我想只保留在底部和顶部百分位数（例如低于20％和高于80％）具有值的行。我希望将我的研究限制在仅具有高方差的基因上，以进行下游分析。在R中是否有常用的基因过滤方法？

我的矩阵有18个样品（列）和47000个探针（行），其值已进行log2转换和归一化处理。我知道quantile()函数可以确定每个样本列中的20％和80％截止值。我无法找出如何在整个矩阵中找到这些值，然后对原始矩阵进行子集操作以删除所有“不变”的行。

以下是示例矩阵，平均值为5.97，因此最后三行应被删除，因为它们包含介于20％和80％截止值之间的值：

> m

                sample1 sample2 sample3 sample4 sample5 sample6
ILMN_1762337    7.86    5.05    4.89    5.74    6.78    6.41
ILMN_2055271    5.72    4.29    4.64    5.00    6.30    8.02
ILMN_1736007    3.82    6.48    6.06    7.13    8.20    4.06
ILMN_2383229    6.34    4.34    6.12    6.83    4.82    5.57
ILMN_1806310    6.15    6.37    5.54    5.22    4.59    6.28
ILMN_1653355    7.01    4.73    6.62    6.27    4.77    6.12
ILMN_1705025    6.09    6.68    6.80    6.85    8.35    4.15
ILMN_1814316    5.77    5.17    5.94    6.51    7.12    7.20
ILMN_1814317    5.97    5.97    5.97    5.97    5.97    5.97
ILMN_1814318    5.97    5.97    5.97    5.97    5.97    5.97
ILMN_1814319    5.97    5.97    5.97    5.97    5.97    5.97

我希望能得到任何建议，或者你认为我应该了解的函数。谢谢！

编辑

抱歉，我的问题可能不太清楚。(1) 我想知道整个矩阵的20％和80％截止值（而不仅仅是每个样本的）。 (2) 然后，如果任何一行包含上限或下限百分位数中的值，R将保留这些行。 如果一行包含接近平均值的值（对所有样本而言），那么这些行将被舍弃。

- Todd

---

感谢澄清。我也已经更新了我的答案，以反映您希望实现的内容。快速问题 - 您是否拥有矩阵或数据框（即ID列是矩阵的行名还是数据框的第一列？）。检查的快速方法是class(m)。 - Simon O'Hanlon

它应该是一个矩阵（仅表达式数据），ID列是我的矩阵的行名（我应该在示例中省略“ID”名称）。 - Todd

好的，太棒了！这正是我在示例中想到的。 - Simon O'Hanlon

最终我使用了下面建议的bioconductor包'genefilter'，使用以下代码：m.var <- varFilter(m, var.func=IQR, var.cutoff=0.6, filterByQuantile=TRUE)，并使用nrow(m)和row(m.var)来比较过滤后剩余探针的数量。 - Todd

Answer 1

好的，假设您有一个矩阵（因此我假设您的ID列实际上是行名），那么这非常容易做到。

#  First find the desired quantile breaks for the entire matrix
qt <- quantile( m , probs = c(0.2,0.8) )
# 20%  80% 
#5.17 6.62 
#  Next get a logical vector of the rows that have any values outside these breaks
rows <- apply( m , 1 , function(x) any( x < qt[1] | x > qt[2] ) )
#  Subset on this vector
m[ rows , ]
#            sample1 sample2 sample3 sample4 sample5 sample6
#ILMN_1762337    7.86    5.05    4.89    5.74    6.78    6.41
#ILMN_2055271    5.72    4.29    4.64    5.00    6.30    8.02
#ILMN_1736007    3.82    6.48    6.06    7.13    8.20    4.06
#ILMN_2383229    6.34    4.34    6.12    6.83    4.82    5.57
#ILMN_1806310    6.15    6.37    5.54    5.22    4.59    6.28
#ILMN_1653355    7.01    4.73    6.62    6.27    4.77    6.12
#ILMN_1705025    6.09    6.68    6.80    6.85    8.35    4.15
#ILMN_1814316    5.77    5.17    5.94    6.51    7.12    7.20

apply函数的any( x < qt[1] | x > qt[2] )部分用于在矩阵的边缘应用函数。如果该行中的任何值超出样本矩阵的20％和80％分位数，它将返回TRUE。根据定义，如果没有值超出这些限制，则返回FALSE，表示我们将在下一行中删除该行。请注意，HTML标签已保留。

- Simon O'Hanlon

---

感谢您的编辑！您的解决方案似乎非常直接。 - Todd

@Todd 太好了！如果这3个解决方案中的任何一个能够充分解决您的问题，您可能需要考虑将其中一个标记为被接受的答案，以便该问题不再显示为未回答（当然，如果您需要更多信息，请告诉我们）。如果您需要有关SO如何工作的更多信息，请参阅关于页面。干杯！ - Simon O'Hanlon

这个解决方案完美地回答了我在我的原始帖子中所问的问题！然而，我还在探索Bioconductor genefilter包中的其他可能非常有用的替代方法。非常感谢大家对我的R和stackoverflow知识有限的耐心支持。 - Todd

@SimonO'Hanlon 我也有同样的问题。在整个矩阵中选择前20%高变基因，我可以这样使用吗？ h.var <-varFilter(h, var.func=IQR, var.cutoff=0.8, filterByQuantile=TRUE) - beginner

Answer 2

Bioconductor的genefilter包提供了与微阵列分析相关的常见过滤器。基于行变异性的典型过滤器如下：

m = matrix(rnorm(47000 * 6), 47000)
varFilter(m)

该软件包的着陆页面引用了细节说明，展示了基本操作和提供过滤使用的诊断指导。

在微阵列分析中的一个基本原则是行内的值是可比较的，但不同行之间的值是不可比较的。这是因为与每一行相关联的探针具有独特的特征，会引入特定于行的偏差 - 第一行中的值可能合理地表明相对于第二行中相同样品的值，基因表达更多、更少或相等。这意味着@Todd希望基于行间比较（整个矩阵中最大和最小值）进行标准化的愿望是不被推荐的。相反，varFilter计算每行的变异性度量（行四分位距），并选择具有最大变异性的一部分（var.cutoff参数）。

快速查看varFilter的定义可以发现，一般来说，这不比对于某些行间变异性的度量var.func和一个（单一的）分位数var.cutoff更棘手。

vars <- apply(m, 1, var.func)
m[vars > quantile(vars, var.cutoff), ]

- Martin Morgan

---

+1！看起来不错。我做了以下操作，mm <- as.matrix(dat[,-1]);rownames(mm) <- dat[,1];rr <- varFilter(mm,var.cutoff=c(0.2,0.8))，但是我收到了警告，我有什么遗漏吗？ - agstudy

谢谢你的建议 - 我会立刻查看，因为可能有标准的方法来处理这个问题。 - Todd

感谢@Martin对基因组阵列过滤的“最佳实践”的解释。我计划在我的数据集中使用您的解决方案。 - Todd

@Martin 我也有同样的问题。为了在整个矩阵中选择前20%高变量基因，我可以这样使用吗？h.var <-varFilter(h, var.func=IQR, var.cutoff=0.8, filterByQuantile=TRUE) - beginner

Answer 3

我不是统计学家，所以不知道是否有一般的方法来解决这个问题。对我来说，如果您将数据重塑为长格式，问题会更简单。

library(reshape2)
dat.m <- melt(dat)
dat.m$value <- as.numeric(dat.m$value)
head(dat.m)
            ID variable value
1 ILMN_1762337  sample1  7.86
2 ILMN_2055271  sample1  5.72
3 ILMN_1736007  sample1  3.82
4 ILMN_2383229  sample1  6.34
5 ILMN_1806310  sample1  6.15
6 ILMN_1653355  sample1  7.01

然后对于每个变量，您需要执行以下操作：

1. 使用分位数计算限制。
2. 删除不符合条件的基因。

例如，您可以使用plyr中的ddply来实现此操作：

res <- ddply(dat.m,.(variable),function(x){
  ## compute limits for each sample
  z <- x$value
  qq <- quantile(z, probs = c(0.2,0.8))
  ## keep only genes with high or low variance
  dd <- x[z < qq[1] | z > qq[2],]
})
## return to the wide format
acast(res,ID~variable)

            sample1 sample2 sample3 sample4 sample5 sample6
ILMN_1653355    7.01      NA    6.62      NA    4.77      NA
ILMN_1705025      NA    6.68    6.80    6.85    8.35    4.15
ILMN_1736007    3.82    6.48      NA    7.13    8.20    4.06
ILMN_1762337    7.86      NA    4.89      NA      NA      NA
ILMN_1806310      NA      NA      NA    5.22    4.59      NA
ILMN_1814316      NA      NA      NA      NA      NA    7.20
ILMN_2055271    5.72    4.29    4.64    5.00      NA    8.02
ILMN_2383229      NA    4.34      NA      NA      NA      NA

编辑 根据问题澄清，如果您想要整个矩阵的20%和80%截止值而不仅仅是每个样本的值，则在ddply之外计算qq。

   qq <- quantile(dat.m$value, probs = c(0.2,0.8))

然后你可以像这样评论对应的行，例如：

res <- ddply(dat.m,.(variable),function(x){
  z <- x$value
  ## keep only genes with high or low variance
  dd <- x[z < qq[1] | z > qq[2],]
})

PS这里的dat是：

dat <- read.table(text='         ID    sample1 sample2 sample3 sample4 sample5 sample6
ILMN_1762337    7.86    5.05    4.89    5.74    6.78    6.41
ILMN_2055271    5.72    4.29    4.64    5.00    6.30    8.02
ILMN_1736007    3.82    6.48    6.06    7.13    8.20    4.06
ILMN_2383229    6.34    4.34    6.12    6.83    4.82    5.57
ILMN_1806310    6.15    6.37    5.54    5.22    4.59    6.28
ILMN_1653355    7.01    4.73    6.62    6.27    4.77    6.12
ILMN_1705025    6.09    6.68    6.80    6.85    8.35    4.15
ILMN_1814316    5.77    5.17    5.94    6.51    7.12    7.20
ILMN_1814317    5.97    5.97    5.97    5.97    5.97    5.97
ILMN_1814318    5.97    5.97    5.97    5.97    5.97    5.97
ILMN_1814319    5.97    5.97    5.97    5.97    5.97    5.97',header=TRUE)

- agstudy

---

+1 我认为我们可能都误解了 OP 的意思。我认为他们想要计算整个矩阵的 20% 和 80% 截止值，然后在这些值上进行子集操作。 - Simon O'Hanlon

是的，正是“我无法找到如何为整个矩阵找到这些值”，然后是“最后三行应该被删除，因为它们包含在20%和80%截止值之间的值”，这让我对我们最初的假设产生了怀疑。也许原帖作者可以澄清一下？ - Simon O'Hanlon

抱歉，我的问题描述不够清晰。(1) 我想知道整个矩阵的20%和80%截止值（而不仅仅是每个样本的）。然后，如果任何一行包含上下百分位数中的值，R将保留这些行。如果一行包含接近平均值的值（对于所有样本），则这些行将被排除。这有帮助吗？非常感谢您的建议 - 我在R方面还很新。 - Todd

@Todd 我编辑了我的答案。实际上，我建议你使用 M.M 的解决方案。 - agstudy

@agstudy 再澄清一下 - 数据框实际上是一个矩阵，而ID列实际上是行名 - OP将ID放在那里只是为了美观，而不是实际数据。 - Simon O'Hanlon