1

我有一个矩阵，其中包含一列不同的基因，以及-log(P-values)每个基因的每个可能 SNP 的对应列。

所以矩阵有 3 列：Gene_lable、 SNP 和minus_logpval。我正在尝试编写一个代码来识别-log(P-value)每个基因最高的 SNP。这是头部（数据）：

  SNP           Gene_label           minus_logpval
1 rs3934834 HES4/ENSG00000188290       14.1031
2 rs3766193 HES4/ENSG00000188290        7.0203
3 rs3766192 HES4/ENSG00000188290       10.7420
4 rs3766191 HES4/ENSG00000188290       10.4323
5 rs9442371 HES4/ENSG00000188290       10.2941
6 rs9442372 HES4/ENSG00000188290        8.4235

这是代码的开始：

for(i in 1:254360) {
max_pval = 0
if(data$Gene_label[i]==data$Gene_label[i+1]) {
    x = array(NA, dim=c(0,2));
    x[i] = data$minus_logpval[i];
    x[i+1] = data$minus_logpval[i+1];
    temp = max(x);
    if (temp>max_pval) {
    max_pval=temp
    line = i
    }

但由于某种原因，R 一直给我错误：Error in is.ordered(x) : argument "x" is missing, with no default.我什至没有使用 is-ordered(x) 函数......我认为错误在于我初始化 x （应该是一个数组）的方式，但我没有不知道如何解决它。

4

2 回答 2

0

您可以使用 tapply 不循环地尝试它

tab <- expand.grid(gene=letters[1:2], SNP=LETTERS[1:3])
tab$minus_logpval <- abs(rnorm(6))*-1
tab <- tab[do.call("order", tab),]
tab$SNP <- as.character(tab$SNP)
with(tab, tapply(minus_logpval, gene, function(x) SNP[which.max(x)]))

高温高压

于 2013-06-18T09:46:36.547 回答
0

完美使用ddplyfrom plyr。将 up 拆分data.frame为子集 (by Gene_label) 并对每一部分进行操作（找到snp与maxof 相关的minus_logpval）：

##  Reproducible example data
set.seed(1234)
df <- data.frame( Gene_label = rep( letters[1:3] , 3 ) , snp = rep( letters[5:7] , each = 3 ) , minus_logpval = rnorm(9) )
df
#  Gene_label snp minus_logpval
#1          a   e    -1.2070657
#2          b   e     0.2774292
#3          c   e     1.0844412
#4          a   f    -2.3456977
#5          b   f     0.4291247
#6          c   f     0.5060559
#7          a   g    -0.5747400
#8          b   g    -0.5466319
#9          c   g    -0.5644520

##  And a single line using 'ddply'
require(plyr)
ddply( df , .(Gene_label) , summarise , SNP = snp[which.max(minus_logpval)] )
#  Gene_label SNP
#1          a   g
#2          b   f
#3          c   e
于 2013-06-18T09:46:48.197 回答