如何按r中其他列中的条件对行进行排序？

Question

我有一个遗传数据集，我希望在其中对样本/基因进行排序，并按基因组中彼此相距一定距离的样本/基因进行分组。

例如，我的数据集如下所示：

#dt1
Gene    chromosome position      CP 
Gene1      1       70000200   1:70000200
Gene2      5       10000476   5:10000476
Gene3      1       70000201   1:70000201
Gene4      5       10000475   5:10000475

我还有一个原点位置数据集：

#dt2
chromosome   position  CP
    1        70005000  1:70005000
    5        10005000  5:10005000

如果基因在我的第二个 dt2 数据集中的任何位置的 +/- 500000 距离内并且在同一条染色体上，我正在尝试对我的第一个数据集中的基因进行分组。我的实际数据中存在一个问题，对于针对多个起源 dt2 位置的基因来说，这可能是正确的，所以我也试图对它最接近的一个进行排序。

输出旨在给出有序组：

Gene   chromosome  position   Group 
Gene1      1       70000200    1
Gene3      1       70000201    1
Gene4      5       10000475    2
Gene2      5       10000476    2

Gene1 和 Gene3 位于原始 dt2 位置的 500000 范围内，并且都在同一条染色体上，因此分组在一起，并且基因 4 和 2 相同

目前我正在尝试这样做：

dt2[, c("low", "high") := .(position - 500000, position  + 500000)]

#find matches on chromosome, with position between low&high
dt1[ dt2, match := i.CP,
     on = .(chromosome, position > low, position < high ) ]

#outputs:
    Gene chromosome position    CP        match
1   Gene1   1   70000200    1:70000200  1:70005000
2   Gene2   5   10000476    5:10000476  5:10005000
3   Gene3   1   70000201    1:70000201  1:70005000
4   Gene4   5   10000475    5:10000475  5:10005000

我在 2 个级别上遇到了这个问题，似乎没有在我的实际数据的匹配列中获得这个输出，所以我想知道是否有其他方法可以尝试对此进行编码。我也在努力将匹配列转换为分组匹配并在我的预期输出中识别我想要的组，我有生物学背景，所以我不确定如何改变这一点 - 任何帮助将不胜感激。

输入数据：

#dt1:
structure(list(Gene = c("Gene1", "Gene2", "Gene3", "Gene4"), 
    chromosome = c(1L, 5L, 1L, 5L), position = c(70000200L, 10000476L, 
    70000201L, 10000475L), CP = c("1:70000200", "5:10000476", 
    "1:70000201", "5:10000475"), match = c("1:70005000", "5:10005000", 
    "1:70005000", "5:10005000")), row.names = c(NA, -4L), class = c("data.table", 
"data.frame"))

#dt2:
structure(list(chromosome = c(1L, 5L), position = c(70005000L, 
10005000L), CP = c("1:70005000", "5:10005000"), low = c(69505000, 
9505000), high = c(70505000, 10505000)), row.names = c(NA, -2L
), class = c("data.table", "data.frame"))

Answer 1

在 Bioconductor 中有专门用于执行此操作的软件包。所以你可以使用的一件事是distanceToNearest()来自GenomicRanges. 首先我们将它们转换为 GRanges 对象：

library(GenomicRanges)
gr1=makeGRangesFromDataFrame(dt1,seqnames.field="chromosome",start.field="position",end.field="position")
values(gr1) = dt1[,c("Gene","CP")]

给 dt2 分组：

dt2$Group = 1:nrow(dt2)
gr2=makeGRangesFromDataFrame(dt2,seqnames.field="chromosome",start.field="position",end.field="position")

此步骤将 gr1 (dt1 GRanges) 中的每一行与 gr2 (dt2) 中最近的范围匹配：

matches = distanceToNearest(gr1,gr2)
Hits object with 4 hits and 1 metadata column:
      queryHits subjectHits |  distance
      <integer>   <integer> | <integer>
  [1]         1           1 |      4799
  [2]         2           2 |      4523
  [3]         3           1 |      4798
  [4]         4           2 |      4524
  -------
  queryLength: 4 / subjectLength: 2

我们将此结果分配回：

dt1$group = NA
dt1$group[queryHits(matches)] = dt2$Group[subjectHits(matches)]
dt1$distance = NA
dt1$distance[queryHits(matches)] = values(matches)$distance[subjectHits(matches)]
dt1

    Gene chromosome position         CP      match group distance
1: Gene1          1 70000200 1:70000200 1:70005000     1     4799
2: Gene2          5 10000476 5:10000476 5:10005000     2     4523
3: Gene3          1 70000201 1:70000201 1:70005000     1     4799
4: Gene4          5 10000475 5:10000475 5:10005000     2     4523

现在你可以过滤掉那些 > 500000

Answer 2

我认为您正在寻找执行非等连接然后通过引用更新的惯用方式，即

dt2[, c("rn", "low", "high") := .(.I, position - 500000L, position  + 500000L)]
#note that you perform the non-equi join first, and then 
#extract the result column before `:=`, which updates by reference
dt1[, Group := dt2[.SD, on=.(chromosome, low<position, high>position), rn]]
dt1

编辑关于多个匹配。在这种情况下，您将需要一个左连接：

dt2[, group := .I][dt1, on=.(chromosome, low<position, high>position)]

Answer 3

如果我理解正确，OP想要

• dt1如果基因位于dt2同一染色体中任何位置的 +/- 500000 距离内，则将它们分组。
• 如果一个基因匹配多个起源dt2位置，则选择最接近的一个。

因此，滚动连接到最近的并随后进行过滤可能是答案。

library(data.table)
dt2[, Group := .I][
  dt1, on = .(chromosome, position), roll = "nearest",
  .(Gene, chromosome, position, CP = i.CP, 
    Group = fifelse(abs(x.position - i.position) <= 500000L, Group, NA_integer_))][
      order(Group, CP)]

     Gene chromosome position         CP Group
1:  Gene1          1 70000200 1:70000200     1
2:  Gene3          1 70000201 1:70000201     1
3: Gene12          1 70005199 1:70005199     1
4: Gene13          1 70005900 1:70005900     2
5:  Gene4          5 10000475 5:10000475     3
6:  Gene2          5 10000476 5:10000476     3
7: Gene11          1 80000200 1:80000200    NA

请注意，此处使用扩展数据集是为了检查该方法是否适用于不同的用例。

Gene11在扩展的数据集中有一个基因 ，它没有被分配到任何组，即Group == NA，因为它在+/- 500000dt2 的距离阈值内没有匹配的位置。

该方法类似于StupidWolf 的GenomicRanges答案，但仅使用data.table功能。

数据

两个数据集都已扩展以涵盖不同的用例：

library(data.table)
dt1 <- fread("Gene    chromosome position      CP 
Gene1      1       70000200   1:70000200
Gene2      5       10000476   5:10000476
Gene3      1       70000201   1:70000201
Gene4      5       10000475   5:10000475
Gene11     1       80000200   1:80000200
Gene12     1       70005199   1:70005199
Gene13     1       70005900   1:70005900")
dt2 <- fread("chromosome   position  CP
    1        70005000  1:70005000
    1        70006000  1:70006000
    5        10005000  5:10005000")