根据SNP位置和基因起始/结束坐标从另一个数据中指定基因名称

Question

我有两个数据:一个是SNPs及其位置的列表，另一个是基因列表及其起始坐标和结束坐标。使用dplyr，我想添加一个列到SNP数据框中，该列具有每个SNP所在的基因的名称(即SNP的位置位于相同的染色体上，位于该基因的起始坐标/结束坐标之间)。

如果一个SNP不属于任何基因坐标，它应该在基因列中得到"NA“。SNP与基因之间的染色体数目必须匹配。例如，即使第二个SNP的位置在Gene4的起始/结束坐标内，但这并不匹配，因为它们位于不同的染色体上。

SNP数据文件：

CHR  POS  REF  ALT
01   5    C    T
01   10   G    A
02   5    G    T
02   15   C    A
02   20   T    C
03   10   A    G
03   20   C    T

基因数据：

CHR  START  END  GENE_NAME
01   2      8    Gene1
01   12     20   Gene2
01   25     30   Gene3
02   10     18   Gene4
02   25     35   Gene5
03   5      15   Gene6

期望产出：

CHR  POS  REF  ALT  GENE_NAME
01   5    C    T    Gene1
01   10   G    A    NA
02   5    G    T    NA
02   15   C    A    Gene4
02   20   T    C    NA
03   10   A    G    Gene6
03   20   C    T    NA

同样，我想使用dplyr来完成这个任务。提前感谢您的帮助！

Answer 1

使用来自purrr的purrr的一个选项是过滤基于POS的GENE数据和SNP中的CHR，并选择相应的GENE_NAME。

library(dplyr)
library(purrr)

SNP %>%
   mutate(GENE_NAME = map2_chr(POS, CHR, function(x, y) {
       inds = x >= GENE$START & x <= GENE$END & y == GENE$CHR
      if (any(inds)) GENE$GENE_NAME[which.max(inds)] else NA
 }))

#  CHR POS REF ALT GENE_NAME
#1   1   5   C   T     Gene1
#2   1  10   G   A      <NA>
#3   2   5   G   T      <NA>
#4   2  15   C   A     Gene4
#5   2  20   T   C      <NA>
#6   3  10   A   G     Gene6
#7   3  20   C   T      <NA>

在基本R中，这可以使用mapply编写

mapply(function(x, y) {
   inds = x >= GENE$START & x <= GENE$END & y == GENE$CHR
   if (any(inds)) GENE$GENE_NAME[which.max(inds)] else NA
}, SNP$POS, SNP$CHR)

#[1] "Gene1" NA      NA      "Gene4" NA      "Gene6" NA    

数据

SNP <- structure(list(CHR = c(1L, 1L, 2L, 2L, 2L, 3L, 3L), POS = c(5L, 
10L, 5L, 15L, 20L, 10L, 20L), REF = c("C", "G", "G", "C", "T", 
"A", "C"), ALT = c("T", "A", "T", "A", "C", "G", "T")), class = 
"data.frame", row.names = c(NA, -7L))

GENE <- structure(list(CHR = c(1L, 1L, 1L, 2L, 2L, 3L), START = c(2L, 
12L, 25L, 10L, 25L, 5L), END = c(8L, 20L, 30L, 18L, 35L, 15L), 
GENE_NAME = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5", 
"Gene6")), class = "data.frame", row.names = c(NA, -6L))

Answer 2

这是dplyr的一种方法。您只需在START和END的基础上扩展left_join数据，然后用snp扩展left_join-

snp %>% 
  left_join(
    gene %>% 
      group_by(CHR, START, END) %>% 
      mutate(
        POS = list(seq(START, END))
      ) %>% 
      unnest(),
    by = c("CHR", "POS")
  ) %>% 
  select(CHR, POS, REF, ALT, GENE_NAME)

  CHR POS REF ALT GENE_NAME
1   1   5   C   T     Gene1
2   1  10   G   A      <NA>
3   2   5   G   T      <NA>
4   2  15   C   A     Gene4
5   2  20   T   C      <NA>
6   3  10   A   G     Gene6
7   3  20   C   T      <NA>

Answer 3

下面是使用non-equi连接和data.table的一个选项

library(data.table)
setDT(snp)[gene, GENE_NAME := GENE_NAME, on = .(CHR, POS >= START, POS <= END)]
snp
#   CHR POS REF ALT GENE_NAME
#1:   1   5   C   T     Gene1
#2:   1  10   G   A      <NA>
#3:   2   5   G   T      <NA>
#4:   2  15   C   A     Gene4
#5:   2  20   T   C      <NA>
#6:   3  10   A   G     Gene6
#7:   3  20   C   T      <NA>

或使用fuzzyjoin

library(fuzzyjoin)
library(dplyr)
fuzzy_left_join(snp, gene, by = c("CHR", "POS" = "START",
   "POS" = "END"), match_fun = list(`==`, `>=`, `<=`)) %>% 
   select(CHR = CHR.x, POS, REF, ALT, GENE_NAME)
#  CHR POS REF ALT GENE_NAME
#1   1   5   C   T     Gene1
#2   1  10   G   A      <NA>
#3   2   5   G   T      <NA>
#4   2  15   C   A     Gene4
#5   2  20   T   C      <NA>
#6   3  10   A   G     Gene6
#7   3  20   C   T      <NA>

或使用rap选项

library(rap)
snp %>% 
   rap(GENE_NAME = ~ filter(gene, CHR ==  !!CHR, START <= POS, END >= POS) %>% 
          pull(GENE_NAME)) %>% 
   mutate(GENE_NAME = replace(GENE_NAME, !lengths(GENE_NAME), NA)) %>% 
   unnest
#  CHR POS REF ALT GENE_NAME
#1   1   5   C   T     Gene1
#2   1  10   G   A      <NA>
#3   2   5   G   T      <NA>
#4   2  15   C   A     Gene4
#5   2  20   T   C      <NA>
#6   3  10   A   G     Gene6
#7   3  20   C   T      <NA>

数据

snp <- structure(list(CHR = c(1L, 1L, 2L, 2L, 2L, 3L, 3L), POS = c(5L, 
10L, 5L, 15L, 20L, 10L, 20L), REF = c("C", "G", "G", "C", "T", 
"A", "C"), ALT = c("T", "A", "T", "A", "C", "G", "T")),
  class = "data.frame", row.names = c(NA, 
-7L))

gene <- structure(list(CHR = c(1L, 1L, 1L, 2L, 2L, 3L), START = c(2L, 
12L, 25L, 10L, 25L, 5L), END = c(8L, 20L, 30L, 18L, 35L, 15L), 
    GENE_NAME = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5", 
    "Gene6")), class = "data.frame", row.names = c(NA, -6L))