计算由多列分组的列中的每对可能的值答案

【问题标题】：Count every possible pair of values in a column grouped by multiple columns计算由多列分组的列中的每对可能的值
【发布时间】：2014-10-07 20:28:06
【问题描述】：

我有一个看起来像这样的数据框（这只是一个子集，实际上数据集有 2724098 行）

> head(dat)

chr   start  end    enhancer motif 
chr10 238000 238600 9_EnhA1  GATA6 
chr10 238000 238600 9_EnhA1  GATA4 
chr10 238000 238600 9_EnhA1    SRF 
chr10 238000 238600 9_EnhA1  MEF2A 
chr10 375200 375400 9_EnhA1  GATA6 
chr10 375200 375400 9_EnhA1  GATA4 
chr10 440400 441000 9_EnhA1  GATA6 
chr10 440400 441000 9_EnhA1  GATA4 
chr10 440400 441000 9_EnhA1    SRF 
chr10 440400 441000 9_EnhA1  MEF2A 
chr10 441600 442000 9_EnhA1    SRF 
chr10 441600 442000 9_EnhA1  MEF2A

我能够将我的数据集转换为这种格式，其中 chr、start、end 和 enhancer 组代表一个 ID：

> dat

 id motif 
 1  GATA6 
 1  GATA4 
 1    SRF 
 1  MEF2A 
 2  GATA6 
 2  GATA4
 3  GATA6 
 3  GATA4 
 3    SRF 
 3  MEF2A 
 4    SRF 
 4  MEF2A

我想找到每对可能的图案的数量，按 id 分组。 所以我想要一个像这样的输出表，

motif1 motif2 count
 GATA6  GATA4     3
 GATA6    SRF     2
 GATA6  MEF2A     2
 ... and so on for each pair of motif

在实际数据集中，有 1716 个独特的图案。有83509个唯一id。

对如何进行有什么建议吗？

【问题讨论】：

您的数据中有多少唯一的ids？
@Gregor 我将在几分钟内编辑我的问题以清晰明了。
我不知道这是否实用，但一种可能的方法是（1）将所有内容转换为整数，将您的分组变量压缩为一个 id，就像以前一样，并且然后（2）创建一个 1716x1716x83509 3-d 稀疏数组（使用slam 包），其中条目是存在/不存在的二进制，（3）您的结果是沿id 维度的colsums。
(2) 显然是最难的部分。或者，也许 SQL 可以处理它，这将是一个相当简单的自连接然后聚合查询。
@Gregor 没有 SQL 的背景，尽管感谢您的建议。

标签： r

【解决方案1】：

更新：这是一个使用data.table的快速和内存效率版本：

第 1 步：大致构建维度的样本数据：

require(data.table) ## 1.9.4+
set.seed(1L)        ## For reproducibility
N = 2724098L
motif = sample(paste("motif", 1:1716, sep="_"), N, TRUE)
id = sample(83509, N, TRUE)
DT = data.table(id, motif)

第 2 步：预处理：

DT = unique(DT) ## IMPORTANT: not to have duplicate motifs within same id
setorder(DT)    ## IMPORTANT: motifs are ordered within id as well
setkey(DT, id)  ## reset key to 'id'. Motifs ordered within id from previous step
DT[, runlen := .I]

第 3 步：解决方案：

ans = DT[DT, {
              tmp = runlen < i.runlen; 
              list(motif[tmp], i.motif[any(tmp)])
             }, 
      by=.EACHI][, .N, by="V1,V2"]

在最后一步 3 中这需要约 27 秒和约 1GB 的内存。

这个想法是执行一个self-join，但是利用data.table 的by=.EACHI 功能，它为每个i 评估j-expression，因此内存效率很高。而j-expression 确保我们只获得条目“motif_a,motif_b”而不是多余的“motif_b,motif_a”。这也节省了计算时间和内存。并且二进制搜索非常快，即使有 87K+ id。最后，我们通过主题组合聚合以获得每个主题组合中的行数 - 这就是您所需要的。

HTH

PS：请参阅旧（+ 较慢）版本的修订版。

【讨论】：

