10

请参阅下面的编辑 使用 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

我将不胜感激任何我应该研究的建议或功能。谢谢！

编辑

对不起，我在OP中不是很清楚。(1) 我想知道整个矩阵的 20% 和 80% 截止值（不仅仅是每个单独的样本）。(2) 然后，如果任何行包含上百分位数或下百分位数的值，R 将保留这些行。如果一行包含接近平均值的值（对于所有样本），则这些行将被丢弃。

4

3 回答 3

9

好的，假设你有一个矩阵（所以我假设你的 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

如果该行中的任何值超出样本矩阵的 20% 和 80% 分位数，则函数的any( x < qt[1] | x > qt[2] )部分apply（旨在将函数应用于矩阵的边距）返回。根据定义，如果没有值超出这些范围，则返回指示我们将在下一行删除该行。TRUEFALSE

于 2013-06-08T21:49:39.723 回答
5

Biocondcutorgenefilter 软件包提供了与微阵列分析相关的常用过滤器。基于逐行可变性的典型过滤器将是

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), ]
于 2013-06-08T22:15:53.323 回答
1

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

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. 删除不满足条件的基因。

例如，您可以使用 ddplyfrom执行此操作plyr：

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

在 OP 澄清后编辑，如果您想要整个矩阵的 20% 和 80% 截止值而不仅仅是每个单独的样本，您可以在 qq 之外计算ddply

   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 <- 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)
于 2013-06-08T22:01:48.023 回答