我想知道为什么 cc 对我来说是一个 12X1 向量而不是你。另外..运行此程序时出现以下错误：Error in bmerge(i <- shallow(i), x, leftcols, rightcols, io <- haskey(i), : x.'motif' is a factor column being joined to i.'V1' which is type 'NULL'. Factor columns must join to factor or character columns.。这是运行data.table 1.9.5`
@Mike.Gahan，我在加载时使用了stringsAsFactors=FALSE。当矩阵列是因素时，似乎as.data.table(.) 返回一个 1 列矩阵。这是您遇到的问题吗？如果您可以看看为什么/发生了什么并提出问题，您会很乐意解决吗？
现在我添加了stringsAsFactors=FALSE，工作完美。何时使用 as.data.table 与 setDT 或仅使用 data.table() 是否有技巧？另外...与foverlaps 一起工作很棒。很棒的新功能。
@Mike.Gahan，谢谢:-)。已用更快的版本更新了答案。尺度非常好。 setDT 仅适用于列表和 data.frames，而不适用于矩阵（由于它的内部结构）。所以，我们不能在这里避免as.data.table()，因为combn 返回矩阵。在其他地方，我更喜欢setDT...
@Arun 这是一个非常聪明的答案。谢谢！我需要提高我的数据表技能。

【解决方案2】：

这里无耻地借用this question的稀疏矩阵技术。

# Create an id
dat$id <- as.factor(paste(dat$chr, dat$start, dat$end, dat$enhancer))

# Create the sparse matrix.
library(Matrix)
s <- sparseMatrix(
      as.numeric(dat$id), 
      as.numeric(dat$motif),
      dimnames = list(levels(dat$id),levels(dat$motif)),
  x = TRUE)

co.oc <- t(s) %*% s # Find co-occurrences.
tab <- summary(co.oc) # Create triplet representation.
tab <- tab[tab$i < tab$j,] # Extract upper triangle of matrix

data.frame(motif1 = levels(dat$motif)[tab$i],
           motif2 = levels(dat$motif)[tab$j],
           number = tab$x)

#    motif1 motif2 number
# 1  GATA4  GATA6      3
# 2  GATA4  MEF2A      2
# 3  GATA6  MEF2A      2
# 4  GATA4    SRF      2
# 5  GATA6    SRF      2
# 6  MEF2A    SRF      3

【讨论】：

@Arun 我怀疑有一种方法可以通过只计算矩阵的“上（或下）三角形”来加快速度。

【解决方案3】：

我认为data.table 包可能是这里最有效的。我们可以计算每个组中的对，然后聚合。与首先计算所有对的总数相比，使用您的大小的数据是一种更有效的方法。

#Bring in data.table and convert data to data.table
require(data.table)
setDT(dat)

#Summarize by two-way pairs
summ <- dat[ , list(motifs=list(combn(unique(as.character(motif)),
   min(2,length(unique(as.character(motif)))), by=list(chr,start,end,enhancer)]

#Transpose and gather data into one table
motifs.table <- rbindlist(lapply(summ$motifs,function(x) data.table(t(x))))

#Summarize table with counts
motifs.table[ , .N, by=list(V1,V2)]

#       V1    V2 N
# 1: GATA6 GATA4 3
# 2: GATA6   SRF 2
# 3: GATA6 MEF2A 2
# 4: GATA4   SRF 2
# 5: GATA4 MEF2A 2
# 6:   SRF MEF2A 3

【讨论】：

谢谢。当我尝试第二个命令 summ
一定有一些组只有 1 个主题。我想这个答案并不能完全解释这一点。我会修改，但@Aruns 的答案要好得多。
好的。不过感谢您的帮助。 +1 试图提供帮助。

【解决方案4】：

如果您可以将数据放入名为dat 的 SQL 表中，则该查询应该可以工作：

select d1.motif m1, d2.motif m2, count(*) count
from dat d1
join dat d2
on d1.chr = d2.chr
  and d1.start = d2.start
  and d1.end = d2.end
  and d1.enhancer = d2.enhancer
  and d1.motif <> d2.motif
group by d1.motif, d2.motif

考虑到您的数据大小，我怀疑 R sqldf 包可以处理它，但是通过免费的 MySQL 安装，您可以使用 RODBC 或 RJDBC 进行 R 和 SQL 对话。

【讨论】：

您几乎可以肯定地用 R 函数模仿这种行为。
是的，你可以。我只是对data.table 感到不舒服，而且我对使用其他任何东西的 RAM 持怀疑态度。我认为这是@Arun 的答案中实现的逻辑。

【解决方案5】：

您可能会从数据语义的正式建模中受益。如果您有基因组范围，请使用 Bioconductor 的 GenomicRanges 包。

library(GenomicRanges)
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)

这是一个 GRanges 对象，它正式理解基因组位置的概念，因此这些操作可以正常工作：

hits <- findMatches(gr, gr)
tab <- table(motif1=gr$motif[queryHits(hits)],
             motif2=gr$motif[subjectHits(hits)])
subset(as.data.frame(tab, responseName="count"), motif1 != motif2)

【讨论】：

【解决方案6】：

这个呢？

res1<- split(dat$motif,dat$id)
res2<- lapply(res1,function(x) combn(x,2))
res3<- apply(do.call(cbind,res2),2,function(x) paste(x[1],x[2],sep="_"))

table(res3)

【讨论】：

【解决方案7】：

...如果这不是你想要的，我就放弃了。显然它没有针对大型数据集进行优化。这只是一个利用 R 的自然优势的通用算法。有几个可能的改进，例如与dplyr 和data.table。后者将大大加快这里的[ 和%in% 操作。

motif_pairs <- combn(unique(dat$motif), 2)
colnames(motif_pairs) <- apply(motif_pairs, 2, paste, collapse = " ")
motif_pair_counts <- apply(motif_pairs, 2, function(motif_pair) {
  sum(daply(dat[dat$motif %in% motif_pair, ], .(id), function(dat_subset){
    all(motif_pair %in% dat_subset$motif)
  }))
})
motif_pair_counts <- as.data.frame(unname(cbind(t(motif_pairs), motif_pair_counts)))
names(motif_pair_counts) <- c("motif1", "motif2", "count")
motif_pair_counts

#   motif1 motif2 count
# 1  GATA6  GATA4     3
# 2  GATA6    SRF     2
# 3  GATA6  MEF2A     2
# 4  GATA4    SRF     2
# 5  GATA4  MEF2A     2
# 6    SRF  MEF2A     3

另一个旧版本。请确保您的问题是明确的！

这正是plyr 旨在实现的目标。试试dlply(dat, .(id), function(x) table(x$motif) )。

但是请不要在不阅读文档的情况下尝试复制和粘贴此解决方案。 plyr 是一个非常强大的包，对你理解它会很有帮助。

回答错误问题的旧帖：

您是否在寻找不相交或重叠的对？

这是使用包zoo 中的函数rollapply 的一种解决方案：

library(zoo)

motif_pairs <- rollapply(dat$motif, 2, c)              # get a matrix of pairs
motif_pairs <- apply(motif_pairs, 1, function(row) {   # for every row...
  paste0(sort(row), collapse = " ")                    #   sort the row, and concatenate it to a single string
                                                       #   (sorting ensures that pairs are not double-counted)
})
table(motif_pairs)                                     # since each pair is now represented by a unique string, just tabulate the string appearances

## if you want disjoint pairs, do `rollapply(dat$motif, 2, c, by = 2)` instead

如果这不是您所需要的，请查看rollapply 的文档。要按其他变量分组，您可以将其与以下之一结合使用：

基本 R 函数 aggregate 或 by（不推荐），或
*ply 函数来自 plyr（更好）

【讨论】：

我已将变量分组到一个 id 中。
不相交和重叠对是什么意思？我只想计算一对在每个组中出现的次数。如何将输出中的数字转换为主题名称？
@KomalRathi 如果这不正确，您需要更清楚地了解“对”的含义。
我不是有意冒犯您，也绝不是在暗示您不正确。我的意思是问你什么是不相交和重叠的对，然后我只是试图解释我需要什么。对，我指的是每对可能的图案。我还根据我显示的数据提供了一个示例，说明所需的输出应该是什么样子。在显示的数据中，GATA6 和 GATA4 是两个基序。它们同时出现在 3 个不同的 id 中。所以计数是 3。同样，GATA6 和 SRF 在 2 个不同的 id 中同时出现，所以计数是 2。
没有冒犯。我假设您的意思是 adjacent 对，但这在您的帖子中并不清楚。